Efficient Computation of NACT Seismograms: Toward Application in Imagining Upper Mantle Discontinuities

Zhao Zheng and Barbara Romanowicz

Introduction

The Berkeley Global Seismology group has been developing global mantle tomographic models (Li and Romanowicz, 1996; Mégnin and Romanowicz, 2000; Panning and Romanowicz, 2006) by inverting a dataset of surface waveforms down to 60 sec and long-period body waveforms down to 32 sec. A normal mode coupling method known as the Non-linear Asymptotic Coupling Theory (NACT; Li and Romanowicz, 1995) has been successfully applied to compute 3D synthetic seismogram and sensitivity kernels. It is an approximate approach which collapses sensitivity over the whole volume onto the vertical plane containing source and receiver. By taking into account the depth variation of sensitivity, NACT is able to bring out the ray character of body waves as well as the finite-frequency behavior of sensitivity kernels. On the other hand, as an approximate approach, it is computationally much faster than purely numerical methods such as the Spectral Element Method (SEM). The frequency band ($\geq$ 32 sec) of the current tomographic dataset is too low to retrieve fine scale structures of upper mantle discontinuities, which place vital constraints on the temperature, composition and dynamics of the mantle. In order to retrieve them, it is necessary to incorporate phases such as the SS precursors. Most SS precursor studies so far measure SS-SdS differential travel time and translate that to discontinuity depth, therefore suffering from the tradeoff between volumetric perturbation and discontinuity topography (for a review, see Deuss, 2009). Modeling the 3D waveforms of precursors will help to resolve this tradeoff. However, these precursor phases are typically peaked at a much shorter period ($\sim$10 sec), therefore presenting a great computational challenge to the current mode coupling scheme of NACT, in which the computational cost grows as $f_{max}^4$ with $f_{max}$ being the cutoff frequency.

An Efficient NACT Formalism

In the NACT theory, the perturbed seismogram consists of an along-branch coupling term, which is computed under the well-known PAVA approximation (Woodhouse and Dziewonski, 1984), and an across-branch coupling term, which is computed under the linear Born approximation. In the classical formalism, the Born part is obtained by a double summation over all pairs of coupling modes (Li and Romanowicz, 1995):

$$\delta u(t)=\sum_{k}\sum_{k'} { A_{kk'} \frac{exp(i\omega_k t)-exp(i\omega_{k'} t)}{\omega_k^2-\omega_{k'}^2} }$$ (19.1)

where $k=(n,l)$ and $k'=(n',l')$ are indices for normal modes ($n$ the radial and $l$ the angular order), and

$$A_{kk'}=\sum_{m}\sum_{m'} R_k^m H_{kk'}^{mm'} S_{k'}^{m'}$$ (19.2)

where $R_k^m$ and $S_{k'}^{m'}$ ($m$ being the azimuthal order of normal mode) are the receiver and source terms as defined in Woodhouse and Girnius (1982), and the coupling matrix $H_{kk'}$ due to perturbation in the elastic tensor ${\mathbf C}$ is:

$$H_{kk'}^{mm'}=\int_V {\mathbf e}_k^{m*}(\mathbf{r};\mathbf... ...hbf{r})} : {\mathbf e}_{k'}^{m'}(\mathbf{r};\mathbf{r}_S) dV$$ (19.3)

with ${\mathbf e}_k^m$ being the strain of mode $(k,m)$, subscript $R$ and $S$ for the location of seismic receiver and source, ``$*$'' denoting complex conjugate, and ``$:$'' double dot product. In practice, the summation is truncated below some cutoff frequency $f_{max}$. Since the number of modes with eigenfrequencies below a certain cutoff frequency is roughly proportional to $f_{max}^2$, the cost of computing one seismogram grows with frequency as $f_{max}^4$. For multiple sources and receivers, the computational cost goes $(N_S*N_R)*f_{max}^4$, where $N_S$ and $N_R$ are the number of sources and receivers, respectively. Here we present a more efficient formalism. From the physical inspiration that the Born approximation is a single-scattering approximation, i.e., the seismic wavefield due to a heterogeneity scatterer is the convolution of a source-to-scatterer term and a scatterer-to-receiver term, we see it is possible to separate the double summation over mode pairs to two single summations ( Capdeville, 2005). In particular, the Born seismogram can be expressed in the frequency domain as (Tanimoto, 1984):

$$\delta u(\omega)=\sum_{k}\sum_{k'} { A_{kk'} \frac{1}{i\omega(\omega^2-\omega_k^2)(\omega^2-\omega_{k'}^2)} }$$ (19.4)

Upon substituting (2) in (4) and separating the terms that are indexed by $(k,m)$ and those by $(k',m')$, one finds

$$\delta u(\omega) = \int_V { {\mathbf R}(\mathbf{r},\omega)... ...\mathbf C}(\mathbf{r}) : {\mathbf S}(\mathbf{r},\omega) } dV$$ (19.5)

with

$${\mathbf R}(\mathbf{r},\omega;\mathbf{r}_R) = \sum_k (\sum_m R_k^m {\mathbf e}_k^{m*}) \frac{1}{\omega^2-\omega_k^2}$$ (19.6)

$${\mathbf S}(\mathbf{r},\omega;\mathbf{r}_S) = \sum_{k'} (\sum_{m'} S_{k'}^{m'} {\mathbf e}_{k'}^{m'}) \frac{1}{i\omega(\omega^2-\omega_{k'}^2)}$$ (19.7)

accounting for the scatterer-to-receiver and the source-to-scatterer contributions, respectively. Returning to the time domain, one can write

$$\delta u(t) = \int_V { {\mathbf R}(\mathbf{r},t) : \delta{\mathbf C}(\mathbf{r}) \star : {\mathbf S}(\mathbf{r},t) } dV$$ (19.8)

with

$${\mathbf R}(\mathbf{r},t;\mathbf{r}_R) =\sum_k (\sum_m R_k^m {\mathbf e}_k^{m*}) \frac{\sin(\omega_k t)}{\omega_k}$$ (19.9)

$${\mathbf S}(\mathbf{r},t;\mathbf{r}_S) =\sum_{k'} (\sum_{m'} S_{k'}^{m'} {\mathbf e}_{k'}^{m'}) \frac{1-\cos(\omega_{k'} t)}{\omega_{k'}^2}$$ (19.10)

and ``$\star$'' in (8) denoting temporal convolution. As commented by Capdeville (2005), this formalism is similar to the adjoint method (e.g. Tarantola, 1984). Since the summations in (9) and (10) can be performed simultaneously, the cost of computing one seismogram now grows as $f_{max}^2$ as opposed to $f_{max}^4$ in the original NACT. For multiple paths, the cost is $(N_S+N_R)*f_{max}^2$, as compared to $(N_S*N_R)*f_{max}^4$. Asymptotic approximation can still be applied to collapse the integral over whole volume in (8) onto the great circle plane, as did the original NACT.

Numerical Validation

A numerical test is designed to validate the new formalism. Synthetic seismograms are computed for a path geometry shown in Figure 2.32 using the old and new formalisms, respectively. The result shows that the two synthetics agree very well in both phase and amplitude.

Figure 2.32: (a) Path geometry of the numerical test. A cylindrical anomaly (blue) between 0 and 1000 km depth with 2% faster $V_S$ than PREM is located halfway between the source and the receiver. (b) Transverse component synthetic seismogram computed up to 60 sec with the original NACT (black line) versus the new formalism (blue dashed line). Major seismic phases are labeled.

Table 2.2 lists the CPU time of computing a single synthetic seismogram with the two formalisms, respectively, for several cutoff frequencies. Time cost of the new formalism scales roughly as $f_{max}^2$, whereas that of the original NACT scales as $f_{max}^3$ (rather than the theoretically $f_{max}^4$, because in practice the mode coupling was restricted to some fixed coupling width $s$ such that $\vert l'- l\vert \leq s$, where $l$ is the angular order of modes). Therefore, the new formalism is computationally preferable when approaching higher frequencies. It is even more advantageous in the case of many sources and receivers, which is the reality of seismic tomography.

Table 2.2: Time cost for computing one synthetic seismogram on a 1 GHz single CPU.

$f_{max}$ (period)	original NACT	new formalism
60 sec	6.9 min	22 min
30 sec	56 min	72 min
10 sec	25 hr	10 hr

References

Capdeville, Y., An efficient Born normal mode method to compute sensitivity kernels and synthetic seismograms in the Earth, Geophys. J. Int., 163, 639-646, 2005. Deuss, A., Global observations of mantle discontinuities using SS and PP precursors, Surv. Geophys., 30, 301-326, 2005. Li, X.D. and B. Romanowicz, Comparison of global waveform inversions with and without considering cross-branch modal coupling, Geophys. J. Int., 121, 695-709, 1995. Li, X.D. and B. Romanowicz, Global mantle shear-velocity model developed using nonlinear asymptotic coupling theory, J. Geophys. Res., 101(B10), 22,245-22,272, 1996. Mégnin, C. and 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(3), 709-728, 2000. Panning, M. and B. Romanowicz, A three dimensional radially anisotropic model of shear velocity in the whole mantle, Geophys. J. Int., 167, 361-379, 2006. Tanimoto, T., A simple derivation of the formula to calculate synthetic long-period seismograms in a heterogeneous earth by normal mode summation, Geophys. J. Int., 77, 275-278, 1984. Tarantola, A., Inversion of seismic reflection data in the acoustic approximation, Geophysics, 49, 1259-1266, 1984. Woodhouse, J.H. and A.M. Dziewonski, Mapping the Upper Mantle: Three-Dimensional Modeling of Earth Structure by Inversion of Seismic Waveforms, J. Geophys. Res., 89(B7), 5953-5986, 1984. Woodhouse, J.H. and T.P. Girnius, Surface waves and free oscillations in a regionalized earth model, Geophys. J. R. Astr. Soc., 68(3), 653-673, 1982.