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.     

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.     

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.     

Chaotic ray behaviour in regional seismology

Rays propagating through strongly laterally varying media exhibit chaotic behaviour. This means that initially close rays diverge exponentially, rather than according to a power law. This chaotic behaviour is especially pronounced if the medium contains laterally varying interfaces. By studying simple 2-D and 3-D versions of models with laterally varying interfaces, the importance of chaotic ray behaviour is determined. A model of the Moho below Germany produces sharp variations with epicentral distance of the number of arrivals. In addition, the number of caustics grows dramatically: up to 1200 caustics are present between a distance of 0 and 800 km. Using the theory of Hamiltonian systems, a more in-depth study of the chaotic character of the ray equations is obtained. It is found that for realistic heterogeneous models most of the relevant rays will exhibit chaotic behaviour. The degree of chaos is quantified in terms of predictability horizons. Beyond the predictability horizons ray tracing cannot be carried out accurately. For the models under consideration, the length from the source to the predictability horizon has an order of magnitude of 1000 km. The chaotic behaviour of the rays makes it necessary to use extensions of asymptotic ray theory, such as Maslov theory, to compute seismic waveforms. It is shown that pseudo-caustics, an important obstacle in computing Maslov synthetics, are a generic feature of the 2-D laterally varying models that are studied. Eventually, the use of asymptotic methods is restricted because of the inaccuracy in the computation of the ray paths.     

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.     

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.     

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.     

Finite difference seismograms for laterally varying marine models

Summary. The method of finite differences is applied to the elastic wave equation to generate synthetic seismograms for laterally varying seafloor structures. The results are compared with borehole seismic data from the Gulf of California (Deep Sea Drilling Project Site 485) in which lines were shot over flat and rough topography. The significant new phenomenon observed in both the synthetic seismograms and the field data is the generation of a 'double head wave' due to the interaction of the incident wavefront with the side of a hill and the flat seafoor adjacent to the hill.
In these models the hills are on the order of a seismic wavelength in height and steep velocity gradients occur over distances comparable to wavelengths. Ray theoretical methods would not be suitable for studying such structures. True amplitude record sections are obtained by the finite difference method, which show for these models that the head wave generated at the flat seafloor adjacent to the hill is lower in amplitude than if the hill were not present and is lower in amplitude than the head wave generated at the hill.
A second feature which is important for borehole receivers is the existence of the 'direct wave root' in the upper basement. This energy occurs below the sharp interface when the direct wave impinges on the interface from above. There is no corresponding Snell's law ray path for this energy and the energy is evanescent with depth in the lower medium.
The properties of both the double head wave and the direct wave root are clearly demonstrated in the finite difference 'snapshot' displays.     

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.     

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.     

Efficient calculation of Gaussian-beam seismograms for two-dimensional inhomogeneous media

Summary. This paper discusses several aspects of the calculation of theoretical seismograms for two-dimensional inhomogeneous media with the method of Gaussian beams. The most important steps of this method, kinematic and dynamic ray tracing, can be performed very efficiently, if the model cross-section is subdivided into triangles with linear velocity laws. Each Gaussian beam is characterized by a complex beam constant ε which determines its width and phase-front curvature. Various possibilities to choose ε are discussed, including cases where beam properties at the beam endpoint (and not at the beginning) are prescribed; for instance, the beam width at the endpoint can be specified. In such cases the beam constant is a function of the radiation angle at the source, and the decomposition of a cylindrical wave into beams has to take this into account by weighting the beams differently, at least in principle. The exact weight function is derived and shown to be reasonably well approximated by the weight function, corresponding to angle-independent ε Theoretical seismograms are presented for a laterally heterogeneous model of the crust–mantle transition which is characterized by complications in the reflection from the transition and in the refraction from below. These complications are modelled by and large with success. The seismograms, however, depend to a certain extent on the choice of the beam constant. Moreover, according to the reciprocity principle calculations with source and receiver interchanged should have the same results as calculations for the original configuration. In practice this is not so, and the difference increases with the strength of lateral heterogeneities. Hence, for a successful application of Gaussian beams the model should not vary too strongly in lateral direction.     

Long-period seismic source inversions using global tomographic models

We investigate the effect of laterally varying earth structure on centroid moment tensor inversions using fundamental mode mantle waves. Theoretical seismograms are calculated using a full formulation of surface wave ray theory. Calculations are made using a variety of global tomographic earth models. Results are compared with those obtained using the so-called great-circle approximation, which assumes that phase corrections are given in terms of mean phase slowness along the great circle, and which neglects amplitude effects of heterogeneity. Synthetic tests suggest that even source parameters which fit the data very well may have large errors due to incomplete knowledge of lateral heterogeneity. The method is applied to 31 shallow, large earthquakes. For a given earthquake, the focal mechanisms calculated using different earth models and different forward modelling techniques can significantly vary. We provide a range of selected solutions based on the fit to the data, rather than one single solution. Difficulties in constraining the dip-slip components of the seismic moment tensor often produce overestimates of seismic moment, leading to near vertical dip-slip mechanisms. This happens more commonly for earth models not fitting the data well, confirming that more accurate modelling of lateral heterogeneity can help to constrain the dip-slip components of the seismic moment tensor.     

Comparison of global synthetic seismograms calculated using the spherical 2.5-D finite-difference method with observed long-period waveforms including data from the intra-Antarctic region

We have been developing an accurate and efficient numerical scheme, which uses the finite-difference method (FDM) in spherical coordinates, for the computation of global seismic wave propagation through laterally heterogeneous realistic Earth models. In the field of global seismology, traditional axisymmetric modeling has been used widely as an efficient approach since it can solve the 3-D elastodynamic equation in spherical coordinates on a 2-D cross-section of the Earth, assuming structures to be invariant with respect to the axis through the seismic source. However, it has the severe disadvantages that asymmetric structures about the axis cannot be incorporated and the source mechanisms with arbitrary shear dislocation have not been attempted for a long time. Our scheme is based on the framework of axisymmetric modeling but has been extended to treat asymmetric structures, arbitrary moment-tensor point sources, anelastic attenuation, and the Earth center which is a singularity of wave equations in spherical coordinates. All these types of schemes which solve 3-D wavefields on a 2-D model cross-section are classified as 2.5-D modeling, so we have named our scheme the spherical 2.5-D FDM. In this study, we compare synthetic seismograms calculated using our FDM scheme with three-component observed long-period seismograms including data from stations newly installed in Antarctica in conjunction with the International Polar Year (IPY) 2007–2008. Seismic data from inland Antarctica are expected to reveal images of the Earth's deep interior with enhanced resolution because of the high signal-to-noise ratio and wide extent of this region, in addition to the rarity of sampling paths along the rotation axis of the Earth. We calculate synthetic seismograms through the preliminary reference earth model (PREM) including attenuation using a moment-tensor point source for the November 9, 2009 Fiji earthquake. Our results show quite good agreement between synthetic and observed seismograms, which indicates the accuracy of observations in the Antarctica, as well as the feasibility of the spherical 2.5-D modeling scheme.     

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.     

The accuracy of models derived by WKBJ waveform inversion

Summary . In this paper the accuracy of velocity-depth profiles derived by matching WKBJ seismograms to observations is quantitatively evaluated. Seismograms computed with the WKBJ method are generally quite reliable but possess predictable, systematic inaccuracies in the presence of strong velocity gradients. The effects of these inaccuracies on models derived through WKBJ waveform inversion are studied, using reflectivity seismograms as 'data'. The velocity structure used is an oceanic lithosphere model that contains several transition regions separated by relatively homogeneous layers, producing partially-reflected reverberations in the reflectivity synthetics that are absent from the WKBJ seismograms. The inversion incorporates the 'jumping' strategy to solve for the smoothest models consistent with the data. We find these solutions to be independent of the starting model and to have a stable basic structure that agrees well with the correct model. The differences, everywhere less than a seismic wavelength, depend on the frequency content of the seismograms. Reverberations in the reflectivity seismograms that are well separated from WKBJ arrivals are treated as 'noise' in the inversion.     

Seismic imaging of the laterally varying D" region beneath the Cocos Plate

We use an axisymmetric, spherical Earth finite difference algorithm to model SH -wave propagation through cross-sections of laterally varying lower mantle models beneath the Cocos Plate derived from recent data analyses. Synthetic seismograms with dominant periods as short as 4 s are computed for several models: (1) a D" reflector 264 km above the core–mantle boundary with laterally varying S -wave velocity increases of 0.9–2.6 per cent, based on localized structures from a 1-D double-array stacking method; (2) an undulating D" reflector with large topography and uniform velocity increase obtained using a 3-D migration method and (3) cross-sections through the 3-D mantle S -wave velocity tomography model TXBW. We apply double-array stacking to assess model predictions of data. Of the models explored, the S -wave tomography model TXBW displays the best overall agreement with data. The undulating reflector produces a double Scd arrival that may be useful in future studies for distinguishing between D" volumetric heterogeneity and D" discontinuity topography. Synthetics for the laterally varying models show waveform variability not observed in 1-D model predictions. It is challenging to predict 3-D structure based on localized 1-D models when lateral structural variations are on the order of a few wavelengths of the energy used, particularly for the grazing geometry of our data. Iterative approaches of computing synthetic seismograms and adjusting model characteristics by considering path integral effects are necessary to accurately model fine-scale D" structure.     

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.     

Numerical modelling of wavefields in three-dimensional inhomogeneous media

Summary. Numerical modelling is one of the most efficient methods for an investigation of the relationship between structural features and peculiarities of observed wavefields. It is practically the only method for 2-D and 3-D inhomogeneous media.
An algorithm based on ray theory has been developed for calculations of travel times and amplitudes of seismic waves in 3-D inhomogeneous media with curved interfaces. It was applied for numerical modelling of kinematic and dynamic characteristics of seismic waves propagating in laterally inhomogeneous media.
Travel-time and amplitude patterns were studied in the 2-D and 3-D models of a geosyncline, in which velocity distribution was given by an analytical function of the coordinates. For a more complicated model representing a subducting high-velocity lithospheric plate in a transition zone between oceanic and continental upper mantle, the velocity distribution was given by discrete values on a 2-D non-rectangular grid. It was shown that when a source was placed above the lithospheric plate, a shadow zone appeared along a strike of the structure, i.e. in the direction which is perpendicular to a strong lateral velocity gradient. Travel-time residuals were calculated along the seismological profile for a 3-D velocity distribution in the upper mantle beneath Central Asia, obtained as a result of inversion of travel times by the Backus-Gilbert method. They were found to be in a good agreement with the observed data.     

Inversion of seismic refraction data in planar dipping structure

Summary. A new method is presented for the direct inversion of seismic refraction data in dipping planar structure. Three recording geometries, each consisting of two common-shot profiles, are considered: reversed, split, and roll-along profiles. Inversion is achieved via slant stacking the common-shot wavefield to obtain a delay time–slowness (tau– p ) wavefield. The tau– p curves from two shotpoints describing the critical raypath of refracted and post-critically reflected arrivals are automatically picked using coherency measurements and the two curves are jointly used to calculate velocity and dip of isovelocity lines iteratively, thereby obtaining the final two-dimensional velocity model.
This procedure has been successfully applied to synthetic seismograms calculated for a dipping structure and to field data from central California. The results indicate that direct inversion of closely-spaced refraction/wide-aperture reflection data can practically be achieved in laterally inhomogeneous structures.     

2-Dreflectivity method and synthetic seismograms for irregularly layered structures -I. SH-wave generation

Summary. The reflectivity method for complete SH seismograms has been extended to two-dimensionally layered structures. The Aki-Larner technique is generalized to solve the integral equations for 2-D boundary conditions, and propagator matrices are enlarged to express a total SH wavefield. Synthetic seismograms in a soft basin are calculated for an incident plane-wave. They compare favourably with the results of the finite-element and finite-difference methods even in the later portion where asymptotic ray and beam theories break down. Synthetic seismograms due to a line force and a point dislocation are also presented.