Parametrizing surface wave tomographic models with harmonic spherical splines

We present a mathematical framework and a new methodology for the parametrization of surface wave phase-speed models, based on traveltime data. Our method is neither purely local, like block-based approaches, nor is it purely global, like those based on spherical harmonic basis functions. Rather, it combines the well-known theory and practical utility of the spherical harmonics with the spatial localization properties of spline basis functions. We derive the theoretical foundations for the application of harmonic spherical splines to surface wave tomography and summarize the results of numerous numerical tests illustrating the performance of a practical inversion scheme based upon them. Our presentation is based on the notion of reproducing-kernel Hilbert spaces, which lends itself to the parametrization of fully 3-D tomographic earth models that include body waves as well.     

Regional tomographic inversion of the amplitude and phase of Rayleigh waves with 2-D sensitivity kernels

In this study, we test the adequacy of 2-D sensitivity kernels for fundamental-mode Rayleigh waves based on the single-scattering (Born) approximation to account for the effects of heterogeneous structure on the wavefield in a regional surface wave study. The calculated phase and amplitude data using the 2-D sensitivity kernels are compared to phase and amplitude data obtained from seismic waveforms synthesized by the pseudo-spectral method for plane Rayleigh waves propagating through heterogeneous structure. We find that the kernels can accurately predict the perturbation of the wavefield even when the size of anomaly is larger than one wavelength. The only exception is a systematic bias in the amplitude within the anomaly itself due to a site response.
An inversion method of surface wave tomography based on the sensitivity kernels is developed and applied to synthesized data obtained from a numerical simulation modelling Rayleigh wave propagation over checkerboard structure. By comparing recovered images to input structure, we illustrate that the method can almost completely recover anomalies within an array of stations when the size of the anomalies is larger than or close to one wavelength of the surface waves. Surface wave amplitude contains important information about Earth structure and should be inverted together with phase data in surface wave tomography.     

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

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

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.     

Non-linear crustal corrections in high-resolution regional waveform seismic tomography

We compare 3-D upper mantle anisotropic structures beneath the North American continent obtained using standard and improved crustal corrections in the framework of Non-linear Asymptotic Coupling Theory (NACT) applied to long period three component fundamental and higher mode surface waveform data. Our improved approach to correct for crustal structure in high-resolution regional waveform tomographic models goes beyond the linear perturbation approximation, and is therefore more accurate in accounting for large variations in Moho topography within short distances as observed, for instance, at ocean–continent margins. This improved methodology decomposes the shallow-layer correction into a linear and non-linear part and makes use of 1-D sensitivity kernels defined according to local tectonic structure, both for the forward computation and for the computation of sensitivity kernels for inversion. The comparison of the 3-D upper mantle anisotropic structures derived using the standard and improved crustal correction approaches shows that the model norm is not strongly affected. However, significant variations are observed in the retrieved 3-D perturbations. The largest differences in the velocity models are present below 250 km depth and not in the uppermost mantle, as would be expected. We suggest that inaccurate crustal corrections preferentially map into the least constrained part of the model and therefore accurate corrections for shallow-layer structure are essential to improve our knowledge of parts of the upper mantle where our data have the smallest sensitivity.     

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

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

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.     

Formalism for traveltime inversion with finite frequency effects

A simple modification of the waveform inversion formula, based on the normal mode perturbation theory, is shown to lead to a formula for traveltime anomalies. The kernel which is derived can be used for traveltime inversion with automatic inclusion of finite frequency effects. Inversion for Earth structure with such kernels will lead to better resolution estimates than ray-theoretical traveltime inversion. Examples of kernels for transverse component seismograms are shown for direct S waves, ScS , Love waves and diffracted S waves. A measure of finite frequency effects is also proposed by comparing our formula with the one from ray theory. A quantity which should be 1 in the case of ray theory is computed for the finite frequency kernels and is shown to have deviations up to about 30 per cent from 1. Therefore, the use of ray theory for long-period body waves applies incorrect weight along a ray path and may introduce a small bias to an earth model.     

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.     

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.     

3-D linearized scattering of surface waves and a formalism for surface wave holography

Summary. Scattering of surface waves by lateral heterogeneities is analysed in the Born approximation. It is assumed that the background medium is either laterally homogeneous, or smoothly varying in the horizontal direction. A dyadic representation of the Green's function simplifies the theory tremendously. Several examples of the theory are presented. The scattering and mode conversion coefficients are shown for scattering of surface waves by the root of an Alpine-like crustal structure. Furthermore a 'great circle theorem'in a plane geometry is derived. A new proof of Snell's law is given for surface wave scattering by a quarter-space. It is shown how a stationary phase approximation can be used to simplify the Fourier synthesis of the scattered wave in the time domain. Finally a procedure is suggested to do 'surface wave holography'.     

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.     

From ACH tomographic models to absolute velocity models

The ACH method, a widely used tomographic inverse method, is characterized by the use of relative residuals in order to avoid possible biases coming from outside the target volume. The ACH method thus does not really retrieve the 3-D structure of the target volume, but instead leads to velocity contrasts relative to the layer average of the velocity, this average value remaining unknown ( Aki et al. 1977 ). Two artefacts derive from this particularity: (1) velocity contrasts are known only in the horizontal direction and it is not possible, in a strict mathematical sense, to estimate the contrasts in the vertical direction with ACH alone; (2) negative anomalies are often interpreted as low velocities, whereas negative anomalies may correspond to high velocities if the average value of the corresponding layer is sufficiently high. The converse is true of positive anomalies. We show with synthetic data how these artefacts can affect the interpretation of tomographic images. We propose to correct the artefacts by reintroducing the 1-D regional average model, and show in synthetic experiments how effective this correction can be.
  The application of this procedure to data recorded in the Kunlun region shows that the retrieval of the absolute values of the 3-D velocity model is helpful for interpreting the tomographic images and better defining which features are anomalous.     

Teleseismic waveform modelling with a one-way wave equation

It is now widely accepted that elastic properties of the continental lithosphere and the underlying sublithospheric mantle are both anisotropic and laterally heterogeneous at a range of scales. To fully exploit modern three-component broad-band array data sets requires the use of comprehensive modelling tools. In this work, we investigate the use of a wide-angle, one-way wave equation to model variations in teleseismic 3-D waveforms due to 2-D elastic heterogeneity and anisotropy. The one-way operators are derived based on a high-frequency approximation of the square-root operator and include the effects of wave propagation as well as multiple scattering. Computational cost is reduced through a number of physically motivated approximations. We present synthetic results from simple 1-D (layer over a half-space) and 2-D (subduction zone) models that are compared with reference solutions. The algorithm is then used to model data from an array of broad-band seismograph stations deployed in northwestern Canada as part of the IRIS-PASSCAL/LITHOPROBE CANOE experiment. In this region radial-component receiver functions show a clear continental Moho and the presence of crustal material dipping into the mantle at the suture of two Palaeo-Proterozoic terranes. The geometry of the suture is better defined on the transverse component where subduction is associated with a ∼10 km thick layer exhibiting strong elastic anisotropy. The modelling reproduces the main features of the receiver functions, including the effects of anisotropy, heterogeneity and finite-frequency scattering.     

Real and complex spectra – a generalization of WKBJ seismograms

Real plane-waves constitute the building blocks for recently developed spectral techniques in synthetic seismology. While providing numerical convenience, real slowness-spectra model certain wave phenomena in a distributed 'unnatural' way, whereas complex spectra model these phenomena in a compact, more 'natural' way. The theory of complex spectra, called by us the 'Spectral Theory of Transients' (STT) and developed elsewhere, is summarized here and contrasted with the real-spectrum approach. Relying strongly on the theory of analytic functions, STT permits the transient responses to be classified and evaluated according to the singularities they introduce in the complex slowness plane. The method is illustrated for a number of 2-D SH -wave model propagation environments, including interface reflection, head waves, multiple encounters with caustics due to concave boundaries or ducting medium inhomogeneities, and diffraction by structures with edges.     

Resolution potential of surface wave phase velocity measurements at small arrays

The deployment of temporary arrays of broadband seismological stations over dedicated targets is common practice. Measurement of surface wave phase velocity across a small array and its depth-inversion gives us information about the structure below the array which is complementary to the information obtained from body-wave analysis. The question is however: what do we actually measure when the array is much smaller than the wave length, and how does the measured phase velocity relates to the real structure below the array? We quantify this relationship by performing a series of numerical simulations of surface wave propagation in 3-D structures and by measuring the apparent phase velocity across the array on the synthetics. A principal conclusion is that heterogeneities located outside the array can map in a complex way onto the phase velocities measured by the array. In order to minimize this effect, it is necessary to have a large number of events and to average measurements from events well-distributed in backazimuth. A second observation is that the period of the wave has a remarkably small influence on the lateral resolution of the measurement, which is dominantly controlled by the size of the array. We analyse if the artefacts created by heterogeneities can be mistaken for azimuthal variations caused by anisotropy. We also show that if the amplitude of the surface waves can be measured precisely enough, phase velocities can be corrected and the artefacts which occur due to reflections and diffractions in 3-D structures greatly reduced.     

Full three-dimensional tomography: a comparison between the scattering-integral and adjoint-wavefield methods

This paper analyses the computational issues of full 3-D tomography, in which the starting model as well as the model perturbation is 3-D and the sensitivity (Fréchet) kernels are calculated using the full physics of 3-D wave propagation. We compare two formulations of the structural inverse problem: the adjoint-wavefield (AW) method, which back-propagates the data from the receivers to image structure, and the scattering-integral (SI) method, which sets up the inverse problem by calculating and storing the Fréchet kernels for each data functional. The two inverse methods are closely related, but which one is more efficient depends on the overall problem geometry, particularly on the ratio of sources to receivers, as well as trade-offs in computational resources, such as the relative costs of compute cycles to data storage. We find that the SI method is computationally more efficient than the AW method in regional waveform tomography using large sets of natural sources, although it requires more storage.     

Windowed f-k spectra of a three-dimensional wavefield excited by a point source in a two-dimensional multilayered elastic medium

When full 3-D modelling is too costly or cumbersome, computations of 3-D elastic wave propagation in laterally heterogeneous, multilayered 2-D geological structures may enhance considerably our ability to predict strong ground motion for seismological and engineering purposes. Towards this goal, we extend the method based on the combination of the thin-layer finite-element and boundary-element methods (TLFE-BEM) and calculate windowend f - k spectra of the 3-D wavefield. The windowed f - k spectra are spatially localized spectra from which the local properties of the wavefield can be extracted. The TLFE-BEM is particularly suited for calculating the complete wavefield where surface waves are dominant in multilayered media. The computations are performed in the frequency domain, providing the f - k spectra directly. From the results for the 3-D wavefield excited by a point source in a 2-D multilayered, sloped structure, it can be said that the phase velocity of the fundamental-mode Rayleigh wave in a laterally heterogeneous multilayered medium, estimated from the windowed f - k spectra, varies with the location of the point source. For the model calculated in this article, the phase velocity varies between the value for the flat layered structure of the thick-layer side and that for the structure just under the centre of the window. The exact subsurface structure just under the centre of an array in a laterally heterogeneous medium cannot be obtained if we use the f - k spectral analysis assuming a flat layered structure.     

Galerkin boundary integral equation method for spontaneous rupture propagation problems: SH-case

We develop a Galerkin finite element boundary integral equation method (GaBIEM) for spontaneous rupture propagation problems for a planar fault embedded in a homogeneous full 2-D space. A 2-D antiplane rupture propagation problem, with a slip-weakening friction law, is simulated by the GaBIEM. This method allows one to eliminate the strong singularities from the integral representation of the traction, and to separate explicitly the expression for the traction into an instantaneous component; static and time-dependent components with weakly (logarithmic) singular kernels; and a dynamic component and a quasi-static component, with continuous, bounded, kernels. Simulated results throw light into the performance of the GaBIEM and highlight differences with respect to that of the traditional, collocation, boundary integral equation method (BIEM). Both methods converge with a power law with respect to grid size, with different exponents. There is no restriction on the CFL stability number for the GaBIEM since an implicit, unconditionally stable method is used for the time integration. The error of the approximation increases with the time step, as expected, and it can remain below that of the BIEM.