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 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.     

A geometrical approach to time evolving wave fronts

The parameter that defines the ray tracing equations in the direct geometrical approach is the product of the radius of curvature of the wave front by the velocity on the wave front ( RV ). To show this, we derive motion equations for the centre and the radius of curvature of an expanding wave front. The continuity of RV along rays implies Snell's Law. For constant velocities the equation for the radius of curvature reduces to the original Huygens' Principle. The variable RV can be computed during ray tracing and used to determine the local radius of curvature, which in turn can be used in geometrical spreading, amplitude corrections and structure interpretation.     

An introduction to Maslov's asymptotic method

Summary. Familiar concepts such as asymptotic ray theory and geometrical spreading are now recognized as an asymptotic form of a more general asymptotic solution to the non-separable wave equation. In seismology, the name Maslov asymptotic theory has been attached to this solution. In its simplest form, it may be thought of as a justification of disc-ray theory and it can be reduced to the WKBJ seismogram. It is a uniformly valid asymptotic solution, though. The method involves properties of the wavefronts and ray paths of the wave equation which have been established for over a century. The integral operators which build on these properties have been investigated only comparatively recently. These operators are introduced very simply by appealing to the asymptotic Fourier transform of Ziolkowski & Deschamps. This leads quite naturally to the result that phase functions in different domains of the spatial Fourier transform are related by a Legendre transformation. The amplitude transformation can also be inferred by this method. Liouville's theorem (the incompressibility of a phase space of position and slowness) ensures that it is always possible to obtain a uniformly asymptotic solution. This theorem can be derived by methods familiar to seismologists and which do not rely on the traditional formalism of classical mechanics. It can also be derived from the sympletic property of the equations of geometrical spreading and canonical transformations in general. The symplectic property plays a central role in the theory of high-frequency beams in inhomogeneous media.     

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.     

Phase shift at caustics along rays in anisotropic media

In isotropic ray tracing, the ray approximation to the wavefield undergoes a phase shift when the ray crosses a caustic. The cumulative number of such phase shifts along a ray is usually called the KMAH index. The sign of these phase shifts is prescribed by the sign of the angular frequency in combination with the sign convention used for the Fourier transformation. In isotropic media the KMAH index always increases by one or by two, depending on the type of caustic crossed. For (quasi-)shear waves in anisotropic media the KMAH index may decrease. This is the case if the associated slowness sheet is locally concave in one or two of its principal directions of curvature.     

Retrieving shear velocity structure from high-frequency toroidal modes of the Earth

Summary. For a smooth earth model, observations of a set of high-frequency toroidal modes at fixed slowness yield only a single piece of information, the tau value for that slowness. In this note, a procedure for obtaining the shear velocity structure from free oscillation data for an earth model with velocity discontinuities is developed, based on the method of tau inversion. The information content of the high-frequency modes is greater in this case, and the nature and depths of the discontinuities may be deduced. It is shown, for the real Earth, that the tau values obtained from free oscillation data are affected significantly by the presence of the Moho, but a simple iterative scheme may be used to remove this contamination. Brune's method of deducing mode frequencies from body wave pulses is shown to produce significant errors for a model with a pronounced Moho discontinuity, and the same iterative scheme may also be employed to correct for this effect.     

Seismic ray theory for lithospheric structures with slight lateral variations

b
In this paper, the method of small perturbations is applied to the ray and energy transport equations in an investigation of the effect of weak inhomogeneities on the propagation of seismic rays through a layer of fixed thickness. Just as Aki et al . have used travel-time residuals to infer first-order velocity perturbations in their block-modelling procedure, it is proposed that surface slowness and amplitude data may be used to give additional information about the structure of the velocity perturbations beneath the observer.     

On the computation of shear waves by the ray method

Summary. The ray series solution of the elastodynamic equation of motion for shear waves propagating through a laterally inhomogeneous three-dimensional medium can be simplified by the use of a particular coordinate system that accompanies the wave front along the ray of investigation. The system is entirely determined by parameters that are obtainable from the ray. The transport equations for the principal shear wave components are then no longer coupled, but reduce to the same type of equation which determines the principal compressional wave component.     

Generalized Born scattering of elastic waves in 3-D media

It is well known that when a seismic wave propagates through an elastic medium with gradients in the parameters which describe it (e.g. slowness and density), energy is scattered from the incident wave generating low-frequency partial reflections. Many approximate solutions to the wave equation, e.g. geometrical ray theory (GRT), Maslov theory and Gaussian beams, do not model these signals. The problem of describing partial reflections in 1-D media has been extensively studied in the seismic literature and considerable progress has been made using iterative techniques based on WKBJ, Airy or Langer type ansätze. In this paper we derive a first-order scattering formalism to describe partial reflections in 3-D media. The correction term describing the scattered energy is developed as a volume integral over terms dependent upon the first spatial derivatives (gradients) of the parameters describing the medium and the solution. The relationship we derive could, in principle, be used as the basis for an iterative scheme but the computational expense, particularly for elastic media, will usually prohibit this approach. The result we obtain is closely related to the usual Born approximation, but differs in that the scattering term is not derived from a perturbation to a background model, but rather from the error in an approximate Green's function. We examine analytically the relationship between the results produced by the new formalism and the usual Born approximation for a medium which has no long-wavelength heterogeneities. We show that in such a case the two methods agree approximately as expected, but that in a media with heterogeneities of all wavelengths the new gradient scattering formalism is superior. We establish analytically the connection between the formalism developed here and the iterative approach based on the WKBJ solution which has been used previously in 1-D media. Numerical examples are shown to illustrate the examples discussed.     

Validity of the great circular average approximation for inversion of normal mode measurements

Summary. Synthetic seismograms based upon first-order perturbation theory are analysed to test the validity of assumptions which form the basis of current velocity inversion procedures. It is found that the lowest order geometrical optics approximation, namely that measured normal mode eigen-frequencies reflect the average structure underlying the source–receiver great circle path, becomes less valid near nodes in the source radiation pattern and near the surface wave foci at the source and its antipode. These failures are a consequence of singlet interference within an isolated normal mode multiplet. The technique of determing frequency by fitting a single resonance peak to a multiplet yields results which agree well with the first-order theory for slow and fast paths where excitation is dominated by one pair of singlets but on intermediate paths where singlet interference is more of a problem, agreement is not as good. Inversion of small data sets is particularly sensitive to frequency fluctuations near radiation nodes, while larger sets are influenced more by antipodal deviations from geometrical optics. The latter leads to inversions which fail to recover the short wavelength structure of the starting model. Basing inversions directly upon first-order theory shows promise of improving recovery of short wavelengths.     

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.     

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.     

A method of comparison of exact and asymptotic wave field computations

Summary. A method of comparison of exact numerical computations with an asymptotic ray series expansion consisting of the two first terms is proposed. The method makes it unnecessary to derive complicated explicit expressions for the second leading term of the ray series.
As a practical example we consider the anomalous PS arrival generated in the case of a near-vertical incidence of a spherical P wave on a solid/solid boundary. The areas in which the PS wave may be described by two leading terms of the ray series expansion are marked and deviations from the ray theory are analysed.     

Some useful approximations to generalized ray theory for regional distance seismograms

Summary . Seismograms recorded at regional distances (2°–12°) are quite complicated due to the waveguide nature of the crust. Generalized ray theory can be used to model the body waves in this distance range but a very large number of rays is required. Here I present a series of approximations to streamline generalized ray theory for the waveguide problem. If a layer over a half-space is used for the structure, then the de Hoop contour for a given ray is most strongly dependent on the fastest velocity of any leg of the ray. This results in analytic approximations to locate the contour. Each ray has two body wave arrivals (a headwave and a reflected arrival) so the displacement response of the ray need only be evaluated at a few points in time about the two arrival times and interpolated in between. A change in structure (increasing crustal thickness or Pn velocity) most strongly affects the relative timing of the headwave and the reflected arrival, so it is possible to 'stretch' or 'squeeze' the waveform of a representative model to simulate a whole suite of models.
Also discussed is the applicability of a single layer over a half-space structure for modelling the observed regional distance waveforms for shallow earthquakes. At periods greater than a few seconds crustal layering can be replaced by a single layer having the appropriate average velocities. Lateral variations in crustal thickness with scale lengths of less than about 100 km can also be modelled with a simple horizontal layer of appropriate average thickness.     

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.     

Born scattering of longperiod body waves

The Born approximation is applied to the modelling of the propagation of deeply turning longperiod body waves through heterogeneities in the lowermost mantle. We use an exact Green's function for a spherically symmetric earth model that also satisfies the appropriate boundary conditions at internal boundaries and the surface of the earth. The scattered displacement field is obtained by a numerical quadrature of the product of the Green's function, the exciting wavefield and structural perturbations. We study three examples: scattering of longperiod P waves from a plume rising from the coremantle boundary (CMB), generation of longperiod precursors to PKIKP by strong, localized scatterers at the CMB, and propagation of corediffracted P waves through largescale heterogeneities in D". The main results are as follows: (1) the signals scattered from a realistic plume are small with relative amplitudes of less than 2 per cent at a period of 20 s, rendering plume detection a fairly difficult task; (2) strong heterogeneities at the CMB of appropriate size may produce observable longperiod precursors to PKIKP in spite of the presence of a diffraction from the PKP B caustic; (3) corediffracted P  waves ( P _diff) are sensitive to structure in D" far off the geometrical ray path and also far beyond the entry and exit points of the ray into and out of D"; sensitivity kernels exhibit ringshaped patterns of alternating sign reminiscent of Fresnel zones; (4) P _diff also shows a nonnegligible sensitivity to shear wave velocity in D"; (5) down to periods of 40 s, the Born approximation is sufficiently accurate to allow waveform modelling of P _diff through largescale heterogeneities in D" of up to 5 per cent.     

Asymptotic eigenfrequencies of the Earth's normal modes

We derive asymptotic formulae for the toroidal and spheroidal eigenfrequencies of a SNREI earth model with two discontinuities, by considering the constructive interference of propagating SH and P-SV body waves. For a model with a smooth solid inner core, fluid outer core and mantle, there are four SH and 10 P-SV ray parameters regimes, each of which must be examined separately. The asymptotic eigenfrequency equations in each of these regimes depend only on the intercept times of the propagating wave types and the reflection and transmission coefficients of the waves at the free surface and the two discontinuities. If the classical geometrical plane-wave reflection and transmission coefficients are used, the final eigenfrequency equations are all real. In general, the asymptotic eigenfrequencies agree extremely well with the exact numerical eigenfrequencies; to illustrate this, we present comparisons for a crustless version of earth model 1066A.     

Slowness-backazimuth weighted migration: A new array approach to a high-resolution image

We have developed a new array method combining conventional migration with a slowness-backazimuth deviation weighting scheme. All seismic traces are shifted based on the theoretical traveltime of the scattered wave from specific gridpoints in a 3-D volume. Observed slowness and backazimuth are calculated for each raypath and compared with theoretical values in order to estimate slowness and backazimuth deviations. Subsequently, stacked energy calculated by a conventional migration method is weighted by the slowness and backazimuth deviations to suppress any arrival energy whose slowness and backazimuth are inconsistent with the expected theoretical values. This new method was applied to two P- wave data sets which comprise (1) underside reflections at the 410 and 660 km mantle discontinuities and (2) D" reflections as well as their corresponding synthetic data sets. The results show that the weighting scheme dramatically increases the resolution of the migrated images and enables us to obtain well-constrained, focused images, making upper-mantle discontinuities and D" reflections more distinct by reducing their surrounding energy.