Global Waveform Tomography with Spectral Element Method

Yann Capdeville, Barbara Romanowicz, Yuang-Cheng Gung

Introduction and Research
Objectives

Because seismogram waveforms contain much more information on the earth structure than body wave time arrivals or surface wave phase velocities, inversion of complete time-domain seismograms should allow much better resolution in global tomography. In order to achieve this, accurate methods for the calculation of forward propagation of waves in a 3D earth need to be utilized, which presents theoretical as well as computational challenges.

In the past 8 years, we have developed several global 3D S velocity models based on long period waveform data, and a normal mode asymptotic perturbation formalism (NACT, Li and Romanowicz, 1996). While this approach is relatively accessible from the computational point of view, it relies on the assumption of smooth heterogeneity in a single scattering framework. Recently, the introduction of the spectral element method (SEM) has been a major step forward in the computation of seismic waveforms in a global 3D earth with no restrictions on the size of heterogeneities (Chaljub, 2003). While this method is computationally heavy when the goal is to compute large numbers of seismograms down to typical body wave periods (1-10 sec), it is much more accessible when restricted to low frequencies (T$>$150sec). When coupled with normal modes (e.g. Capdeville et al., 2000), the numerical computation can be restricted to a spherical shell within which heterogeneity is considered, further reducing the computational time.

Here, we present a tomographic method based on the non linear least square inversion of time domain seismograms using the coupled method of spectral elements and modal solution. SEM/modes are used for both the forward modeling and to compute partial derivatives. The parameterization of the model is also based on the spectral element mesh, the "cubed sphere", which leads to a 3D local polynomial parameterization. This parameterization, combined with the excellent earth coverage resulting from the full 3D theory used for the forward modeling, leads to a very stable inversion scheme. Synthetic tests show that, with a limited number of events (between 50 and 100), using long period records ($>$150s) and representative background seismic noise, it is possible to recover both amplitude and phase of the earth model with high accuracy.

Method

Our aim is to find the ``best'' Earth model that explain our seismic data set. Let assume we wish the solve the inverse problem with a classical least square inversion with a complete modeling theory that is the Spectral Element Method (SEM) applied to the wave equation. Classical inversion processes require to compute the partial derivative matrix

$$\mathbf{A}=\left[\frac{\partial\mathbf{f}(\mathbf{x})}{\partial\mathbf{x}}\right]_{\mathbf{x}} .$$ (32.1)

where $\mathbf{x}$ the set of parameters which describe the model, and $\mathbf{f}$ is the ``function'' that produce a synthetic data set $\mathbf{d}$ for a given set of model parameters $\mathbf{x}$.

In most of the classical tomographic methods, the forward problem is solved using the Born approximation within the normal modes framework which leads to a linear relation between the set of parameters and the synthetic data, and this relation matrix is $\mathbf{A}$. Our case is different since there is no way to compute directly the partial derivative matrix or kernel with the SEM but using brute force with a finite differences formula:

$$\left[\frac{\partial \mathbf{f}(\mathbf{x})}{\partial \mathb... ...ta\mathbf{x})-\mathbf{f}(\mathbf{x}_0)}{\delta \mathbf{x} } .$$ (32.2)

Therefore, in order to compute the partial derivative matrix, one need to compute the whole data set, that is as many runs as the number of sources, for each parameter of the model. This will be obviously the most expensive part of the inversion. An estimation of the computing time shows that computing the partial derivative matrix might takes years, even on large computers.

Here we used a scheme that reduce the computation by a factor equal to the number of sources and make the process possible within a reasonable amount of time.

Figure 32.1: Distribution of the sources and stations used in this experiment.

We now show a synthetic inversion to test the scheme. The ``data'' are produced with the SEM in a known model that we try to retrieve. The source-receiver distribution is shown on Figure 32.1. The parameterization used here to describe the model is based of the spectral element mesh. We use here a mesh roughly equivalent to a degree 8 in spherical harmonics as shown on Figure 32.3 . We show results of the 2 first iterations of the inversion on Figure 32.2. The model is represented in a linear way: the velocity contrast are shown as a function of the parameter number. It is sorted such that the left part is the lower mantle and the right part is the upper mantle. The residual between the wanted model an the obtained model shows an excellent agreement.

Figure 32.2: Mesh used for the parameterization of the experiment

Figure 32.3: Synthetic inversion model results of the two first iterations in a linear representation (see text). The residual between the wanted model an the obtained model shows an excellent agreement

A description of the coupled method with illustrations can be found on http://seismo.berkeley.edu/~yann

Acknowledgements

Most of the computation were made using the computational resources of the NERSC, especially the IBM SP, under repo mp342.

References

Capdeville, Y., B. Romanowicz and A. To (2003). Coupling Spectral Elements and Modes in a spherical earth: an extension to the ``sandwich'' case. Geophys. J. Int. Vol 154, pp 44-57
Chaljub, E. Capdeville, Y. and Vilotte, J.P. (2003) Solving elastodynamics in a solid heterogeneous 3D-sphere: a parallel spectral element approximation on non-conforming grids. J. Comp Phy. Vol. 183 pp. 457-491
Capdeville Y., E. Chaljub, J.P. Vilotte and J.P. Montagner (2003). Coupling the Spectral Element Method with a modal solution for Elastic Wave Propagation in Global Earth Models. Geophys. J. Int., Vol 152, pp 34-66
Capdeville Y., C. Larmat, J.P. Vilotte and J.P. Montagner (2002) Numerical simulation of the scattering induced by a localized plume-like anomaly using a coupled spectral element and modal solution method. Geoph. Res. Lett. VOL. 29, No. 9, 10.1029/2001GL013747
Capdeville Y., E. Stutzmann and J. P. Montagner (2000). Effect of a plume on long period surface waves computed with normal modes coupling. Phys. Earth Planet. Inter., Vol 119 pp. 57-74.
Li, X. B. and Romanowicz, B. (1996) . Global mantle shear velocity model developed using nonlinear asymptotic coupling theory J. Geophys. Res. Vol 101, pp 11,245-22,271.