A simple method for improving crustal corrections in waveform tomography

Vedran Lekic and Barbara Romanowicz

Introduction

Recordings of surface waves and overtones provide excellent constraints on the structure of the Earth's crust, upper mantle, and transition zone. This is because they provide good global coverage, and are sensitive to elastic and anelastic structure in both the crust and the mantle. Yet, in order to determine the seismic velocities and anisotropy in the mantle, we must separate the effects of the crust from those due to the sought-after mantle structure. Despite recent improvements in the global mapping of crustal structure (e.g. CRUST2.0: Bassin et al., 2000), the difficulties associated with accurately modeling the effects of the crust seismic on waves mean that these improvements do not automatically translate into better corrections for crustal effects.

Crustal corrections in long period waveform modeling were first applied in a linear fashion, by calculating the effects of perturbations in Mohorovicic depth and surface topography on the (eigen)frequencies of Earth's free oscillations (Woodhouse and Dziewonski, 1984). However, variations in crustal thickness are often large enough to produce non-linear effects on the eigenfrequencies. In order to account for this non-linearity, Montagner and Jobert (1988) proposed a two-step approach in which the eigenfunctions and eigenfrequencies are calculated exactly for a set of tectonic settings (thereby capturing the non-linear effects), while perturbations away from these canonical 1D profiles are handled using linear corrections. This approach has recently been implemented in full-waveform analyses (Marone and Romanowicz, 2007). However, when applied to higher modes and to high frequencies, these non-linear crustal corrections can be computationally very expensive.

Here, we present an alternative method for performing crustal corrections. Like the aforedescribed methods, we calculate exactly the eigenfunctions and eigenfrequencies for a set of tectonic settings, but instead of using these directly, we solve for scaling coefficients, which, when applied to standard linear crustal corrections, mimic the non-linear effect. The main advantage of this approach is that, once the correction factors have been calculated, it requires no additional computational costs aside from those associated with linear corrections. This allows it to be more easily applied to overtones and to higher frequencies than the standard quasi-non-linear approach.

Figure 2.29: Map showing geographical distribution of the 7 crustal types used in this study.

Modified linear crustal corrections

The effects of lateral heterogeneity, such as variations in crustal velocities and thicknesses, can be approximated by considering how this heterogeneity shifts the frequencies of Earth's free oscillations, compared to their frequencies in a spherically symmetric reference model such as PREM (Dziewonski and Anderson, 1981). For the case relevant to this study, in which only the radii of discontinuities in the Earth are perturbed by $\delta r_d$, for an isolated mode branch, local frequency shifts can be calculated in a linear fashion through the use of sensitivity kernels $H_k$, the expressions for which can be found in Woodhouse and Dahlen (1978). Note that these kernels are calculated for the spherically symmetric reference model:

$$\delta\omega_{k}^{2} = \sum_{d}r^{2}_{d}\delta r_{d}H^{d}_{k}.$$ (2.4)

In order to quantify the inadequacy of linear crustal corrections, we divide the Earth's surface into 7 regions with similar crustal thicknesses and ocean depths. We base this regionalization on Mohorovicic depth from CRUST2.0, with the first 6 regions characterized by Moho depth range of 10-25km, 25-40km, 40-50km, 50-60km, and $>$60km. The 7th region is introduced to capture the strong effect of a shallow ($<$ 2km) ocean layer that characterizes the continental shelves. Figure 2.29 shows the geographical extent of the 7 regions. For each region, we calculate an average radial profile of density ($\rho$) and shear ($V_S$) and compressional ($V_P$) wave velocity. Armed with a set of radial models that define 7 canonical crustal types, we proceed to calculate the frequencies $\omega_k$ of the fundamental modes for each model. These frequencies are then compared with those of PREM, and frequency shifts between PREM and each of the 7 regional models are calculated as

$$\delta\omega_k^{NL} = \omega_k^{(i)}-\omega_k^{PREM},$$ (2.5)

where the superscript $i$ is an index representing the frequency of mode k in the canonical crustal model i. Since these frequency shifts capture some of the non-linear effects of crustal structure, we identify them with a superscript $NL$. Figure 2.30 shows the $\delta\omega_k^{NL}$ plotted (solid lines) for spheroidal and toroidal fundamental modes.

Figure 2.30: Frequency shifts of the fundamental toroidal (grey) and spheroidal (black) modes with respect to PREM anisotropic due to differences in crustal structure between each of the canonical crustal types shown in Figure 2.29. Solid lines denote non-linear corrections ( $\delta\omega_k^{NL}$), dotted lines indicate linear corrections( $\delta\omega_k^{SL}$), and the dashed lines indicate linear corrections improved using the method outlined in this paper. Only Mohorovicic corrections are applied in the upper row, while corrections for both surface and Mohorovicic topography are required by the broader frequency range of the bottom row.

We can use the kernels $H_k^d$ that we obtained for the reference spherically symmetric model, in this case PREM, to predict the linearized effects of the canonical crustal structures on the normal mode frequencies. To do this, we only consider the differences in the radii of the discontinuities between each canonical crustal model and PREM, neglecting the differences in crustal velocities and density. This is an often used approximation of the true linear crustal effect, and is appropriate because crustal velocities have been shown to have minimal effect on long period waves (e.g. Stutzman and Montagner, 1994). Henceforth, which shall refer to the frequency shifts calculated in this standard, linear fashion as $\omega_k^{SL}$. The dotted lines in Figure 2.30 show the $\delta\omega_k^{SL}$ for each of the canonical crustal structures.

A comparison of the approximate terms $\delta\omega_k^{SL}$ with the $\delta\omega_{k}^{NL}$ calculated before (and displayed as solid lines) confirms that linear crustal corrections are inadequate, even at long periods. Therefore, we are interested in ways of correcting the $\delta\omega_k^{SL}$ so that they better track $\delta\omega_k^{NL}$. In order to accomplish this task, we are confronted with a crucial choice.

We must decide which term or terms in equation 2.4 to correct. Since $H_k^d$ needs to be calculated for each mode, correcting this term can be computationally expensive. This is what is done in the aforedescribed non-linear corrections. Correcting $\delta r_d$, on the other hand, does not increase computational costs. The gradual change with frequency of the differences between $\delta\omega_k^{NL}$ and $\delta\omega_k^{SL}$ change gives us hope that modifying $\delta r_d$ might significantly improve the accuracy of $\delta\omega_k^{SL}$.

We start the procedure by rewriting equation 2.4 in matrix notation, where we only consider N fundamental modes and identify perturbations relating to the Mohorovicic with a subscript $m$ and those pertaining to the surface with $t$:

$$\left( \begin{array}{c} \delta\omega_1^{SL} \\ \delta\omega_... ... \begin{array}{c} \delta r_m \\ \delta r_t \end{array} \right)$$ (2.6)

We attempt to improve standard linear corrections by introducing factors $c_{m,t}$ , calculated for each canonical crustal type and mode type, that are added to $\delta r_{m,t}$ before being multiplied by the kernel matrix (relabeled $\mathbf{H}$). Written in vector notation, we seek $\mathbf{c}$ that minimises:

$$\mathbf{w}-\mathbf{H}(\mathbf{\delta r_{m,t} + c_{m,t}}),$$ (2.7)

where the vector $\mathbf{w}$ contains the non-linear frequency shifts $\delta\omega_k^{NL}$. The least-squares solution to this minimisation problem is given by:

$$\mathbf{c_{m,t}} = (\mathbf{H'}\mathbf{H})^{-1}\mathbf{H'}(\mathbf{w - H \delta r_{m,t}}),$$ (2.8)

where the apostrophe indicates the transpose.

We could have introduced a multiplicative correction term, instead of the additive one described above. However, solving for such a term becomes unstable when the $\delta r_d$'s are small. Given that discontinuity topography is likely to vary both above and below its depth in the reference model, the accompanying zero-crossings of $\delta r_d$ might have adverse effects.

Because the non-linearity of crustal effects depends strongly on both crustal and mode type, we perform the minimisation in equation 2.8 for spheroidal and toroidal modes separately for each crustal type. Once the set of factors $\mathbf{c_{m,t}}$ appropriate for a given mode type are obtained, we modify the surface and Mohorovicic topography of CRUST2.0 at each point on the surface by the correction factor appropriate for the relevant crustal type (obtained from Figure 2.29). Therefore, the crustal type and correction factor information is fused into a single file that specifies a modified discontinuity topography for each mode type.

The dashed lines in Figure 2.30 show the frequency shifts predicted by our modified discontinuity radii. Henceforth, we label them $\delta\omega_k^{CL}$. The improvement in fit to $\delta\omega_k^{NL}$ is significant, and good for a large frequency range. When only long period waves ($T>60s$) are considered, excellent agreement between $\delta\omega_k^{CL}$ and $\delta\omega_k^{NL}$ can be achieved by only correcting the Mohorovicic topography.

Conclusions

We propose and validate a new method for improving linear crustal corrections. By considering a set of 7 crustal types, we quantify the inadequacy of standard linear corrections at accounting for the effects of the crust on the fundamental mode surface waves. Then, we improve the accuracy of linear corrections by introducing additive factors to the discontinuity topographies. Incorporating an additive correction factor to the discontinuity topography as opposed to the kernels, results in no additional computation costs, compared to standard linear corrections. The correction factors depend on the local crustal type, on the discontinuity considered, and on the reference model used for calculating the sensitivity kernels, as well as mode type.

Acknowledgements

This research was supported by National Science Foundation grant NSF/EAR-0308750 and an National Science Foundation Graduate Fellowship held by VL.

References

Bassin, C., Laske, G. and G. Masters, The Current Limits of Resolution for Surface Wave Tomography in North America. EOS Trans AGU 81, F897, 2000.

Dziewonski, A.M. and D.L. Anderson, Preliminary reference Earth model. Phys. Earth and Planetary Int., 25, 297-356, 1981.

Marone, F. and B. Romanowicz, Non-linear crustal corrections in high-resolution regional waveform seismic tomography. Geophys. J. Int., 101 22245-22272, 2007.

Montagner, J.-P. and N. Jobert, Vectorial tomography, II, Application to the Indian Ocean. Geophys. J. R. Astron. Soc. 94, 309-344, 1988.

Stutzmann, E. and J.P. Montagner, Tomography of the transition zone from the inversion of higher-mode surface waves. Phys. Earth and Planetary Int., 86, 99-115, 1994.

Woodhouse, J.H. and F.A., The effect of a general aspherical perturbation on the free oscillations of the earth, Geophys. J. Int., 53, 335-354, 1978.

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, 5953-5986, 1984.