Coupled Method of Spectral Element and Modal solution for Wave Propagation in Global Earth Model

Contents

• Introduction
• Principle of the method
• Technical Aspect
  • Spectral Element Method
  • Coupling
  • Code portability
  
  
• Validations
  • Wave propagation in a homogeneous sphere
  • Wave propagation in PREM
  
  
• Some 3D runs
  • SAW12D
  • Plume
  
  
• Bibliography

Introduction

The coupled method of Spectral Elements and Modal Solution has been the main work of my Ph.D. in Global Seismology and an important part of my work as a Miller fellow at the Seismological Lab. of UC Berkeley. This work is related to the ones of Dimitri Komatitsch and done with Emmanuel Chaljub on the Spectral Element Method applied to the wave equation developed at the IPG de Paris during their Ph.D. with Jean Pierre Vilotte (Dimitri is now at Pau and Emmanuel at Grenoble)

This quick and dirty introduction to the method is out of date since 2001

Solving the wave equation in 3D Global Earth Models is a challenging problem. Usually, in global seismology, it is solve by at least assuming that velocity contrasts are small compare to a spherically symmetric reference Earth model. Unfortunately, in order to go further in our understanding of the Earth interior, new numerical methods that are not based on this assumption are required for at least two reasons:

• the assumption of small velocity contrast is not valid in at least in the crust and the in D'' layer above the Core Mantle Boundary.
• The actual method based upon this approximation are not numerically efficient for higher frequencies than 50 s which limitate considerably the range of application of these methods.

Recently a new numerical method, the Spectral Element Method (SEM), introduce by Komatitsch and Vilotte (1998) for regional applications and by Chaljub & all (2003) for its applications to spherical geometry, has shown to be very accurate and efficient to solve the wave equation in the general case. This method is efficient for surface waves, body waves, Stoneley waves ... in media with large velocity contrasts, topography, liquide, gravity, anisotropy, anelasticity, etc. Unfortunately, its numerical cost is still to high to access the frequency range necessary for a large number of applications. The main idea of the method presented here is to allow higher frequency simulation by limiting the use of the SEM to region of the Earth Model that has to be 3D for a particular study and using an other method, a Modal Solution (MS), where the model can be assumed to be spherically symmetric.

Principle of the method

The numerical method (Capdeville, 2000; Capdeville & all, 2003) is based on the partition of the earth in two domains : an outer shell, in which 3D heterogeneities are confined, and an inner spherically symmetric sphere, see Fig. 1.

Figure 1: The Earth model is partitioned in two parts, an external shell SEM and an internal sphere MS connected by a spherical boundary. Lateral heterogeneities are confined within the outer shell.

The two domains are connected through a spherical interface ($\Gamma$ ). Depending on the problem, the outer heterogeneous shell can be either the whole mantle or a portion of the upper mantle or the crust. In the outer shell, the solution is sought in terms of the SEM, based on a high order variational formulation in space and a second order explicit time scheme (Komatitsch and Vilotte, 1998). Using a central projection transformation (Chaljub, 2000; Chaljub & all, 2003) the outer shell can be discretized almost uniformly by hexahedra elements avoiding the pole problems associated with the spherical coordinates system (Fig. 2).

Figure: Left: a splitted view of the six regions produced by the gnomic projection which map the six faces of a cube inscribed inside the sphere to the surface of the sphere. Right: a gathered view of the six regions of the sphere. Each region has its own coordinates system and the coordinate transformations can be calculated analytically.
(click on image to zoom it)

High resolution can be achieved, for realistic Earth models, using non-conforming mesh refinements (Chaljub & all, 2003) (Fig. 3).

Figure: Radial view of one region of the SEM mesh. Two different mesh are used in and connected through the non-conforming interface $\Sigma _{12}$ . $\Gamma$ is the coupling interface between MS and SEM, and is the MS area.

In the inner sphere, the solution is based upon a modal method in frequency after a spherical harmonic expansion. The modal solution allows to build explicitly the interface operator (known as a Dirichlet to Neumann (DtN) operator) that relates the traction along the coupling interface in terms of the displacement along the same interface. Using the continuity conditions of both traction and displacement, this operator fully characterizes the coupling between the two methods and provides, for the outer shell, the traction along the interface as a response of the inner sphere to a surface displacement distribution. In the SEM, this is introduced as a dynamic boundary condition and significantly reduces the computational cost. However, the construction of the interface operator is only explicit in the frequency and spherical harmonics domains. The back transformation in the space and time domain requires special attention and an asymptotic regularization (Capdeville & all, 2003)

Technical Aspect

The program is mainly written in Fortran90 with some Fortran77 and C. The MPI library is used for the parallelization.

Spectral Element Method

For Large problem, the critical part is the SEM which uses the largest amount of memory and computing time. The number of double precision variables required is roughly where the number of scalar fields used by the SEM (currently 36) the number of elements (proportional the frequency ) and the polynomial degree (usually 8). The memory becomes rapidly a problem as the frequency increase (each time the maximum frequency is multiply by 2, the memory needed is multiplied by 8). To give an idea, for the run in PREM presented here, the memory required for the SEM part is about 6 G-bytes.

As for all finite element method, the SEM can be well parallelized. All computations are performed locally on processors and only the assembly phase require communication between processors. Each one of the 6 regions are distributed into domains of elements. The number of domains is equal to number of processors and each domain has the same number of elements. This distribution is only performed in two dimensions in the horizontal directions. The coloring (numbering) of domains is performed such as no processor has to communicate to himself. An example with 16 processors is given Fig. 4. Practically, due to this distribution, only number of processors that are a power of two or twice a power of two are allowed.

Figure: Example of disctribution of the spectral element mesh of 16 processors. Each number represent a processor number.
(click on image to zoom it)

The vertical assembly do not required any communication, but the horizontal one requiered 6 communication phases described on Fig. 5.

Figure: The 6 communication phases requiered for the assembly phase of the SEM.
(click on image to zoom it)

Coupling

The coupling part is usually small in term of memory and computation time face to the SEM. Nevertheless, in some cases (very small case, or a very small area of spectral elements compare to the volume of the Earth), the two parts can be equivalent.

The coupling is performed in the spherical harmonic domain (spectral in two dimensions, (l,m) space) which require two FLT (Fast Legendre Transform) per SEM time step. The FLT integrations are performed on a different mesh that the SEM mesh which require two interpolations (one forward and one backward). These interpolations can be performed locally on processors and require a neglectable number of communications. They nevertheless require to store the interpolation matrixes which, depending on the case, can be significant (up to 25% of the whole memory in bad configuration). The FLT integrations use 1D FFT in one direction and matrix products in the other one. Each direction can be performed locally on processors, but a 'ALL TO ALL' communication type is needed between the two directions. As this step is only 2D (the coupling is perform on a 2D interface), the time spend in this communication phase is neglectable. The FLT process is summarized on Fig. 6.

Figure: Parallel FLT on 4 processors.
(click on image to zoom it)

In the spectral domain, the computation of traction using the DtN operator can be fully parallelize without any communication. The DtN operator is not localized in time and all previous time steps to the current one are required. To bypass this difficulty, only a small part of the time steps are stored in the spectral domain, the rest is stored in the frequency domain where only a small amount of eigenfrequencies are allowed. Nevertheless, in bad cases (2 DtN, very small area of spectral elements) the cost in memory can be up to 35% of the whole problem.

Code portability

Due to the choice of the MPI library and Fortran90 language, the code is relativelly easy to port to a new parallel machine.

• The SEM and the Coupled Method have been developed in Paris at the DMPN on a Digital-Alpha EV6 cluster (32 processors) and on a Linux PC cluster (8 AMD Athlon 800 MHz processors);
• It runs on the CINES's SGI-CRAY ORIGIN 2000;
• It has recently ported on the NERSC's IBM SP;
• It also runs on a 4 processors SUN machine at Berkeley.

Validations

To validate the coupled method, synthetic seismograms for a simple homogeneous sphere and the spherical earth model PREM are computed and compared with the solution obtained by the summation of the free oscillations in order to test the accuracy of the method. In a homogeneous sphere, normal modes radial eigenfunctions are known analytically and eigenfrequencies are known up to the computer accuracy. In a spherically symmetric earth model, normal modes are known with a very good precision. Therefore, in both cases, the normal mode solution is a very accurate reference solution which allows for the validation of the coupled method.

The homogeneous sphere test only incorporates a solid-solid interface,while for the more realistic PREM test, the coupling interface is set to be the CMB and includes a solid-liquid coupling.

Wave propagation in a homogeneous sphere

We consider here a homogeneous solid sphere with a radius km, a density kg/m , a P wave speed km/s and a S wave speed km/s. The position of the solid-solid coupling interface is set to km, and therefore the overall thickness of the outer spherical shell is km. Each of the 6 regions of are discretized with elements horizontally and 2 elements vertically (see Fig. 7), with a total number of elements of .

Figure: The spectral element mesh used for the homogeneous sphere model. The regions 1 and 4 have been removed for visibility.
(click on image to zoom it)

In each element, the polynomial degree is , which gives a total number of grid points in of 559872. For such a mesh, the minimal distance between two grid points is km, which, for a Courant number of , leads to a maximal time step of s. In the numerical experiment, the time step was set to s, which is quite conservative.

For the source, the lowest period of the wavelet time function has been set to s which, for a wave speed of km/s (the approximate Rayleigh wave group velocity in this model), insures less than two wavelengths per element. The wavelet time function, is a Ricker with a central frequency Hz and a corner frequency of Hz.

Snapshots of the wave propagation in are shown in crossection, Fig. 8, for a source at 1050 km depth. At s, the central time of the source time function is s, the P-wave begins to be absorbed by the DtN operator, without any visible spurious reflection. At s, the pP wave start to be absorbed by the DtN operator. At s the P-wave emerges from into . The X phase, a superposition of higher spheroidal modes, is strongly excited by the source B. The Rayleigh wave is poorly excited here, due to the depth of the source, but nevertheless is still visible on the snapshots.

Figure: The wave field propagation in a crossection of the homogeneous sphere model for an explosive source. The snapshots represent the -component of the displacement field inside , the spherical shell domain of the spectral element method, at different time steps of the propagation. The displacements in , the inner sphere of the modal summation method, are not shown here. The accuracy of the coupling method can be assessed here looking at the transmitted wavefield.
(click on image to zoom it)

Traces, of the vertical component of the displacement, are also recorded on the surface as a function of the epicentral distance, Fig. 9. Despite the depth of the source, the Rayleigh wave is clearly visible. As expected, the Rayleigh wave is almost non dispersive for this frequency source range, a small dispersion can only be observed,as expected, for the lowest frequencies. The X phase, which corresponds in terms of body waves to a superposition of P-wave subsurface reflections, is strongly excited. This phase is also, as expected, quite dispersive. Finally the P-wave is clearly observed as well as some multiples PP, PPP, ... It is worth to note here that this multiples tend slowly towards the X phase.

Figure: The vertical component of seismograms, recorded at the surface, as a function of the epicentral distance from the source B. Different phases are pointed out here: Rayleigh (R), X phase (X) and P wave (P).
(click on image to zoom it)

No spurious phase can be observed on Fig. 8 and 9. However, in order to assess the accuracy of the simulation, a direct comparison between the results obtained with the coupled method, for the three source positions, and the reference solution obtained by the normal modes summation is shown on Fig. 10. In each case the residual, i.e. the difference between the two solutions, has been multiplied by a factor ten. The agreement between the two solutions is indeed very good, with less than of error. Surface waves, which are strongly excited in the case of the the source A, are very accurately modeled by the spectral element method, as well as the body waves recorded directly on the coupling interface , which is an especially difficult test for this method.

Figure: The vertical component of the displacement, at an epicentral distance of , recorded on the surface (left) and on the coupling interface (right). The reference normal mode solution is drawn with the solid line, the coupled method solution is displayed with the dotted line and the residual between the two methods, amplified by a factor 10, is displayed with the bold dotted line. The method is very accurate and the maximum relative error is less than a few percents.
(click on image to zoom it)

Wave propagation in PREM

We consider now the spherically symmetric reference earth model PREM. The coupling interface is set at the CMB and the spectral element outer shell includes the whole mantle and the crust. This is in fact a quite challenging problem for the spectral element method due to the crustal structure inherent to the PREM model. The crust is characterized by slow wave velocities and the presence of very thin piecewise homogeneous layers with sharp elastic property contrasts: 15 and 10 km for the two uppermost layers. In the vertical direction, each of these concentric layers has to match with an element boundary in order to correctly approximate those contrasts. This drastically reduce the minimum size of the elements, in the vertical direction, and therefore put stringent constraints on the time step which must fulfill the CFL condition. Since each of these layers are piecewise homogeneous, the problem can be partly improved when using a very low vertical polynomial degree, a degree 2, within these crustal layers while retaining a degree 8 in the lateral directions (see Fig. 11).

Figure: The non conforming mesh used for the PREM example: below the km interface, the mesh has been derefined by a factor two in the horizontal direction. The color indicates the density structure of the model. The regions 1 and 4 have been removed here for visibility.
(click on image to zoom it)

Nevertheless, the time step imposed by the vertical structure of the crust in PREM is still of the order of s. Moreover the slow wave velocities in the crust imply laterally a good spatial resolution in order to accurately represent the surface waves.

The mesh used in this example has been built in order to match the nine surfaces of discontinuity in the mantle that are included in PREM. The mesh is non conforming in the vertical direction. In the upper mantle, the number of elements in the lateral direction is 32 elements while in the lower crust it is only of 16 elements allowing a rather uniform CFL condition all through the mantle. The non conforming interface is set at the km transition zone. This mesh could theoretically support a corner frequency of Hz, but in the results shown below, the actual corner frequency is Hz. The source is still explosive, for sake of simplicity, and its location is at Km.

The vertical component of the displacement recorded at the surface, for PREM, is shown on Fig. 12 as a function of time and of the epicentral distance. In order to identify some of the body wave phases, time arrivals of various phases have been computed by an asymptotic ray tracing, as shown on Fig. 13 and superimposed on the synthetic seismograms of Fig. 12.

Figure: The vertical component of the seismograms computed in PREM and recorded at the surface, as a function of the epicentral distance. Time arrivals of some body waves computed by a ray tracing has been superposed and can be identified using Fig. 13.
(click on image to zoom it)

Figure: Time arrivals of some body waves in PREM computed with a ray tracing.
(click on image to zoom it)

The comparison between the time arrivals, derived from a high frequency approximation and a full wave form modeling, derived from a direct numerical simulation, is not that obvious due to finite frequency effects on the whole wave form, especially for body waves. Actually the notion of the ``arrival time'' only makes sense within the high frequency approximation. As a matter of fact, as for the P-wave, the ray tracing arrival time points sometimes to the maximum of the phase amplitude, in which case the agreement appears quite good, while sometimes it points to the minimum of the phase amplitude, or to change of sign ... Obviously full waveforms obtained by direct numerical simulation does contain much more informations than traditional arrival time analysis based on asymptotic theories. However for the subject of this paper, such a comparison does not allow to assess the accuracy of the coupled method. A more thorough discussion of these results in the case of PREM and other laterally heterogeneous earth models will part of a forth coming paper.

The accuracy can be assessed by a direct comparison between the synthetic seismograms obtained by the proposed direct numerical simulation and the classical normal modes summation method. Such a comparison is shown on Fig. 14 and show a very good agreement.

Figure: Comparision between results obtained by the coupled method and modal solution method in PREM
(click on image to zoom it)

Some 3D runs

I briefly present here one 3D example done during my Ph.D.. It is perform in SAW12D Tomographic model (Li & Romanowicz 1996). A second one in a spherically symmetric model with a plume inclusion is available (ref).

SAW12D

Fig. 15 present the mesh used to perform this test in SAW12D. The S velocity contrast to PREM has been superimpose in color to the mesh. Deplacement recorded on four location are shown Fig. 16. To know more about this, you will have to read my Ph.D. Thesis (in French).

Figure: Mesh used for the test in SAW12D. In color is shown the S velocity contrast face to PREM.
(click on image to zoom it)

Figure: Radial componant of the deplacement recorded on four locations. Line 1 and line 2 are defined on Fig. 17
(click on image to zoom it)

Figure: Definition of line 1 (equator) and 2 (North-South) used to locate receivers and the crossection of the model beneath them.
(click on image to zoom it)

Plume

To be written ....: see Capdeville & all (2002)

Bibliography

1

Y. Capdeville.
Méthode couplée éléments spectraux - solution modale pour la propagation d'ondes dans la Terre à l'échelle globale.
PhD thesis, Université Paris 7, 2000.

2

Y. Capdeville, E. Chaljub, J. P. Vilotte, and J. P. Montagner.
Coupling the spectral element method with a modal solution for elastic wave propgation in global earth models.
Geophys. J. Int., 152:34-66, 2003.

3

Y. Capdeville, C. Larmat, J. P. Vilotte, and J. P. Montagner.
Direct numerical simulation of the scattering induced by a localized plume like using a coupled spectral element and modal solution.
Geophys. Res. Lett., 29, No. 9:10.1029/2001GL013747, 2002.

4

E. Chaljub.
Modèlisation numérique de la propagation d'ondes sismiques à l'échelle du globe.
Thèse de doctorat de l'Université Paris 7, 2000.
Available on http://www.ipgp.jussieu.fr/dept-sce/dmpn/Chaljub/index.html.

5

E. Chaljub, Y. Capdeville, and J. Vilotte.
Solving elastodynamics in a solid heterogeneous 3-Sphere: a spectral element approximation on geometrically non-conforming grids.
J. Comp. Physics, 183:457-491, 2003.

6

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

7

X. B. Li and B. Romanowicz.
Global mantle shear velocity model developed using nonlinear asymptotic coupling theory.
J. Geophys. Res., 101:11,245-22,271, 1996.