Towards Regional tomography using the Spectral Element Method

Barbara Romanowicz, Aimin Cao, Federica Marone, Mark Panning, Paul Cupillard


We are developing an approach which relies on a cascade of increasingly accurate theoretical approximations for the computation of the seismic wavefield, with the goal to develop a model of regional structure for a subregion of Southeast Asia (longitude 75 to 150 degrees and latitude 0 to 45 degrees). The selected area is highly heterogeneous, but is well surrounded by earthquake sources and includes high quality broadband digital stations. In previous years, we developed preliminary models based on time domain inversion of long period seismograms, in the framework of normal mode theory: 1) a 3D model based on the Nonlinear Asymptotic Coupling Theory (NACT, Li and Romanowicz, 1995), which includes the consideration of 2D kernels in the vertical plane containing source and receiver. This model was developed for a larger region (longitude 30-150 degrees and latitude -10 to 60 degrees), starting from a global 3D model (Panning and Romanowicz, 2006); 2) a 3D model based on the “NBorn” approximation (Panning et al., 2008). This approach combines the 3D Born approximation with the path average approximation (PAVA) and allows us to accurately account for large accumulated phase delays on paths that sample large scale smooth anomalies. The resulting model has a horizontal resolution of about 200 km.

In parallel, a regional version of the Spectral Element Method (SEM) code, in spherical geometry, RegSEM, was completed (Cupillard, 2008). This code accepts a non-conformal grid, uses PML (Perfectly-Matched Layers) at the borders of the region, and includes general 3D anisotropy, Moho and surface topography, ocean bathymetry, attenuation, and ellipticity. Because each SEM run (i.e. for one event) is time consuming, we proposed to implement an approach in which the wavefields for several events are computed simultaneously. This approach was introduced by Capdeville et al. (2005) and tested on synthetic data at the global scale, but never applied to real data or at a regional scale.

Towards regional tomography using summed waveforms

In the past year, in collaboration with researchers at the Institut de Physique du Globe in Paris, we completed the codes and procedures that allow us to invert a collection of summed seismograms over the sub-region of Eurasia already considered for the NBORN inversion.

The computation of the forward wavefield using SEM is very accurate, but heavy computationally. It takes 2 hours to compute the wavefield of a single seismic event down to 60 sec (i.e., at long periods) on a modest computer cluster (32 cpu's). To speed up computations and develop the capability of performing many iterations and reaching higher frequencies, we have implemented the approach of Capdeville et al. (2005), in which the wavefield corresponding to many events is computed simultaneously. To do so, all the seismic sources considered are shifted in time to a common origin time, and the wavefield generated by this composite source is calculated once. The observed waveforms, at each station, thus consist of the summed observed seismograms for all events, each of them appropriately shifted in time, and are compared to the predicted synthetic wavefield. While these ``summed'' waveforms do not have a physical meaning, i.e. there are no identifiable seismic phases, they retain sufficient information to allow retrieval of 3D structure, even in the presence of realistic background noise, at least at long wavelengths. This approach can reduce computation time by one to several orders of magnitude, depending on the implementation. Capdeville et al. (2005) tested this approach in the framework of global tomography and a synthetic dataset. They showed that, even with realistic noise added to the synthetics, they can retrieve 3D upper mantle structure accurately from a single forward computation involving  50 events observed simultaneously at  150 stations worldwide.

In the adaptation of this method to our regional case, we have proceeded with the RegSEM code (rather than C-SEM) and have relaxed the restriction of “all events observed by all stations, replacing the single summation over all events by a series of runs, each of which considers the sum of all events observed at a given station. The number of RegSEM runs is then equal to the number of stations considered in the region, which, in general, is significantly smaller than the number of events. Contrary to the NACT and NBorn cases, for which we collected teleseismic data, RegSEM, as currently implemented, requires sources and receivers to be included in the region of study. We thus collected an entirely new dataset, consisting of 96 events with 6.0$<M_w<$7.0 observed simultaneously at 6 broadband GSN stations during the time period 1995 to 2007, with an epicentral distance between $5^o$ and $45^o$.

Figure 2.37 shows the region of study and the ray coverage for our initial summed event experiment. Here, we have relaxed the restriction of teleseismic observations. The epicentral distances range from 5 to $40^o$. In this case, we cannot perform any comparisons with the NACT or N-Born inversions, for which the theory breaks down at short distances; however, we can still perform comparisons with inversions obtained using the simpler Path Average (PAVA) approximation.

Figure 2.38 shows an example of comparison between ``observed'' versus RegSEM computed summed waveforms at station XAN, showing that, at these long periods, the starting model predicts the observed waveforms already quite well.

In order to reduce the computational time in the inverse step of the procedure, we also have opted, at least at this early stage, to compute partial derivatives approximately, using PAVA, in the inversion step of our procedure. We are assuming that, if our starting 3D model (the model derived using N-Born) is sufficiently accurate, these approximate partial derivatives will point us in the right direction and we will at the worst, need to compute a few more iterations to converge to the final model. Eventually, we can implement accurate numerical partial derivatives, whose calculation will also be faster, compared to conventional inversions, due to summation over events.

In order to test the summed event approach, we have first compared the results of an inversion obtained using, on the one hand side, the ``summed event'' seismograms at the 6 stations considered, and on the other, a ``standard'' approach, in which we have separately calculated the wavefield using RegSEM for all the events considered, and inverted the corresponding perturbations to the time domain waveforms, as we would do conventionally, using any of our previous inversion methods (PAVA, NACT, N-Born). The only difference in the latter is that now, the forward part of the modelling is computed using RegSEM, rather than mode perturbation theory. The synthetics have been computed down to 60s period, and the observed seismograms have been filtered accordingly. Figure 2.39 shows a comparison of the isotropic models obtained using the conventional approach (single event-station paths) versus the summed approach. The two RegSEM models are very similar indicating that the summed event approach holds promise. A striking feature of the RegSEM models is the fast velocity anomaly south of Korea shallow depth (30 km), which is not present in any of the models developed with asymptotic mode perturbation theory, which shows low velocities throughout the ocean basins at this depth. The fast velocity seen in the RegSEM models is likely an artifact, due to the inaccurate crustal model, which does not incorporate the thick sediments present in that region. Since RegSEM models the crust much more accurately than our usual inversions based on mode perturbation theory, these kind of details start to matter, even at periods as long as this test case.

Figure 2.37: Raypath coverage achieved in the ``summed seismogram'' experiment. The stations are indicated by triangles.
\epsfig{, bbllx=119,bblly=185,bburx=508


This research was supported by NNSA contract No. DE-FC52-04NA25543 (BAA).


Capdeville, Y., E. Chaljub, J. P. Vilotte, and 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-66, 2002.

Capdeville, Y., B. Romanowicz, and Y. Gung, Towards global earth tomography using the spectral element method: a technique based on source stacking, Geophys. J. Int., 162, 541-544, 2005.

Cupillard, P. , Simulation par la méthode des éléments spectraux des formes d'onde obtenues par corrélation de bruit sismique, These de Doctorat, Institut de Physique du Globe de Paris, 2008.

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

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.

Panning, M., Y. Capdeville and B. Romanowicz , Do first order 3D Born finite-frequency kernels improve modeling of surface waveforms?, Geophys. J. Int., in press, 2008.

Figure 2.38: Summed seismogram on the vertical component at station XAN. Black: observed trace, red: synthetic trace computed using RegSEM for the starting 3D model (NBORN) shown in Figure 2.37. Both observed and synthetic traces have been filtered with a cut-off period of 60 s.
\epsfig{, bbllx=0,bblly=301,bburx=624

Figure 2.39: Comparison of upper mantle models obtained after one iteration from the starting NBORN model, using waveform data on the path collection shown in Figure 2.38, RegSEM for the forward part of the computation, and PAVA for the calculation of partial derivatives. Left: conventional computation using individual source-station paths. Right: summed-event computation. The models differ in some details, which is not surprising given that only 6 stations have been used.
\epsfig{, bbllx=65,bblly=114,bburx=573

\epsfig{, bbllx=65,bblly=114,bburx=573

Berkeley Seismological Laboratory
215 McCone Hall, UC Berkeley, Berkeley, CA 94720-4760
Questions or comments? Send e-mail:
© 2007, The Regents of the University of California