Paraxial ray approximations in the computation of seismic wavefields in inhomogeneous media

Summary. Several important applications of the paraxial ray approximation (PRA) to numerical modelling of high-frequency seismic body wavefields are discussed. The PRA can be used to evaluate the displacement vector not only directly on the ray, as in the standard ray method. but also approximately in the vicinity of this ray. The PRA also offers simple ways of approximate evaluation of paraxial rays, situated in the vicinity of the central ray, and of two-point ray tracing. A very important application of the PRA consists in a simple, fast and effective Computation of body-wave synthetic seismograms in general, 3-D, laterally inhomogeneous, layered structures. Examples of synthetic seismograms for 3-D structures, computed using the PRA, are presented.     

Investigation of 3D wavefields of reflected shear-waves and converted waves: mathematical modelling and reflection data processing

Summary. Seismic investigations using shear-wave and converted wave techniques show that very often reflected PS - and SS -waves have anomalous polarizations ( accessory components ). This phenomenon cannot be explained in terms of isotropic models with dipping boundaries. Computations of synthetic seismograms of reflected PS - and SS -waves were made for different models of transversely isotropic media with dipping anisotropic symmetry axes not normal to the boundaries. Synthetic seismograms were computed by ray techniques using an optimization algorithm to construct all rays arriving at a given receiver. These computations indicate that accessory components arise when the medium above the boundary is anisotropic, where they are caused by the constructive interference of qSV - and qSH -waves. If a low-velocity layer is present, displacement vectors of both waves have horizontal projections which are approximately orthogonal. The algorithm for wave separation is presented and some results of its use are given.     

A method for the computation of finite frequency body wave synthetic seismograms in laterally varying media

Summary. Body wave synthetic siesmograms for laterally varying media are computed by means of a slowness implementation of the extended WKBJ (EWKBJ) theory of Frazer & Phinney. An EWKBJ seismogram is computed by first tracing rays through a particular model to obtain conventional ray information (travel time, ray end point, ray slowness) and then using these data in the finite frequency integral expression for the EWKBJ seismogram. The EWKBJ seismograms compare favourably to geometrical ray theory (GRT) seismograms but are significantly better because of the finite frequency nature of the EWKBJ calculation. More realistic behaviour is obtained with EWKBJ seismograms at normal seismic frequencies near caustics, where the GRT amplitude is infinite, and within geometrical shadow zones where GRT predicts zero amplitudes. In addition the EWKBJ calculation is more sensitive than GRT to focuses and defocuses in the ray field. The major disadvantage of the EWKBJ calculation is the additional computer time over that of GRT, necessary to calculate one seismogram although an EWKBJ seismogram costs much less to compute than a reflectivity seismogram. Another disadvantage of EWKBJ theory is the generation of spurious, non-geometrical phases that are associated with rapidly varying lateral inhomogeneities. Fortunately the amplitudes of these spurious phases are usually much lower than that of neighbouring geometrical phases so that the spurious phases can usually be ignored. When this observation is combined with the moderately increased computational time of the EWKBJ calculation then the gain in finite frequency character significantly outweighs any disadvantages.     

Three-dimensional geometrical ray theory and modelling of transmitted seismic energy of data from the Nevada Test Site

This paper presents a geometrically based algorithm for computing synthetic seismograms for energy transmitted through a 3-D velocity distribution. 3-D ray tracing is performed to compute the traveltimes and geometrical spreading (amplitude). The formulations of both kinematic and dynamic ray-tracing systems are presented. The two-point ray-tracing problem is solved by systematically updating the initial conditions and adjusting the ray direction until the ray intersects the specified endpoint. The amount of adjustment required depends on the derivatives of the position with respect to the given starting angles between consecutive rays. The algorithm uses derivatives to define the steepest-descent direction and to update the initial directions. The convergence rate depends on the complexity of the model.
Test seismograms compare favourably with those from a 2-D asymptotic ray theory algorithm and a 3-D Gaussian-beam algorithm. The algorithm is flexible in modelling arbitrary source and recorder geometries for various smoothly varying 3-D velocity distributions. The algorithm is further tested by simulating surface-to-tunnel vibroseis field data. Shear waves as well as compressional waves may be approximately included. Application of the algorithm to a data set from the Rainier Mesa of the Nevada Test Site produced a good fit to the transmitted (first arrival) traveltimes and amplitudes, with approximately 15 per cent variation in the local 3-D velocity.     

A method for calculating synthetic seismograms in laterally varying media

Summary An effective algorithm for computing synthetic seismograms in laterally inhomogeneous media has been developed. The method, based on zero-order asymptotic ray theory, is primarily intended for use in refraction and reflection studies and provides an economical means of seismic modelling.
A given smoothed velocity-depth-distance model is divided into small squares with constant seismic parameters and first-order interfaces are represented by an arbitrary number of dipping linear segments. The computation of ray propagation and amplitudes through such a model does not involve complicated analytic expressions and therefore minimizes computer time.
Amplitudes are determined by geometrical spreading of spherical wave-fronts and energy partitioning at interfaces. Synthetic seismograms calculated for laterally homogeneous models are in good agreement with those obtained by the Reflectivity Method.     

Ray-reflectivity method for SH-waves in stacks of thin and thick layers

Summary Reflectivity and ray theories are united to produce a hybrid technique of computing synthetic seismograms for a plane layered medium in subcritical regions. Numerical experiments have indicated that this technique is useful when the depth structure is one composed of thick layers separated by finely layered zones. As the theory for wave propagation in a plane layered medium is well known, the simple SH case is investigated so that the basic idea of the method may be conveyed without an excess of mathematics that would be necessitated if the P-SV problem were considered.
In computing the ray-reflectivity seismogram, the thick layers are treated using asymptotic ray theory while the thin-layered zones are treated as quasiinterfaces where analogues of reflection and transmission coefficients called reflectivities and transmittivities are calculated utilizing a Thomson-Haskell formulation. A stationary phase approximation is employed when evaluating the integral which gives the displacement due to an arbitrary ray propagating in the thick layers of the above-mentioned medium, and the validity of this approximation is discussed.
A comparison of ray, numerical integration (reflectivity) and ray-reflectivity synthetic sections indicates that this method yields quite acceptable results for subcritical reflection work and is suitable for application in seismic interpretation as individual arrivals associated with ray-paths in the thick layers may be identified. Furthermore, the method is quite cost efficient and may be extended to a medium where the thick layers are non-planar using asymptotic ray theory in these layers.     

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.     

Ray tracing and inhomogeneous dynamic ray tracing for anisotropy specified in curvilinear coordinates

Ray tracing has recently been expressed for anisotropy specified in a local Cartesian coordinate system, which may vary continuously in a model specified by elastic parameters. It takes advantage of the fact that anisotropy is often of a simpler nature locally (and is thus specified by a smaller number of elastic parameters) and that the orientation of its symmetry elements may vary. Here we extend this approach by replacing the local Cartesian coordinate system with a curvilinear coordinate system of global extent and by applying the new approach to ray tracing and inhomogeneous dynamic ray tracing. The curvilinear coordinate system is orthogonal and is constructed so that the coordinate axes are consistent with the considered anisotropy of the medium. Our formulation allows for computation of ray attributes (e.g. ray velocity vector and paraxial ray attributes) in the curvilinear coordinate system, while rays are computed in global Cartesian coordinates. Compared to the classic formulation in terms of 21 elastic moduli in global Cartesian coordinates, the main advantages are improved efficiency, lower computer-memory requirements, and conservation of anisotropic symmetry throughout the model.     

Kirchhoff–Helmholtz reflection seismograms in a laterally inhomogeneous multi-layered elastic medium – II. Computations

Summary. High-frequency reflection and refraction seismograms for laterally variable multi-layered elastic media are computed by using the frequency domain elastic Kirchhoff–Helmholtz (KH) theory of Frazer and Sen. Both source and receiver wavefields are expanded in series of generalized rays and then elastic (KH) theory is applied to determine the coupling between each source ray and each receiver ray at each interface. The motion at the receiver is given as a series of integrals, one for each generalized ray. We use geometrical optics and plane wave reflection and transmission coefficients for rapid evaluation of the integrand. When the source or the receiver ray field has caustics on the surface of integration geometrical ray theory breaks down and this gives rise to singularities in the KH integrand. We repair this using methods suggested by Frazer and Sen.
Examples of reflection seismograms for 2-D structures computed by elastic KH theory are shown. Those for a vertical fault scarp structure are compared with the seismograms obtained by physical modelling. Then OBS data obtained from the mid-America trench offshore Guatemala area are analysed by computing KH synthetics for a velocity model that has been proposed for that area. Our analysis indicates the existence of a small low-velocity zone off the trench axis.
No head wave arrivals are obtained in our KH synthetics since we do not consider multiple interactions of a ray with an interface. The nearly discontinuous behaviour of elastic R/T coefficients near the critical angle causes small spurious phases which arrive later than the correct arrivals.     

Synthetic body wave seismograms for laterally varying layered structures by the Gaussian beam method

Summary. Several approaches to computing body wave seismograms in 2–D and 3–D laterally inhomogeneous layered structures are suggested. They are based on the Gaussian beam method, which has been recently applied to the evaluation of time-harmonic high-frequency wavefields in inhomogeneous media. Three variants are discussed in some detail: the spectral method, the convolutory method and the wave-packet method. The most promising seems to be the wave-packet approach. In this approach, the wavefield, generated by a source, is expanded into a system of wave packets, which propagate along rays from the source in all directions. The wave packets change their properties due to diffusion, spreading, reflections/transmissions, etc. The resulting seismogram at any point of the medium is then obtained as a superposition of those packets which propagate close to the point. The final expressions in all the three methods are regular even in regions, in which the ray method fails, e.g. in the vicinity of caustics, in the critical region, at boundaries between shadow and illuminated regions, etc. Moreover, they are not as sensitive to the minor details of the medium as the ray method and, what is more, they remove the time-consuming two-point ray tracing from computations. Numerical examples of synthetic seismograms computed by the wave-packet approach are presented.     

Synthetic body wave seismograms for three-dimensional laterally varying media

Summary. Two methods of computing body wave synthetic seismograms in three-dimensional laterally varying media are discussed. Both these methods are based on the summation of Gaussian beams. In the first, the initial beam parameters are chosen at the source, in the second at the beam endpoints. Both these variants eliminate the ray method singularities. The expansion of the wavefield into plane waves may be considered as the limiting case of the first approach and the Chapman–Maslov method as the limiting case of the second approach. Computer algorithms are briefly described and numerical examples presented. In the first numerical example, the comparisons of the two approaches, based on summing Gaussian beams, with the reflectivity method indicate that the computed synthetic seismograms are satisfactorily accurate even in the caustic region. The next example suggests that the two methods discussed can be simply and effectively applied to 3-D laterally inhomogeneous structures.     

Backscattering from spherical elastic inclusions and accuracy of the Kirchhoff approximation for curved interfaces

The Kirchhoff (or tangent plane) approximation, derived from the theoretically complete Kirchhoff–Helmholtz integral representation for the seismic wavefield, has been used extensively for the analysis of seismic-wave scattering from irregular interfaces; however, the accuracy of this method for curved interfaces has not been rigorously established. This paper describes an efficient Kirchhoff algorithm to simulate scattered waves from an arbitrarily curved interface in an elastic medium. Synthetic seismograms computed using this algorithm are compared with exact synthetics computed using analytical formulae for scattering of plane P waves by a spherical elastic inclusion. A windowing technique is used to remove strong internal reverberations from the analytical solution. Although the Kirchhoff method tends to underestimate the total scattering intensity, the accuracy of the approximation improves with increasing value of the wavenumber-radius product, kR . The arrival times and pulse shapes of primary reflections from the sphere are well approximated using the Kirchhoff approach regardless of curvature of the scattering surface, but the amplitudes are significantly underestimated for kR ≤ 5. The results of this work provide some new guidelines to assess the accuracy of Kirchhoff-synthetic seismograms for curved interfaces.     

Quantized Ray Theory for Non-Zero Focal Depths

The relation between p-Δcurves for surface and deep focus sources is investigated in order to construct synthetic body wave seismograms for non-zero focal depths by the quantized ray theory algorithm. The transformation of a surface focus p-Δ curve into a deep focus p-Δ curve is denned in terms of that curve which corresponds to surface focus rays reflected from the depth at which the deep focus is located. By analogy with the geometry of the surface focus formulation, paths of integration to obtain absolute travel-time and velocity-depth curves can be denned in the p-Δ plane. Explicit inversion from deep focus data is possible only when the velocity-depth structure above the depth of focus is known. Through a comparison of short period quantized ray theory synthetic seismograms with similar Cagniard-de Hoop computations, it is shown that quantized ray theory can be used for accurate predictions of body wave amplitude behaviour corresponding to a wide range of focal depths.     

Analytical computation of coupling coefficients in non-Poissonian media

As indicated by numerical tests, the computation of synthetic seismograms based on the modal summation technique, for 2-D or 3-D media, requires the computation of the transmission and reflection coupling coefficients, at sharp lateral heterogeneities, with a formulation that takes into account the possible non-Poissonian character of the media.     

Wide-band coupling of Earth's normal modes due to anisotropic inner core structure

We investigate the importance of wide-band coupling of normal modes due to inner core anisotropy. We compare four different seismic models of inner core anisotropy, which were obtained by others using the splitting of Earth's normal modes. These models have been developed using a self-coupling (SC) approximation, which assumes that coupling between nearby modes through anisotropic inner core structure is negligible. We test the SC approximation by comparing the frequencies and quality factors of 90 inner core sensitive modes, computed for these models using either the SC approximation or full-coupling (FC) among large groups of modes. We find significant shifts in the quality factors and frequencies for some modes. Groups of modes which significantly couple together are constructed for six target modes. These groups are model dependent and in some cases contain large numbers of modes. Synthetic seismograms are calculated to show that the difference between SC and FC is observable on the scale of seismograms and of the same order of magnitude as the difference between synthetic and observed seismograms. Thus, future models of inner core anisotropy should take cross-coupling between large groups of modes into account.     

P-SH conversions in a flat-layered medium with anisotropy of arbitrary orientation

P-SH conversion is commonly observed in teleseismic P waves, and is often attributed to dipping interfaces beneath the receiver. Our modelling suggests an alternative explanation in terms of flat-layered anisotropy. We use reflectivity techniques to compute three-component synthetic seismograms in a 1-D anisotropic layered medium. For each layer of the medium, we prescribe values of seismic velocities and hexagonally symmetric anisotropy about a common symmetry axis of arbitrary orientation. A compressional wave in an anisotropic velocity structure suffers conversion to both SV -and SH -polarized shear waves, unless the axis of symmetry is everywhere vertical or the wave travels parallel to all symmetry axes. The P-SV conversion forms the basis of the widely used 'receiver function' technique. The P-SH conversion occurs at interfaces where one or both layers are anisotropic. A tilted axis of symmetry and a dipping interface in isotropic media produce similar amplitudes of both direct ( P ) and converted ( Ps ) phases, leaving the backazimuth variation of the P-Ps delay as the main discriminant. Seismic anisotropy with a tilted symmetry axis leads to complex synthetic seismograms in velocity models composed of just a few flat homogeneous layers. It is possible therefore to model observations of P coda with prominent transverse components with relatively simple 1-D velocity structures. Successful retrieval of salient model characteristics appears possible using multiple realizations of a genetic-algorithm (GA) inversion of P coda from several backazimuths. Using GA inversion, we determine that six P coda recorded at station ARU in central Russia are consistent with models that possess strong (> 10 per cent) anisotropy in the top 5 km and between 30 and 43 km depth. The symmetry axes are tilted, and appear aligned with the seismic anisotropy orientation in the mantle under ARU suggested by SKS splitting.     

Quasi-isotropic approximation of ray theory for anisotropic media

Wave propagation in weakly anisotropic inhomogeneous media is studied by the quasi-isotropic approximation of ray theory. The approach is based on the ray-tracing and dynamic ray-tracing differential equations for an isotropic background medium. In addition, it requires the integration of a system of two complex coupled differential equations along the isotropic ray.
The interference of the qS waves is described by traveltime and polarization corrections of interacting isotropic S waves. For qP waves the approach leads to a correction of the traveltime of the P wave in the isotropic background medium.
Seismograms and particle-motion diagrams obtained from numerical computations are presented for models with different strengths of anisotropy.
The equivalence of the quasi-isotropic approximation and the quasi-shear-wave coupling theory is demonstrated. The quasi-isotropic approximation allows for a consideration of the limit from weak anisotropy to isotropy, especially in the case of qS waves, where the usual ray theory for anisotropic media fails.     

A fast converging and anti-aliasing algorithm for Green's functions in terms of spherical or cylindrical harmonics

Geophysical observables are generally related to earth structure and source parameters in a complicated non-linear way. Consequently, a large number of forward modelling processes are commonly necessary to obtain a satisfactory estimate of such parameters from observed data. The most time-consuming part of the forward modelling is the computation of the Green's functions of the different earth models to be tested. In this study, we present a fast converging algorithm: the differential transform method for the computation of Green's functions in terms of spherical or cylindrical harmonics. In this method, a deconvolutable high-pass filter is used to enhance the numerical significance of the far-field spectrum of Green's functions. Compared with existing fast converging algorithms such as the Kummer's transformation and the disc factor method, the differential transform method is more efficient except for the extremely near-source region. The new method can be used to suppress numerical phases (non-physical seismic signals) associated with the aliasing effect that may arise in synthetic seismograms when the latter are computed from a windowed wavenumber (or slowness) spectrum. The numerical efficiency of the new method is demonstrated via two representative tests.     

Time-domain computation of scattering from a heterogeneous layer

Adopting Born and ray approximations, time-domain synthetic seismograms for P-P and P-S scattering from a plane wave incident on a thin, laterally heterogeneous layer are presented in this paper. The time-domain P coda is a convolution between a structure function and the second-order derivative of the time function of the incident P wave. Examples of synthetic seismograms are given using a time function from a computed short-period seismogram for a point explosive source in a half-space. These show that it is impossible, with realistic values of the parameters involved, to generate significant codas when only single scattering is involved.     

Gaussian beams in elastic 2-D laterally varying layered structures

Summary. In a paper by Červený & Pšenčik, high-frequency Gaussian beams in elastic 2-D, laterally inhomogeneous, smooth media were investigated as asymptotic high-frequency solutions of elastodynamic equations, concentrated close to rays of P - and S -waves. This paper generalizes the above results for 2-D, laterally inhomogeneous, layered structures. Gaussian beams concentrated close to any multiply-reflected, possibly converted, ray are investigated. Gaussian beams are regular everywhere, including caustic regions. The paraxial ray approximation, which allows the wavefield in the zero-order ray approximation to be evaluated not only directly on the ray, but also in its vicinity, is derived as a limiting case of the Gaussian beams.