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.     

Linear programming approach to moment tensor inversion of earthquake sources and some tests on the three-dimensional structure of the upper mantle

Summary. A new method of moment tensor inversion is developed, which combines surface wave data and P -wave first motion data in a linear programming approach. Once surface wave spectra and first motion data are given, the method automatically obtains the solution that satisfies first motion data and minimizes the L1 norm of the surface wave spectra. We show the results of eight events in which the method works and is stable even for shallow events. We also show one event in which surface wave data and P -wave first motion data seem to be incompatible. In such cases, our method does not converge or converges to a solution which has a large minor (second) double couple component. It is an advantage that the method can determine the compatibility of two data sets without trial and error.
Laterally heterogeneous phase velocity corrections are used to obtain spectra at the source. The method is also applied to invert moment tensors of eight events in two recent three-dimensional (3-D) upper mantle structures. In both 3-D models, variances of spectra are smaller than those in a laterally homogeneous model at 256 s. Statistical tests show that those reductions are significant at a high confidence level for five events out of eight examined. For three events, we examined those reductions at shorter periods, 197 and 151 s. The reduction of variances is comparable to the results at 256 s and is again statistically significant at a high confidence level. Orientation of fault planes does not change very much by incorporation of lateral variations of phase velocity or by doing inversions at different periods. This is mainly because of the constraints from P -wave first motion data. Scatter of phase spectra at shorter periods, especially at 151 s, is great and suggests that surface wave ray paths deviate from great circle paths substantially and these effects cannot be ignored.     

Determination of three-dimensional velocity structure from observations of refracted body waves

Summary. The computational effectiveness of travel-time inversion methods depends on the parameterization of a 3-D velocity structure. We divide a region of interest into a few layers and represent the perturbation of wave slowness in each layer by a series of Chebyshev polynomials. Then a relatively complex velocity structure can be dcscribed by a small set of parameters that can be accurately evaluated by a linearized inversion of travel-time residuals. This method has been applied to artificial and real data at small epicentral distances and in the teleseismic distance range. The corresponding matrix equations were solved using singular value decomposition. The results suggest that the method combines resolution with computational convenience.     

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.     

A WKBJ spectral method for computation of SV synthetic seismograms in a cylindrically symmetric medium

Summary. A method of synthetic seismogram computation for teleseismic SV -waves is developed in order to treat quantitatively SV -waves in problems of body wave source inversion and source—receiver structure studies. The method employs WKBJ theory for a generalized ray in a vertically inhomogeneous half-space and the propagator matrix technique for waves in near-surface homogeneous layers. Wavenumber integration is done along the real axis of the wavenumber plane and anelasticity is included by using complex velocity in all regions of the earth model. The near-surface source structure is taken into account in the computation for the case of the shallow source by allowing a point source to be located in the homogeneous layers. Source and receiver area structures are also allowed to differ. A general moment tensor point source is considered.     

Surface wave tomography of the Barents Sea and surrounding regions

The goal of this study is to refine knowledge of the structure and tectonic history of the European Arctic using the combination of all available seismological surface wave data, including historical data that were not used before for this purpose. We demonstrate how the improved data coverage leads to better depth and spatial resolution of the seismological model and discovery of intriguing features of upper-mantle structure. To improve the surface wave data set in the European Arctic, we extensively searched for broad-band data from stations in the area from the beginning of the 1970s until 2005. We were able to retrieve surface wave observations from regional data archives in Norway, Finland, Denmark and Russia in addition to data from the data centres of IRIS and GEOFON. Rayleigh and Love wave group velocity measurements between 10 and 150 s period were combined with existing data provided by the University of Colorado at Boulder. This new data set was inverted for maps showing the 2-D group-velocity distribution of Love and Rayleigh waves for specific periods. Using Monte Carlo inversion, we constructed a new 3-D shear velocity model of the crust and upper mantle beneath the European Arctic which provides higher resolution and accuracy than previous models. A new crustal model of the Barents Sea and surrounding areas, published recently by a collaboration between the University of Oslo, NORSAR and the USGS, constrains the 3-D inversion of the surface wave data in the shallow lithosphere. The new 3-D model, BARMOD, reveals substantial variations in shear wave speeds in the upper mantle across the region with a nominal resolution of 1°× 1°. Of particular note are clarified images of the mantle expression of the continent-ocean transition in the Norwegian Sea and a deep, high wave speed lithospheric root beneath the Eastern Barents Sea, which presumably is the remnant of several Palaeozoic collisions.     

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.     

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.     

Seismic tomographic imaging of P- and S-waves velocity perturbations in the upper mantle beneath Iran

The inverse tomography method has been used to study the P - and S -waves velocity structure of the crust and upper mantle underneath Iran. The method, based on the principle of source–receiver reciprocity, allows for tomographic studies of regions with sparse distribution of seismic stations if the region has sufficient seismicity. The arrival times of body waves from earthquakes in the study area as reported in the ISC catalogue (1964–1996) at all available epicentral distances are used for calculation of residual arrival times. Prior to inversion we have relocated hypocentres based on a 1-D spherical earth's model taking into account variable crustal thickness and surface topography. During the inversion seismic sources are further relocated simultaneously with the calculation of velocity perturbations. With a series of synthetic tests we demonstrate the power of the algorithm and the data to reconstruct introduced anomalies using the ray paths of the real data set and taking into account the measurement errors and outliers. The velocity anomalies show that the crust and upper mantle beneath the Iranian Plateau comprises a low velocity domain between the Arabian Plate and the Caspian Block. This is in agreement with global tomographic models, and also tectonic models, in which active Iranian plateau is trapped between the stable Turan plate in the north and the Arabian shield in the south. Our results show clear evidence of the mainly aseismic subduction of the oceanic crust of the Oman Sea underneath the Iranian Plateau. However, along the Zagros suture zone, the subduction pattern is more complex than at Makran where the collision of the two plates is highly seismic.     

Sequential integrated inversion of refraction and wide-angle reflection traveltimes and gravity data for two-dimensional velocity structures

A new algorithm is presented for the integrated 2-D inversion of seismic traveltime and gravity data. The algorithm adopts the 'maximum likelihood' regularization scheme. We construct a 'probability density function' which includes three kinds of information: information derived from gravity measurements; information derived from the seismic traveltime inversion procedure applied to the model; and information on the physical correlation among the density and the velocity parameters. We assume a linear relation between density and velocity, which can be node-dependent; that is, we can choose different relationships for different parts of the velocity–density grid. In addition, our procedure allows us to consider a covariance matrix related to the error propagation in linking density to velocity. We use seismic data to estimate starting velocity values and the position of boundary nodes. Subsequently, the sequential integrated inversion (SII) optimizes the layer velocities and densities for our models. The procedure is applicable, as an additional step, to any type of seismic tomographic inversion.
We illustrate the method by comparing the velocity models recovered from a standard seismic traveltime inversion with those retrieved using our algorithm. The inversion of synthetic data calculated for a 2-D isotropic, laterally inhomogeneous model shows the stability and accuracy of this procedure, demonstrates the improvements to the recovery of true velocity anomalies, and proves that this technique can efficiently overcome some of the limitations of both gravity and seismic traveltime inversions, when they are used independently.
An interpretation of field data from the 1994 Vesuvius test experiment is also presented. At depths down to 4.5 km, the model retrieved after a SII shows a more detailed structure than the model obtained from an interpretation of seismic traveltime only, and yields additional information for a further study of the area.     

Inversion for interface structure using teleseismic traveltime residuals

An inversion method is presented for the reconstruction of interface geometry between two or more crustal layers from teleseismic traveltime residuals. The method is applied to 2-D models consisting of continuous interfaces separating constant-velocity layers. The forward problem of determining ray paths and traveltimes between incident wave fronts below the structure and receivers located on the Earth's surface is solved by an efficient and robust shooting method. A conjugate gradient method is employed to solve the inverse problem of minimizing a least-squares type objective function based on the difference between observed and calculated traveltimes. Teleseismic data do not accurately constrain average vertical structure, so a priori information in the form of layer velocities and average layer thicknesses is required. Synthetic tests show that the method can be used to reconstruct interface geometry accurately, even in the presence of data noise. Tests also show that, if layer velocities and initial interface positions are poorly chosen, lateral structure is still recoverable. The inversion method was applied to previously published teleseismic data recorded by an in-line array of portable seismographs that traversed the northern margin of the Musgrave Block, central Australia. The solution based on interface parametrization is consistent with models given by other studies that used the same data but different methods, most notably the standard tomographic approach that inverts for velocity rather than interface structure.     

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.     

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.     

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.     

Three-dimensional image of the Moho undulations beneath the Gulf of Bothnia using wide-angle seismic data

The BABEL marine seismic experiment has been carried out to investigate the lithospheric structure and antecedent tectonic signatures of the Baltic Shield, including the Archaean-Proterozoic collisional structure in the northern part of the Gulf of Bothnia.
Lithospheric seismic-reflection streamer data and simultaneously recorded wide-angle reflection and refraction data collected in the Gulf of Bothnia as part of the BABEL project have been used for 3-D modelling. The distribution of land stations around the Gulf provides a good 3-D ray coverage of the PMP reflection data recorded at the eight stations in the area and allows an estimation of strikes and dips of the Moho boundary in the area. The traveltimes of reflected phases are calculated using a method that utilizes the finite-difference solution of the eikonal equation. The Moho wide-angle-reflection (PMP) traveltimes are modelled using an inversion method. A 2-D model from the Gulf of Bothnia extended into the third dimension is used as an initial model. During the inversion the velocity is kept constant and only the Moho boundary is allowed to vary. To estimate the strike of the Moho boundary and the stability of the inversion, two initial models with different strikes are examined.
The results indicate that the Moho depth in the Gulf of Bothnia undulates and has a maximum depth of 55 km in the south, rising to 42 km in the north. The Moho depth variations seem to be step-like. This change in the Moho depth coincides with the location of the presumed fossil subduction zone in the area. The crustal-thickness variations seem to be well approximated by a nearly 2-D structure striking parallel to a postulated subduction zone immediately to the south of the Skellefte area. The presence of the step at the crust/mantle boundary can be interpreted as a result of a plate-collision event at about 2 Ga.     

P- and S-velocity images of the lithosphere–asthenosphere system in the Central Andes from local-source tomographic inversion

About 50 000 P and S arrival times and 25 000 values of t * recorded at seismic arrays operated in the Central Andes between 20°S and 25°S in the time period from 1994 to 1997 have been used for locating more than 1500 deep and crustal earthquakes and creating 3-D P , S velocity and Qp models. The study volume in the reference model is subdivided into three domains: slab, continental crust and mantle wedge. A starting velocity distribution in each domain is set from a priori information: in the crust it is based on the controlled sources seismic studies; in slab and mantle wedge it is defined using relations between P and S velocities, temperature and composition given by mineral physics. Each iteration of tomographic inversion consists of the following steps: (1) absolute location of sources in 3-D velocity model using P and S arrival times; (2) double-difference relocation of the sources and (3) simultaneous determination of P and S velocity anomalies, P and S station corrections and source parameters by inverting one matrix. Velocity parameters are computed in a mesh with the density of nodes proportional to the ray density with double-sided nodes at the domain boundaries. The next iteration is repeated with the updated velocity model and source parameters obtained at the previous step. Different tests aimed at checking the reliability of the obtained velocity models are presented. In addition, we present the results of inversion for Vp and Vp/Vs parameters, which appear to be practically equivalent to Vp and Vs inversion. A separate inversion for Qp has been performed using the ray paths and source locations in the final velocity model. The resulting Vp , Vs and Qp distributions show complicated, essentially 3-D structure in the lithosphere and asthenosphere. P and S velocities appear to be well correlated, suggesting the important role of variations of composition, temperature, water content and degree of partial melting.     

Three-dimensional massively parallel electromagnetic inversion—I. Theory

An iterative solution to the non-linear 3-D electromagnetic inverse problem is obtained by successive linearized model updates using the method of conjugate gradients. Full wave equation modelling for controlled sources is employed to compute model sensitivities and predicted data in the frequency domain with an efficient 3-D finite-difference algorithm. Necessity dictates that the inverse be underdetermined, since realistic reconstructions require the solution for tens of thousands of parameters. In addition, large-scale 3-D forward modelling is required and this can easily involve the solution of over several million electric field unknowns per solve. A massively parallel computing platform has therefore been utilized to obtain reasonable execution times, and results are given for the 1840-node Intel Paragon. The solution is demonstrated with a synthetic example with added Gaussian noise, where the data were produced from an integral equation forward-modelling code, and is different from the finite difference code embedded in the inversion algorithm     

Seismic reflection traveltimes in two-dimensional statistically anisotropic random media

Velocity estimation remains one of the main problems when imaging the subsurface with seismic reflection data. Traveltime inversion enables us to obtain large-scale structures of the velocity field and the position of seismic reflectors. However, as the media currently under study are becoming more and more complex, we need to know the finer-scale structures. The problem is that below a certain range of velocity heterogeneities, deterministic methods become difficult to use, so we turn to a probabilistic approach. With this in view, we characterize the velocity field as a random field defined by its first and second statistical moments. Usually, a seismic random medium is defined as a homogeneous velocity background perturbed by a small random field that is assumed to be stationary. Thus, we make a link between such a random velocity medium (together with a simple reflector) and seismic reflection traveltimes. Assuming that the traveltimes are ergodic, we use 2-D seismic reflection geometry to study the decrease in the statistical traveltime fluctuations as a function of the offset (the source–receiver distance). Our formulae are based on the Rytov approximation and the parabolic approximation for acoustic waves. The validity and the limits are established for both of these approximations in statistically anisotropic random media. Finally, theoretical inversion procedures are developed for the horizontal correlation structure of the velocity heterogeneities for the simplest case of a horizontal reflector. Synthetic seismograms are then computed (on particular realizations of random media) by simulating scalar wave propagation via finite difference algorithms. There is good agreement between the theoretical and experimental results.     

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.