Pseudo-bending method for three-dimensional seismic ray tracing in a spherical earth with discontinuities

We present a 'pseudo-bending' approach to 3-D ray tracing in a spherical earth with discontinuities. This method is based on a three-point perturbation associated with a first-order approximation, while Snell's law in curvilinear coordinates is applied at the discontinuities. We demonstrate the computational accuracy and efficiency of the pseudo-bending method in tracing rays for various velocity models by comparing results with analytical solutions and with results from the bending method. The improvement of efficiency is significant, but is reduced as the number of discontinuities increases. Since the bending approach may be computationally unstable in some situations, even though it is exact, the pseudo-bending approach is preferable for automatic calculation of rays.     

Identification of symmetry planes in weakly anisotropic elastic media

A procedure is proposed to obtain symmetry properties of weakly anisotropic (WA) elastic media by giving 15 WA parameters of qP wave in an arbitrary Cartesian coordinate system. The 15 WA parameters, which linearly depend on 21 elastic elements, form a complete set to determine the symmetry planes in WA materials. The procedure is based on the eigenvalue problems of two matrices. One of the matrices consists of the Voigt and dilatational matrices, and the other is an acoustic tensor defined by an irreducible, completely symmetric and traceless, fourth-rank tensor resulting from decomposition of the fourth-rank elastic tensor. If the eigenvectors are taken as the axes of a new coordinate system (called symmetry Cartesian coordinate system), the transformation of WA parameters from an arbitrary Cartesian coordinate system to a symmetry Cartesian coordinate system can reduce the number of distinct WA parameters of elastic materials except in triclinic medium.     

Computation of qP-wave rays, traveltimes and slowness vectors in orthorhombic media

The eikonal equation is the equation of the phase slowness surface for isotropic and anisotropic media. In general anisotropic media, there is no simple explicit expression for the phase slowness surface. An approximate expression of the eikonal equation may be obtained in weakly anisotropic media. In orthorhombic media, the approximate eikonal equation of the qP wave is the sum of an ellipsoidal form and a more complicated term. The ellipsoidal form corresponds to what we call ellipsoidal anisotropy. Ray equations written in the Hamiltonian formulation are characteristics of the eikonal equation. Ray perturbation theory may be used to compute changes in ray paths and physical attributes (traveltime, polarization, amplitude) due to changes in the medium with respect to a reference medium. Examples obtained in homogeneous orthorhombic media show that a reference medium with ellipsoidal anisotropy is a better choice to develop the perturbation approach than an isotropic reference medium. Models with strong anisotropy can be considered. The comparison with results obtained by an exact ray program shows a relative traveltime error of less than 0.5 per cent for a model with relatively strong anisotropy. We propose a finite element approach in which the medium is divided into a set of elements with polynomial elastic parameter distributions. Inside each element, using a perturbation approach, analytical expressions for rays and traveltimes are obtained Ray tracing reduces to connecting these analytical solutions at the vertices of the cells.     

Dynamic ray tracing on Lagrangian manifolds

Summary. We show that Maslov's extension of the WKBJ method allows an extension of the dynamic ray tracing to wavefields involving caustics of arbitrary form. If the receiver lies off the caustics, then the synthetic seismogram can be obtained by integrating the DRT system along a single ray joining the receiver to the source which may touch caustics. If the receiver-lies in the vicinity of a caustic then DRT has to be carried out along a bunch of rays covering a neighbourhood of the receiver. Our approach encompasses pre-stressed and/or anisotropic media. Initial boundary conditions for a point source embedded in an anisotropic elastic medium are also presented.     

A matrix representation of the potential second-rank gradient tensor for local modelling

Summary. This paper introduces a new matrix computational approach to the local determination of gravity gradients, convenient for comparing with gradient signals from moving base gradiometer systems or calculating topographic effects at instrument heights. The method represents a practical alternative to the more conventional spherical harmonics formulation, primarily global in nature, and it may be considered as an extension to other previously used local representations, such as point masses. Important characteristics of the analytical development outlined herein are its conceptual simplicity and the possibility of obtaining at once, up to a certain order n , and in an arbitrary Cartesian coordinate system, the symmetric point gradient tensor of second rank.     

Invariant expressions for elastic radiation from a fault with prescribed slip

Summary New formulae are obtained for the displacement potentials and displacements due to a point source with moment tensor, and for a fault with prescribed slip. These formulae, unlike previous formulae, are invariant, i.e. they are valid in any coordinate system, not just Cartesian coordinates or orthogonal curvilinear coordinates.
Apart from their invariance, these formulae have other advantages: they are exceedingly simple; the expressions for P -motion are nearly the same as those for S -motion; and the separation of far field motion from final static displacement is automatic.     

Gaussian beams for surface waves in laterally slowly-varying media

Summary. Asymptotic ray theory is applied to surface waves in a medium where the lateral variations of structure are very smooth. Using ray-centred coordinates, parabolic equations are obtained for lateral variations while vertical structural variations at a given point are specified by eigenfunctions of normal mode theory as for the laterally homogeneous case. Final results on wavefields close to a ray can be expressed by formulations similar to those for elastic body waves in 2-D laterally heterogeneous media, except that the vertical dependence is described by eigenfunctions of 'local' Love or Rayleigh waves. The transport equation is written in terms of geometrical-ray spreading, group velocity and an energy integral. For the horizontal components there are both principal and additional components to describe the curvature of rays along the surface, as in the case of elastic body waves. The vertical component is decoupled from the horizontal components. With complex parameters the solutions for the dynamic ray tracing system correspond to Gaussian beams: the amplitude distribution is bell-shaped along the direction perpendicular to the ray and the solution is regular everywhere, even at caustics. Most of the characteristics of Gaussian beams for 2-D elastic body waves are also applicable to the surface wave case. At each frequency the solution may be regarded as a set of eigenfunctions propagating over a 2-D surface according to the phase velocity mapping.     

Ray amplitudes of seismic body waves in laterally inhomogeneous media

Summary . The most complicated part in the computation of ray amplitudes of seismic body waves in laterally inhomogeneous media with curved interfaces lies in the evaluation of the geometrical spreading. Geometrical spreading can be simply expressed in terms of the Jacobian J of the transformation from the Cartesian into ray coordinates. Several systems of ordinary differential equations to compute the function J are suggested. For general three-dimensional media, in which the velocity changes with all the three spatial coordinates, a system of three non-linear ordinary differential equations of the first order is derived. If the velocity does not depend on one coordinate, the system of equations reduces to only one non-linear differential equation. The initial conditions for these differential equations at point (or line) source and at points of intersection of the ray with curved interfaces are presented.     

First-order P-wave ray synthetic seismograms in inhomogeneous weakly anisotropic media

We propose approximate equations for P -wave ray theory Green's function for smooth inhomogeneous weakly anisotropic media. Equations are based on perturbation theory, in which deviations of anisotropy from isotropy are considered to be the first-order quantities. For evaluation of the approximate Green's function, earlier derived first-order ray tracing equations and in this paper derived first-order dynamic ray tracing equations are used.
The first-order ray theory P -wave Green's function for inhomogeneous, weakly anisotropic media of arbitrary symmetry depends, at most, on 15 weak-anisotropy parameters. For anisotropic media of higher-symmetry than monoclinic, all equations involved differ only slightly from the corresponding equations for isotropic media. For vanishing anisotropy, the equations reduce to equations for computation of standard ray theory Green's function for isotropic media. These properties make the proposed approximate Green's function an easy and natural substitute of traditional Green's function for isotropic media.
Numerical tests for configuration and models used in seismic prospecting indicate negligible dependence of accuracy of the approximate Green's function on inhomogeneity of the medium. Accuracy depends more strongly on strength of anisotropy in general and on angular variation of phase velocity due to anisotropy in particular. For example, for anisotropy of about 8 per cent, considered in the examples presented, the relative errors of the geometrical spreading are usually under 1 per cent; for anisotropy of about 20 per cent, however, they may locally reach as much as 20 per cent.     

A simple notation for differential vector expressions in orthogonal curvilinear coordinates

We present a simple notation for performing differential vector operations in orthogonal curvilinear coordinates, and for easily obtaining partial differential expressions in terms of the physical components. We express n th-order tensors as the summed products of the physical components and n th-order polyads of unit vectors (an extension of Gibbs dyadic notation convenient for a summation convention). By defining a gradient operator with partial derivatives balanced by the inverse scale factors, differential vector (or tensor) operations in orthogonal coordinates do not require the covariant/contravariant notation. Our primary focus is on spherical-polar coordinates, but we also derive formulae which may be applied to arbitrary orthogonal coordinate systems. The simpler case of cylindrical-polar coordinates is briefly discussed. We also offer a compact form for the gradient and divergence of general second-order tensors in orthogonal curvilinear coordinates, which are generally unavailable in standard handbooks. We show how our notation relates to that of tensor analysis/differential geometry. As the analysis is not restricted to Euclidean geometry, our notation may be extended to Riemannian surfaces, such as spherical surfaces, so long as an orthogonal coordinate system is utilized. We discuss the Navier-Stokes equation for the case of spatially variable viscosity coefficients.     

Transient SH-waves in a wedge-shaped layer

Summary. A generalized ray theory for transient SH -waves in a wedge-shaped layer over an elastic half-space is developed in this paper. The ray integrals for multiply reflected waves in the layer are derived in terms of two systems of coordinates and two sets of local wavenumbers, one along the free surface and the other along the sloped interface. All local wavenumbers are then transformed to a common wavenumber in all ray integrals which are evaluated by the Cagniard method. Results for the first motion approximation are in agreement with previous investigations.     

Synthetic seismograms for complex three-dimensional geometries using an analytical–numerical algorithm

Summary. An algorithm which is part analytical and part numerical is suggested for the computation of complete synthetic seismograms for complex three-dimensional geological structures with radial symmetry. A partial separation of variables based on the combination of a finite Fourier integral transform with respect to the spatial coordinate z together with the finite difference method is the essence of the algorithm. Upon application of the finite transform the problem reduces to solving a system of equations containing only partial derivatives with respect to one spatial coordinate ( r ) and time. As radial symmetry is assumed, there is no functional dependence on φ in the cylindrical system of coordinates ( r , φ, z ). The coefficients of the transformed equations may contain finite Fourier integrals of the z dependence of the elastic parameters. Several examples of synthetic seismograms computed for both SH - and P – SV -waves propagating in complex subsurface geometries are presented and their interpretation discussed.     

Computation of high-frequency seismic wavefields in 3-D laterally inhomogeneous anisotropic media

Summary. An algorithm for the computation of travel times, ray amplitudes and ray synthetic seismograms in 3-D laterally inhomogeneous media composed of isotropic and anisotropic layers is described. All 21 independent elastic parameters may vary within the anisotropic layers. Rays and travel times are evaluated by numerical solution of the ray tracing equations. Ray amplitudes are determined by evaluating reflection/ transmission coefficients and the geometrical spreading along individual rays. The geometrical spreading is computed approximately by numerical measurement of the cross-sectional area of the ray tube formed by three neighbouring rays. A similar approximate procedure is used for the determination of the coefficients of the paraxial ray approximation. The ray paraxial approximation makes computation of synthetic seismograms on the surface of the model very efficient. Examples of ray synthetic seismograms computed with a program package based on the described algorithm are presented.     

Weak-contrast reflection/transmission coefficients in weakly anisotropic elastic media: P-wave incidence

We present approximate displacement and energy PP and PS reflection/transmission coefficients for weak-contrast interfaces in general weakly anisotropic elastic media. The coefficients were obtained by applying first-order perturbation theory and then expressed in a compact and relatively simple form. The formulae can be used for arbitrary orientations of the incidence plane and interface, without the need to transform the elasticity parameters to a local Cartesian coordinate system. The accuracy of the approximate formulae is illustrated for the PS reflection coefficient for two synthetic models. For these models, we also study the possibility of using the approximate PP reflection coefficient in the inverse problem.     

Surface wave tomography for azimuthal anisotropy in a strongly reduced parameter space

Large scale seismic anisotropy in the Earth's mantle is likely dynamically supported by the mantle's deformation; therefore, tomographic imaging of 3-D anisotropic mantle seismic velocity structure is an important tool to understand the dynamics of the mantle. While many previous studies have focused on special cases of symmetry of the elastic properties, it would be desirable for evaluation of dynamic models to allow more general axis orientation. In this study, we derive 3-D finite-frequency surface wave sensitivity kernels based on the Born approximation using a general expression for a hexagonal medium with an arbitrarily oriented symmetry axis. This results in kernels for two isotropic elastic coefficients, three coefficients that define the strength of anisotropy, and two angles that define the symmetry axis. The particular parametrization is chosen to allow for a physically meaningful method for reducing the number of parameters considered in an inversion, while allowing for straightforward integration with existing approaches for modelling body wave splitting intensity measurements. Example kernels calculated with this method reveal physical interpretations of how surface waveforms are affected by 3-D velocity perturbations, while also demonstrating the non-linearity of the problem as a function of symmetry axis orientation. The expressions are numerically validated using the spectral element method. While challenges remain in determining the best inversion scheme to appropriately handle the non-linearity, the approach derived here has great promise in allowing large scale models with resolution of both the strength and orientation of anisotropy.     

About the influence of pre-stress upon adiabatic perturbations of the Earth

Summary. In this paper we examine the influence of the state of stress in the equilibrium configuration of the Earth (i.e. the pre-stress) upon its adiabatic perturbations. The equations governing these perturbations to the first order (Woodhouse & Dahlen; Dahlen) are re-derived using a Lagrangian approach. Different expressions of the sesquilinear form associated to the elastic-gravitational operator are given. One of these provides a way to extend to hydrostatically pre-stressed solids the criterion of local stability given by Friedman & Schutz for uniformly rotating fluids. Then the propagation in the Earth of seismic wavefronts is considered. It is shown that the nature of these different wavefronts is entirely determined by the quadratic coefficients of the development of the specific internal energy variation, corresponding to isentropic evolution, with respect to the Lagrangian finite deformation tensor. Expressions for the velocities of the various waves are given as functions of incidence angle and pre-stress for orthotropic elastic material. In the particular case where the elastic parameters depend only on one coordinate of a curvilinear system and the axis of orthotropy of the material coincides with the corresponding natural base vector, the elastodynamic equations are reduced to a simple system for a displacement stress vector, using surface operators. In particular for spherical geometry, equations are obtained which generalize to orthotropic pre-stress those given by Alterman et al. and Takeuchi & Saito.     

Existence of transverse waves in anisotropic poroelastic media

Out of the four waves in an anisotropic poroelastic medium, two are termed as quasi-transverse waves. The prefix 'quasi' refers to their polarizations being nearly, but not exactly, perpendicular to direction of propagation. In this composite medium, unlike perfectly elastic medium, the propagation of a longitudinal wave along a phase direction may not be accompanied by transverse waves. The existence of a transverse wave in anisotropic poroelastic media is ensured by the two equations restricting the choice of elastic coefficients of porous aggregate as well as fluid–solid coupling. Necessary and sufficient conditions for the existence of transverse waves along the coordinate axes and in the coordinate planes for general anisotropy are discussed. The discussion is extended to the case of orthotropic materials and existence for few specific phase directions is also explored. The conditions for the transverse waves decided on the basis of their apparent polarizations, that is, particle motion being perpendicular to ray direction, are also discussed. For a particular numerical model, the existence of these apparent transverse waves is solved numerically for phase directions in coordinate planes. For general directions of phase propagation, the existence of these transverse waves is checked graphically for the chosen numerical model.     

Seismic body waves in anisotropic media: reflection and refraction at a plane interface

Summary. A formulation is derived for calculating the energy division among waves generated by plane waves incident on a boundary between generally anisotropic media. A comprehensive account is presented for P, SV and SH waves incident from an isotropic half-space on an orthorhombic olivine half-space, where the interface is parallel to a plane of elastic symmetry. For comparison, a less anisotropic medium having transverse isotropy with a horizontal axis of symmetry is also considered. The particle motion polarizations of waves in anisotropic medium differ greatly from the polarizations in isotropic media, and are an important diagnostic of the presence of anisotropy. Incident P and SV waves generate quasi- SH waves, and incident SH waves generate quasi- P and quasi- SV waves, often of considerable relative magnitude. The direction of energy transport diverges from the propagation direction.     

Characterization of crack distribution: fabric analysis versus ultrasonic inversion

We analyse the relation between rock fabric, expressed by the preferred orientation of rock-forming minerals and microcracks, and elastic anisotropy of crystalline rock from the KTB pilot well. Detailed analyses of mineralogical composition, textures and microcrack fabrics were performed. In addition, ultrasonic velocity measurements of spherical samples in several directions were carried out at various confining pressures, and inverted in terms of the complete set of 21 elastic constants. By comparing the elastic tensors of the rocks at the final confining pressure (at which most of the microcracks are closed) with those at a lower pressure level, it is possible to separate the anisotropy induced by microcracks from that caused by mineral alignment. In contrast to previous work, no a priori knowledge of the type of anisotropy (triclinic, monoclinic, orthotropic etc.), or of the spatial orientation of the symmetry elements (planes, axes) of the cracked rock or of the intact rock is assumed. Furthermore, no restrictive assumptions on the orientation distribution function and the shape of the cracks are needed.
The results show that the elastic anisotropy characteristics, whether they are related to the microcracks or to the rock-forming minerals, are clearly correlated with the directly observed rock fabrics. We show that the symmetry directions of the mineral fabric and of microcrack fabric agree. A further result is that the microcrack-induced anisotropy dominates the other causes of anisotropy at confining pressures smaller than a few tens of megapascals, the situation being reversed at higher pressures. The laboratory data are quantitatively compared with sonic log data from the KTB well, showing the influence of pore fluids, effective pressure and crack density reduction on the anisotropy in situ .     

Finite-frequency vectorial tomography: a new method for high-resolution imaging of upper mantle anisotropy

Amplitude measurements of the transverse component of SKS waves, the so-called splitting intensity, can be used to formulate a non-linear inverse problem to image the 3-D variations of upper mantle anisotropy. Assuming transverse isotropy (or hexagonal symmetry), one can parametrize anisotropy by two anisotropic parameters and two angles describing the orientation of the symmetry axis. These can also be written as two collinear pseudo-vectors. The tomographic process consists of retrieving the spatial distribution of these pseudo-vectors, and thus resembles surface wave vectorial tomography. Spatial resolution results from the sensitivity of low-frequency SKS waves to seismic anisotropy off the ray path. The expressions for the 3-D sensitivity kernels for splitting intensity are derived, including the near-field contributions, and validated by comparison with a full wave equation solution based upon the finite element method. These sensitivity kernels are valid for any orientation of the symmetry axis, and thus generalize previous results that were only valid for a horizontal symmetry axis. It is shown that both lateral and vertical subwavelength variations of anisotropy can be retrieved with a dense array of broad-band stations, even in the case of vertically propagating SKS waves.