Full waveform tomography with the spectral element method

Map of lateral variations in the velocity of seismic shear waves (S-waves) at a depth of 125 km. Warm colors denote slower wavespeeds, and are thought to indicate warmer temperatures and/or the presence of melt or water. Cool colors denote faster wavespeeds, and indicate cool regions, such as the mantle beneath the ancient, stable portions of continents.

Mapping the elastic and anelastic structure of the Earth's mantle is crucial for understanding the temperature, composition and dynamics of our planet. In the past quarter century, global tomography based on ray theory and first-order perturbation methods has imaged long-wavelength elastic velocity heterogeneities of the Earth's mantle. However, the approximate techniques upon which global tomographers have traditionally relied become inadequate when dealing with crustal structure, as well as short-wavelength or large amplitude mantle heterogeneity. The spectral element method, on the other hand, permits accurate calculation of wave propagation through highly heterogeneous structures, and is computationally economical when coupled with a normal mode solution and applied to a restricted region of the earth such as the upper mantle (SEM: Capdeville et al., 2003). Importantly, SEM allows a dramatic improvement in accounting for the effects of crustal structure.

 

We developed and applied a new hybrid method of tomography, which allows us to leverage the accuracy of SEM to model fundamental and higher-mode long period (>60s) waveforms. We then present the first global model of upper mantle velocity and radial anisotropy developed using SEM. Our model, SEMum, confirms that the long-wavelength mantle structure imaged using approximate semi-analytic techniques is robust and representative of the Earth's true structure. Furthermore, it reveals structures in the upper mantle that were not clearly seen in previous global tomographic models, providing new constraints on the temperature, composition as well as flow in the mantle. We show that SEMum favorably compares to and rivals the resolving power of continental-scale studies. This new hybrid approach to tomography can be applied to a larger and higher-frequency dataset in order to gain new insights into the structure of the lower mantle and more robustly map seismic structure at the regional and smaller scales.

 

SEMum is now available for download by clicking here.

 

We developed SEMum, our high resolution model of upper mantle structure, using a fully numerical wave propagation code C-SEM (1) which minimizes forward-modeling errors (2) and is capable of accurately representing both the scattering and (de)focusing of seismic waves by elastic heterogeneity, and the effects of the oceans, topography/bathymetry, ellipticity, gravity, rotation and anelasticity. We optimized data utilization through the use of full-waveform modeling, and minimized crustal contamination by including both long period waveforms and higher-frequency group velocity dispersion maps. We keep computational costs reasonable by calculating partial derivatives that relate structure perturbations to waveform perturbations using approximate, 2D finite-frequency kernels based on normal modes (3).

 

Data – Our dataset includes minor- and major-arc fundamental mode and overtone wavepackets from 200 earthquakes recorded at more than 100 broadband stations around the globe.  In order to minimize the contamination from source complexity, all earthquakes have moment magnitude <7.0 (moment tensors and centroid locations are from the Harvard Centroid Moment Tensor Project). The complete seismic accelerograms for all three components of motion are bandpass filtered with corners at 80 and 250 s, and are divided into wavepackets, which separate, in the time domain, the larger fundamental mode surface waves from the smaller overtones (4). Finally, an automatic, but user-reviewed, picking scheme ensures that only well-recorded data are included (5). In order to improve constraints on crustal structure, this waveform dataset is supplemented with Love and Rayleigh group velocity dispersion maps between 25s and 150s period, provided by M. Ritzwoller.

 

Parameterization – In SEMum, we model the Earth as a transversely isotropic medium, which is described by five elastic parameters: Voigt average Vp, Vs, and anisotropic parameters  ξ = VSH²/VSV² φ = VPV²/VPH², and η. Because long period Love and Rayleigh waves are largely insensitive to Vp, we scale deviations of Vp to Vs, and of j and h to those of ξ (5).

 

In depth, SEMum is expressed on 21 cubic splines defined in Megnin and Romanowicz (7), whose knot locations are at radii: 3480, 3600, 3775, 4000, 4275, 4550, 4850, 5150, 5375, 5575, 5750, 5900, 6050, 6100, 6150, 6200, 6250, 6300, 6346, 6361km and the surface.  Because our dataset is mostly insensitive to lower mantle structure, we keep structure below 800 km fixed to that in SAW24B16 (6). Lateral variations of Vs and x are parameterized with 2562 and 642 spherical splines, respectively (equivalent to degree 48 and 24 spherical harmonic expansions).  The local spline parameterization reduces both aliasing (7) and contamination from distant structure.

 

In order to avoid biasing SEMum toward previous 3D models of Vs structure, we choose to start with a 1D model. Our starting model is identical to PREM (8) below 400 km depth; at shallower depths, its isotropic profile is obtained by inversion of long period waveforms from a physical reference model of Cammarano et al. (9). A starting 1D model of x is found through a grid search to fit group velocity dispersion maps at periods shorter than 80s.

 

Forward-modeling and Inversion Scheme – Because SEM-based calculation of the forward wavefield is very computationally expensive, we use approximate NACT partial derivatives relating model perturbations to waveform perturbations. The use of SEM-based adjoint methods (10) would result in a many-fold increase in our computational costs. NACT kernels, while approximate, do account for finite-frequency effects in the vertical plane defined by the great circle path from source to receiver, and are more accurate (3) in representing the sensitivity of overtones compared to traditional means of calculating this sensitivity (e.g. PAVA: 11). Furthermore, the kernels need only be accurate enough to give the correct sign of the partial derivatives when the kernels for all available waveforms are summed at a given model parameter. For the group velocity dataset, both the forward modeling and sensitivity kernels are calculated by perturbation theory. We minimize the effects of non-linearity in the sensitivity of group velocity to structure by calculating 5 representative profiles of structure at each iteration, so that the perturbations between the closest of these profiles and the structure beneath any given point on the Earth is within the linear regime.

 

Because waveform tomography is highly non-linear, an iterative inversion procedure is necessary (12). The scheme is initialized using our starting 1D model and CRUST2 crustal velocities and Mohorovicic topography. Misfits between observed waveforms and synthetics calculated using SEM are calculated point-by-point in the time domain. In order to avoid cycle-skipping, only waveforms that are sufficiently similar to the synthetics are used in each iteration. As a result, the model parameterization at the start only allows long wavelength structure to be imaged; with successive iterations, the synthetics better predict the observed waveforms, and the effective size of the dataset increases allowing for the parameterization to be refined accordingly. The fact that fits improve even for waveforms that are not included in the inversion at a given iteration independently confirms the validity of the inversion scheme and justifies the use of the NACT kernels.

 

Inaccurate calculation of the effects of crustal structure on waveforms has the potential to contaminate the retrieved mantle isotropic and especially anisotropic structure. Therefore, accurate treatment of the effects of the crust is crucial for developing a robust and reliable mantle model. Unlike wave propagation calculations carried out using perturbation theory, the SEM can accurately represent the effects of the crust on the seismic waveforms (13). The presence of a thin, low velocity crust in oceanic areas imposes a small timestep in the SEM computation, slowing it down considerably. After 4 iterations of our inversion with CRUST2 as a crustal model, we replaced CRUST2 with a smooth crustal model which allowed us to both dramatically speed up SEM calculations and improve the waveform and group-velocity fits compared to CRUST2. The smooth starting crustal model is obtained by generating 21,000 models of smooth crustal structure and choosing among these the model that best fits the Rayleigh and Love group velocity dispersion beneath each point on the surface.

 

A total of 10 iterations are carried out before convergence is achieved. The SEMum model thus obtained provides >70% variance reduction with respect to the starting model for fundamental mode surface waves on all three components, and ~50% variance reduction for overtone waveforms. SEMum also explains ~60% of the variance in the group velocity dataset. An analysis of the phase-alignment and amplitude fits confirms that while phase-alignment is the dominant source of the variance reductions provided by SEMum, amplitude fits are also improved.

 

 

References

1. Y. Capdeville, E. Chaljub, J. P. Vilotte, J. P. Montagner, Coupling the spectral element method with a modal solution for elastic wave propagation in global earth models. Geophys. J. Int. 152, 34–67 (2003).

2. D. Komatitsch, J. Vilotte, The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures, Bull. Seismol. Soc. Am. 88(2), 368-392 (1998).

3. X-D. Li, B. Romanowicz, Comparison of global waveform inversions with and without considering cross-branch modal coupling. Geophys. J. Int. 121, 695-709 (1995).

4. X-D. Li, B. Romanowicz, Global mantle shear velocity model developed using non-linear asymptotic coupling theory, J. Geophys. Res. 101, 22,245-22,272 (1996).

5. see Appendices of M. Panning, B. Romanowicz, Romanowicz, A three-dimensional radially anisotropic model of shear velocity in the whole mantle, Geophys. J. Int. 167, 361-379 (2006).

6. Mégnin, B. Romanowicz, The three-dimensional shear velocity structure of the mantle from the inversion of body, surface and higher-mode waveforms. Geophys. J. Int. 143, 709–728 (2000).

7. L. Chiao, B. Kuo, Multiscale seismic tomography, Geophys. J. Int. 145(2), 517-527 (2001).

8. A. Dziewonski, D. Anderson, Preliminary reference earth model, Phys. Earth Planet. Inter. 25, 297-356 (1981).

9. F. Cammarano, A. Deuss, S. Goes, D. Giardini, One-dimensional physical reference models for the upper mantle and transition zone: Combining seismic and mineral physics constraints, J. Geophys. Res. 110, B01306, doi:10.1029/2004JB003272 (2005).

10. A. Tarantola, Inversion of seismic reflection data in the acoustic approximation, Geophysics 49(8), 1259-1266 (1984).

11. J. Woodhouse, A. Dziewonski, Mapping the upper mantle: Three dimensional modeling of earth structure by inversion of seismic waveforms, J. Geophys. Res. 89, 5953-5986 (1984).

12. A. Tarantola, B. Valette, Generalized nonlinear inverse problems solved using the least squares criterion, Rev. Geophys. Space Phys. 20, 219-232 (1982).

13. D. Komatitsch, J. Tromp, Spectral-element simulations of global seismi wave propagation – ii. Three-dimensional models, oceans, rotation and self-gravitation, Geophys. J. Int. 150, 303-318 (2002).