Finite-frequency sensitivity kernels for head waves

Head waves are extremely important in determining the structure of the predominantly layered Earth. While several recent studies have shown the diffractive nature and the 3-D Fréchet kernels of finite-frequency turning waves, analogues of head waves in a continuous velocity structure, the finite-frequency effects and sensitivity kernels of head waves are yet to be carefully examined. We present the results of a numerical study focusing on the finite-frequency effects of head waves. Our model has a low-velocity layer over a high-velocity half-space and a cylindrical-shaped velocity perturbation placed beneath the interface at different locations. A 3-D finite-difference method is used to calculate synthetic waveforms. Traveltime and amplitude anomalies are measured by the cross-correlation of synthetic seismograms from models with and without the velocity perturbation and are compared to the 3-D sensitivity kernels constructed from full waveform simulations. The results show that the head wave arrival-time and amplitude are influenced by the velocity structure surrounding the ray path in a pattern that is consistent with the Fresnel zones. Unlike the 'banana–doughnut' traveltime sensitivity kernels of turning waves, the traveltime sensitivity of the head wave along the ray path below the interface is weak, but non-zero. Below the ray path, the traveltime sensitivity reaches the maximum (absolute value) at a depth that depends on the wavelength and propagation distance. The sensitivity kernels vary with the vertical velocity gradient in the lower layer, but the variation is relatively small at short propagation distances when the vertical velocity gradient is within the range of the commonly accepted values. Finally, the depression or shoaling of the interface results in increased or decreased sensitivities, respectively, beneath the interface topography.     

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.     

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.     

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.     

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.     

Effects of slight anisotropy on surface waves

We present a complete ray theory for the calculation of surface-wave observables from anisotropic phase-velocity maps. Starting with the surface-wave dispersion relation in an anisotropic earth model, we derive practical dynamical ray-tracing equations. These equations allow calculation of the observables phase, arrival-angle and amplitude in a ray theoretical framework. Using perturbation theory, we also obtain approximate expressions for these observables. We assess the accuracy of the first-order approximations by using both theories to make predictions on a sample anisotropic phase-velocity map. A comparison of the two methods illustrates the size and type of errors which are introduced by perturbation theory. Perturbation theory phase and arrival-angle predictions agree well with the exact calculation, but amplitude predictions are poor. Many previous studies have modelled surface-wave propagation using only isotropic structure, not allowing for anisotropy. We present hypothetical examples to simulate isotropic modelling of surface waves which pass through anisotropic material. Synthetic data sets of phase and arrival angle are produced by ray tracing with exact ray theory on anisotropic phase-velocity maps. The isotropic models obtained by inverting synthetic anisotropic phase data sets produce deceptively high variance reductions because the effects of anisotropy are mapped into short-wavelength isotropic structure. Inversion of synthetic arrival-angle data sets for isotropic models results in poor variance reductions and poor recovery of the isotropic part of the anisotropic input map. Therefore, successful anisotropic phase-velocity inversions of real data require the inclusion of both phase and arrival-angle measurements.     

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.     

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.     

Transient and impulse responses of a one-dimensional linearly attenuating medium — I. Analytical results

Summary. The transient and impulse responses (Green's function) for onedimensional wave propagation in a standard linear solid are calculated using a Laplace Transform method. The spectrum of relaxation times is chosen so as to model a constant Q medium within an absorption band covering a broad frequency range which may be chosen so as to include the seismic frequencies. The inverse transform may be evaluated asymptotically in the limit of very long propagation times using the saddle point method. For shorter propagation times the method of steepest descent may be modified so as to yield an accurate first motion approximation. The character of the small amplitude precursor to the large amplitude Visible' signal is investigated analytically. It is shown that the signal velocity is intermediate between the high-frequency ('unrelaxed') and the low-frequency ('relaxed') limits of the phase velocity.     

Global multiresolution models of surface wave propagation: comparing equivalently regularized Born and ray theoretical solutions

I invert a large set of teleseismic phase-anomaly observations, to derive tomographic maps of fundamental-mode surface wave phase velocity, first via ray theory, then accounting for finite-frequency effects through scattering theory, in the far-field approximation and neglecting mode coupling. I make use of a multiple-resolution pixel parametrization which, in the assumption of sufficient data coverage, should be adequate to represent strongly oscillatory Fréchet kernels. The parametrization is finer over North America, a region particularly well covered by the data. For each surface-wave mode where phase-anomaly observations are available, I derive a wide spectrum of plausible, differently damped solutions; I then conduct a trade-off analysis, and select as optimal solution model the one associated with the point of maximum curvature on the trade-off curve. I repeat this exercise in both theoretical frameworks, to find that selected scattering and ray theoretical phase-velocity maps are coincident in pattern, and differ only slightly in amplitude.     

Simultaneous velocity and interface tomography of normal-incidence and wide-aperture seismic traveltime data

A tomographic inversion technique that inverts traveltimes to obtain a model of the subsurface in terms of velocities and interfaces is presented. It uses a combination of refraction, wide-angle reflection and normal-incidence data, it simultaneously inverts for velocities and interface depths, and it is able to quantify the errors and trade-offs in the final model. The technique uses an iterative linearized approach to the non-linear traveltime inversion problem. The subsurface is represented as a set of layers separated by interfaces, across which the velocity may be discontinuous. Within each layer the velocity varies in two dimensions and has a continuous first derivative. Rays are traced in this medium using a technique based on ray perturbation theory, and two-point ray tracing is avoided by interpolating the traveltimes to the receivers from a roughly equidistant fan of rays. The calculated traveltimes are inverted by simultaneously minimizing the misfit between the data and calculated traveltimes, and the roughness of the model. This 'smoothing regularization' stabilizes the solution of the inverse problem. In practice, the first iterations are performed with a high level of smoothing. As the inversion proceeds, the level of smoothing is gradually reduced until the traveltime residual is at the estimated level of noise in the data. At this point, a minimum-feature solution is obtained, which should contain only those features discernible over the noise.
The technique is tested on a synthetic data set, demonstrating its accuracy and stability and also illustrating the desirability of including a large number of different ray types in an inversion.     

Plane wave propagation in nonlinear-elastic anisotropic media

Summary. Kelvin-Christoffel equations describing plane wave propagation in anisotropic media are generalized to account for the effects of nonlinear elasticity. The polarization and waveform of nonlinear distortions of a transient plane wave are investigated by means of perturbation theory. Detailed analysis for an anisotropic medium with hexagonal symmetry shows that for "pure" shear-waves the polarization vector of the nonlinear component is always perpendicular to that of the linear wave. In the case of a high-amplitude excitation (for instance, in the vicinity of large earthquakes) the influence of nonlinearity may cause distortions of shear-wave polarization, which contains the most reliable information on the presence and characteristics of anisotropy. The solutions presented in this paper make it possible to solve reflection-transmission problems in nonlinear-elastic anisotropic media.     

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.     

On crustal corrections in surface wave tomography

Mantle models from surface waves rely on good crustal corrections. We investigated how far ray theoretical and finite frequency approximations can predict crustal corrections for fundamental mode surface waves. Using a spectral element method, we calculated synthetic seismograms in transversely isotropic PREM and in the 3-D crustal model Crust2.0 on top of PREM, and measured the corresponding time-shifts as a function of period. We then applied phase corrections to the PREM seismograms using ray theory and finite frequency theory with exact local phase velocity perturbations from Crust2.0 and looked at the residual time-shifts. After crustal corrections, residuals fall within the uncertainty of measured phase velocities for periods longer than 60 and 80 s for Rayleigh and Love waves, respectively. Rayleigh and Love waves are affected in a highly non-linear way by the crustal type. Oceanic crust affects Love waves stronger, while Rayleigh waves change most in continental crust. As a consequence, we find that the imperfect crustal corrections could have a large impact on our inferences of radial anisotropy. If we want to map anisotropy correctly, we should invert simultaneously for mantle and crust. The latter can only be achieved by using perturbation theory from a good 3-D starting model, or implementing full non-linearity from a 1-D starting model.     

Computation of Large Anisotropic Seismic Heterogeneities (CLASH)

A general tomographic technique is designed in order (i) to operate in anisotropic media; (ii) to account for the uneven seismic sampling and (iii) to handle massive data sets in a reasonable computing time. One modus operandi to compute a 3-D body wave velocity model relies on surface wave phase velocity measurements. An intermediate step, shared by other approaches, consists in translating, for each period of a given mode branch, the phase velocities integrated along ray paths into local velocity perturbations. To this end, we develop a method, which accounts for the azimuthal anisotropy in its comprehensive form. The weakly non-linear forward problem allows to use a conjugate gradient optimization. The Earth's surface is regularly discretized and the partial derivatives are assigned to the individual grid points. Possible lack of lateral resolution, due to the inescapable uneven ray path coverage, is taken into account through the a priori covariances on parameters with laterally variable correlation lengths. This method allows to efficiently separate the 2ψ and the 4ψ anisotropic effects from the isotropic perturbations. Fundamental mode and overtone phase velocity maps, derived with real Rayleigh wave data sets, are presented and compared with previous maps. The isotropic models concur well with the results of Trampert & Woodhouse. Large 4ψ heterogeneities are located in the tectonically active regions and over the continental lithospheres such as North America, Antarctica or Australia. At various periods, a significant 4ψ signature is correlated with the Hawaii hotspot track. Finally, concurring with the conclusions of Trampert & Woodhouse, our phase velocity maps show that Rayleigh wave data sets do need both 2ψ and 4ψ anisotropic terms.     

Influence of the Interface Fresnel zone on the reflected P-wave amplitude modelling

The aim of the paper is to emphasize the importance of accounting for the Fresnel volume and for the Interface Fresnel zone (IFZ) for calculating the amplitude of the P wave emanating from a point source and recorded at a receiver after its specular reflection on a smooth homogeneous interface between elastic media. For this purpose, by considering the problem of interest as a problem of diffraction by the IFZ, that is, the physically relevant part of the interface which actually affects the reflected wavefield, we have developed a method which combines the Angular Spectrum Approach (ASA) with the IFZ concept to get the 3-D analytical solution. The variation in the reflected P -wave amplitude evaluated with the ASA, as a function of the incidence angle, is compared with the plane wave (PW) reflection coefficient and with the exact solution provided by the 3-D code OASES, for one solid/solid configuration and two dominant frequencies of the source. For subcritical incidence angles the geometrical spreading compensation is mostly quite sufficient to reduce the point-source amplitudes to the PW amplitudes. On the contrary, for specific regions of incidence angles for which the geometrical spreading compensation is not sufficient anymore, that is, near the critical region and in the post-critical domain, the ASA combined with the IFZ concept yields better results than the PW theory whatever the dominant frequency of the source, which suggests that the additional application of the IFZ concept is necessary to obtain the reflected P -wave amplitude. Nevertheless, as the ASA combined with the IFZ has been used only for evaluating the contribution of the reflected wavefield at the receiver, its predictions fail when the interference between the reflected wave and the head wave becomes predominant.     

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.     

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.     

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.