RegSEM: a versatile code based on the spectral element method to compute seismic wave propagation at the regional scale Free

Summary

1 Introduction

Solving the wave equation in realistic geological media is a crucial issue to properly model and study the propagation of seismic waves. At large scale, the effect of both crust and upper mantle 3-D structures on seismograms has been known for a long time (e.g. Montagner & Tanimoto 1991). At local scale, geological site conditions are now recognized as one of the dominant factors controlling the variations in ground motion (e.g. Olsen 2000). The accurate incorporation of geological structures in wave propagation modelling would therefore greatly improve the knowledge in fields such as tomography and site effects estimation.

Numerous techniques, like ray tracing or normal mode summation, have been developed for decades to compute synthetic seismograms. Ray tracing assumes that the seismic wavefield can be modelled as a large number of very narrow beams (e.g. Cervený 2001). For many purposes, this technique is very convenient, but it relies on a high frequency approximation, which means that it is only valid when seismic wavelength is much smaller than the scale of heterogeneity. In the case of low-frequency waves with large Fresnel zones, the ray theory no longer holds. The normal mode summation technique (e.g. Gilbert 1971), in addition with a high-order perturbation theory (Lognonné & Romanowicz 1990; Lognonné 1991; Clévédé & Lognonné 1996), is able to model waves with large Fresnel zones in 3-D Earth models, but the computation cost of such a technique is quickly prohibitive as the number of modes to couple increases with frequency. Moreover, normal mode perturbation methods are limited to weak lateral heterogeneities.

Direct numerical solutions have also been investigated to solve the wave equation. Finite differences have been applied (e.g. Alterman & Kara 1968; Boore 1972; Kelly et al. 1976; Virieux 1984, 1986; Olsen & Archuleta 1996; Moczo et al. 2007) but they present intrinsic problems in dealing with strong and deformed interfaces like basin edges, the Moho, the free surface and solid-fluid discontinuities. Such limitations do not exist in finite element methods, but the low polynomial order classically used in this kind of approaches make them inaccurate and dispersive when applied to elastodynamic problems (Lysmer & Drake 1972; Dupond 1973; Backer 1976; Marfurt 1984; Toshinawa & Ohmachi 1992). Since the 1990s, efforts have been focused on developing higher-order numerical modelling of seismic wave propagation. An important result from these efforts is the discontinuous Galerkin method (e.g. Dumbser & Käser 2006). Another major result, which is used all along this paper, is the spectral element method (SEM). Initially introduced in fluid mechanics (Patera 1984; Maday & Patera 1989), this method has been successfully applied to elastodynamics with the increasing concern of developing numerical techniques ensuring both a great precision and a reasonable numerical cost (Seriani & Priolo 1994; Faccioli et al. 1997; Komatitsch & Vilotte 1998; Seriani 1998; Komatitsch & Tromp 1999). The SEM was first applied at the global scale (Chaljub 2000; Komatitsch & Tromp 2002a,b; Komatitsch et al. 2002; Capdeville et al. 2003; Chaljub et al. 2003; Chaljub & Valette 2004). In the more recent years, applications to local (Komatitsch et al. 2004; Delavaud et al. 2006; Lee et al. 2008; Stupazzini et al. 2009; Chaljub et al. 2010; Peter et al. 2011) and continental scales (Chen et al. 2007a; Fichtner et al. 2009a) appeared. Such applications of the SEM proved that a great precision and a weak numerical dispersion can be obtained.

The SEM is a major contribution to seismology because it allows to compute the whole seismic wavefield propagating in a 3-D Earth model with no approximation on the wave equation (except minor numerical approximations). In the recent years, it was used to solve the inverse problem and get images of the Earth interior: Tape et al. (2009) and Tape et al. (2010) developed a model in southern California by inverting traveltime measurements down to 2 s period; Fichtner et al. (2009b) and Fichtner et al. (2010) obtained a model of the Australian region using a full waveform inversion down to 30 s period. To avoid the large number of simulations classically required to compute all the Fréchet derivatives, these authors implemented the adjoint technique (Tarantola 1984; Tromp et al. 2005; Fichtner et al. 2006a,b). This makes the computation of the gradient of the misfit function independent of the number of stations and parameters of the model. When dealing with a large set of sources, an alternative to the adjoint technique is the scattering-integral approach (Chen et al. 2007b). In any case, with the current computational power, one cannot solve the inverse problem using the SEM with a classic procedure.

In this paper, we focus our attention on the forward modelling. When using the SEM to compute 3-D seismic wave propagation, the main practical issue consists in meshing the medium using hexahedra. Indeed, the mesh has to honour the discontinuities of the geological model under study to fully benefit from the accuracy of the method and properly model effects associated with these discontinuities such as wave diffraction. Because realistic models of the Earth often have complex geometry, lots of efforts and time are usually needed to build an appropriate mesh. Casarotti et al. (2008b) developed automatic procedures to create 3-D unstructured hexahedral meshes, but it is still not possible to generate in a fully automatic way meshes that would honour detailed geological discontinuities such as realistic sediment-rock interfaces. The goal of this paper is to promote a code, RegSEM, that can accurately take into account 3-D discontinuities in regional meshes and then compute seismic waves within them using the SEM. More precisely, RegSEM has two versions: a continental version which is able to generate structured meshes of crustal and mantle structures separated by a 3-D Moho in spherical geometry, and a local version which uses an external mesh generator, CUBIT (), to produce 3-D unstructured meshes essentially designed to study the seismic response of sedimentary basins.

In a first part, RegSEM′s features are described, including the different kinds of meshes the code can provide and/or handle. In a second part, waveforms computed with RegSEM in simple elastic models are compared to analytical or semi-analytical solutions. These comparisons allow to validate the code. The last part presents simulations of a real earthquake in different tomographic models of Europe. A comparison between synthetic seismograms and real data enables to investigate the improvement of the fit when using 3-D models and point out the great influence of the shallow layers. The large number of simulations shown in this work serves as a set of examples to put in evidence the capability and the great versatility of RegSEM. The code sources, in addition with a manual and several examples, can be downloaded at .

2 Regsem′s features

2.1 The spectral element method

The SEM was developed in fluid dynamics in the 1980s (Patera 1984; Maday & Patera 1989) and was adapted to elastodynamics in the 1990s (Seriani & Priolo 1994; Faccioli et al. 1997; Komatitsch & Vilotte 1998; Seriani 1998; Komatitsch & Tromp 1999). The SEM is similar to a finite element method. It is based on a primal variational formulation of the equations of motion. This formulation enables to naturally take into account both interface and free boundary surface conditions, allowing a good resolution of evanescent interface and surface waves.

The discretization process implies the decomposition of the spatial domain into non-overlapping elements. Classical implementations of the SEM in computational seismology are based on hexahedral elements to benefit from advantageous properties of tensorization. Although hexahedra are less favourable than tetrahedra for meshing geometrically complex structures, a certain flexibility is ensured by a local geometrical transformation from a reference element (unit cube) to any deformed element, as detailed in -1">Section 2.2.1. Unstructured meshes offer additional possibilities, as seen in -2">Section 2.2.2.

Associated with the domain decomposition, the functional discretization is based on a piecewise high-order polynomial approximation. The specificity of the SEM holds in the choice of basis functions intimately related to the Gauss-Lobatto-Legendre (GLL) quadrature used to evaluate the integrals in the variational formulation. The basis is obtained from the orthogonal Lagrange polynomials associated with (N+ 1) interpolation nodes (where N is the polynomial order). These nodes are chosen to be the nodes of the GLL quadrature. Such nodes define a tensor product grid where the displacement, its spatial derivatives and products encountered in the variational formulation are evaluated. The choice of a Lagrangian interpolation associated with the GLL nodes gives the SEM a very interesting convergence property: an increase of the polynomial order leads to an exponential diminution of the aliasing error. This property, called spectral precision, gives its name to the method.

Inserting the polynomial interpolation and quadrature rules into the variational form of the equations of elastodynamics leads to a system of ordinary differential equations governing the evolution at the global nodal position, which can be written as follows:

1

2

where U, V and are vectors containing the components of the displacement, velocity and traction at the global nodes, respectively. is the mass matrix. The vectors and contain the external and internal forces, respectively, and corresponds to the traction forces. The use of an orthogonal basis defined as the Lagrangian functions associated with the GLL nodes leads to a second interesting characteristic of the hexahedral version of the SEM: the mass matrix is diagonal. This property enables to use an explicit time stepping in which the inverse mass matrix can be exactly computed. In RegSEM, like in most of the SEM implementations, this time stepping is a second-order finite difference scheme.

2.2 Meshing a chunk of the Earth with hexahedra

When using the SEM for 3-D complex geological models, the main difficulty consists in meshing the model under study. Indeed, to benefit from the high accuracy of the method and properly model the effects linked with the geology, the mesh has to adapt the velocity structure of the model, in particular the zeroth-order discontinuities such as the Moho and the sediment-rock interface. This task often requires significant efforts and time.

In this section, we describe two kinds of meshes in which RegSEM can simulate seismic wave propagation. The first kind corresponds to regular meshes of crustal and mantle structures in spherical geometry. These meshes are suitable to compute seismic waves for source-receiver distances ranging from 1 to 90°. RegSEM not only handles such meshes but it can also create them in a versatile way: the size and the location of the chunk, the spherical discontinuities and possible 3-D Moho and 3-D free surface can be defined by the user. The second kind of meshes requires more efforts from the user because it has to be generated externally. It corresponds to unstructured meshes. Such meshes can deal with more complex geometries. In the following, they are used to study the seismic response of sedimentary basins, but they could be designed for other applications.

Before giving more details on the capabilities of RegSEM, it must be pointed out two important conditions that all meshes have to fulfill to make the SEM accurate and stable:

1. (i)
   In classical applications of the SEM, 4 ≤ N ≤ 8. For such values of N, at least five GLL nodes per wavelength are needed everywhere in the region to properly describe the seismic wavefield (e.g. Komatitsch & Vilotte 1998). This means that the size of the elements d and the polynomial order N are both constrained by the shortest wavelength λ_min propagated in the medium. This condition can be summarized by the following relation:

   

   3
2. (ii)
   To ensure the stability of the time-marching, the time step Δt of the finite difference scheme has to verify the Courant-Friedrichs-Lewy (CFL) condition:

   

   4

   where C denotes the Courant number, usually chosen between 0.3 and 0.4, and the minimum ratio of grid spacing Δx (distance between two GLL nodes) and P-wave speed a.

2.2.1 Regular meshes of crustal and mantle structures

RegSEM can provide a regular mesh of any chunk of the Earth whose lateral size is smaller than 90°. To do so, the code uses the so-called cubed sphere mapping (Sadourny 1972; Ronchi et al. 1996). For each element, this mapping allows to define the Cartesian coordinates of 27 control points. Using the Lagrange polynomials of degree 2 associated with these control points, the unit cube can be deformed and the shape and position of each element in the chunk can be defined. Such a classical procedure enables to easily design a structured and conformal mesh for any section of the Earth (Fig. 1). Moreover, ellipticity can be taken into account using the Clairaut′s equation (Dahlen & Tromp 1998).

Figure 1

A chunk of the Earth meshed by hexahedral elements. It is 20°× 40° large and 1400 km thick. The elements are gathered under different colours. Each colour represents a subdomain. Here, the chunk is divided in eight subdomains. The star and the triangle correspond to the source and the station used in Section 3.1.

Open in new tabDownload slide

As mentioned above, the seismic discontinuities in the velocity model have to be honoured by the mesh. If all the discontinuities are spherical, then RegSEM is quite versatile: one just needs to introduce the radius of each discontinuity, and then the code fills the seismic layers with the appropriate number of elements. Of course, this number depends on the vertical size of the elements. This size is first equal to the horizontal size d_h introduced by the user, and then it is adjusted in each seismic layer to fit the thickness of the layer. In the case of PREM (Dziewonski & Anderson 1981), there is one more level of sophistication because not only the thickness of the layers but also the seismic velocities in the layers are used by the code to constrain and optimize the vertical size of the elements. Fig. 2 shows an example of a mesh of PREM. Some elements within this mesh appear to have a large aspect ratio (up to 5). In the context of the SEM, this is not a problem: accuracy is preserved thanks to the high spatial degree of the method (Oliveira & Seriani 2011) and stability is kept up because the shear deformation of the elements is small. Examples of simulations in PREM are shown in Section 3.2.

Figure 2

A mesh for PREM. The shallowest fluid layer has been replaced by the underlying solid. The chunk is 20°× 40° large and 1400 km thick. The horizontal size of the elements is d_h= 0.44°. This allows the propagation of a 20 s period wavefield using a polynomial order N= 4. The S-wave speed at each GLL node on the vertical border of the domain is plotted. A zoom into the upper part of the model shows the discontinuities. The star and the triangle correspond to the source and the station used in Section 3.2.

Open in new tabDownload slide

RegSEM can also mesh any surface and Moho topography (Fig. 3). This is an important feature because the crust has significant effects on surface waves (Montagner & Tanimoto 1991; Curtis et al. 1998; Komatitsch et al. 2002; Shapiro & Ritzwoller 2002; Marone & Romanowicz 2007; Ferreira et al. 2010), even at relatively long period (up to about 60 s). The capability to consider any model with a realistic crust is therefore a major benefit. To do so, the code uses only one layer of elements in the crust. This means that discontinuities within the crust, such as the sediment-rock interface and the upper-lower crust interface, cannot be taken into account. Moreover, the fact that only one layer of elements is used to mesh the crust limits the frequency content that the simulation can handle. For example, when performing a simulation in Tibet (where the crust is more than 70 km thick) with a polynomial order N= 8, the highest frequency to be propagated will be approximately 0.1 Hz. Examples of simulations in 3-D crustal models are shown in Section 4.

Figure 3

A mesh for CUB in the Atlantic-European region. The chunk is 30°× 70° large and 1500 km thick. The horizontal size of the elements is d_h= 0.35°. This allows the propagation of a 20 s period wavefield using a polynomial order N= 4. The S-wave speed at each GLL node on the border of the domain is plotted. A zoom into the upper part of the southern side of the chunk shows the discontinuities. The transition between a thin oceanic crust and a thick continental crust is visible. A map of the Moho corresponding to the present chunk is shown in Fig. 12.

Open in new tabDownload slide

Realistic models of the Earth all have thin shallow layers. This is the case in Figs 2 and 3. Because of the CFL condition 4, these thin layers imply a very small time step, which makes the computation cost high:

1. (i)
   The mesh shown in Fig. 2 is designed to propagate a wavefield with a minimum wavelength λ_min= 60 km using a polynomial order N= 4. The Earth model is PREM. In this case, ratio is minimum in a thin layer defined by two discontinuities at 15 and 24.4 km depth. Indeed, the elements used to mesh this layer have a vertical size of 9.4 km, which is very small compared to the size of the other elements. The layer has a P-wave speed α = 6.8 km s^-1, so we find

   

   5

   where γ_N is a coefficient which depends on N and which comes from the fact that the GLL nodes are non-evenly spaced in an element. For N= 4, γ_N ≃ 6. Combining eqs 4 and 5 and taking C= 0.35, we obtain Δt ≥ 0.080637 s, which is small and makes the computation cost high. Note that in practice, RegSEM computes the time step automatically: it first finds by a grid search and then determines Δt using 4 with C= 0.35.
2. (ii)

   The mesh shown in Fig. 3 is designed to propagate a wavefield with a minimum wavelength λ_min= 50 km using a polynomial order N= 4. The region under study is the Atlantic ocean and Europe. The velocity model is CUB (Shapiro & Ritzwoller 2002), which has a realistic crust. Such a crust is thin below the ocean: the minimum thickness is about 7 km, so there are some elements that have a very small vertical size compared to the size of the other elements. Ratio is minimum in these elements. Assuming a P-wave velocity α = 5 km s^-1 in the oceanic crust, we find the following condition on the time step: Δt= 0.081667 s.

To avoid the small time step induced by a thin layer, recent works developed techniques to replace the layer by a thicker effective medium (Capdeville & Marigo 2008; Fichtner & Igel 2008; Lekic 2010). These techniques yield a new layer with a large and constant thickness at the top of the Earth. Although this kind of effective layers is easy to mesh, we do not show examples of computation in such media in this work.

2.2.2 Unstructured meshes of sedimentary basins

Basin effects are characterized by scattering, focusing and basin-edge induced surface waves which are closely associated with the geometry of the basin. These effects are recognized to be responsible for a long duration of the seismic signal in the basin and especially for large local amplifications. To accurately model and study these effects, especially at high frequencies, it is important that the sediment-rock interface is honoured. Considering the complexity of most of these discontinuities, 3-D unstructured meshes are necessary to achieve this goal. The local version of RegSEM has the ability to handle such unstructured meshes where topology is totally arbitrary. To do so, the code is written according to a strategy of independence against the Cartesian coordinates. This enables to handle the random orientation of the four different objects (elements, faces, edges and vertices) which compose the mesh. Defining such objects allows to assign specific actions to each of them, such as Neumann conditions (Section 2.5).

The creation of the mesh is not performed by RegSEM; it is done externally using the CUBIT mesh generation tool kit (). Input mesh files for RegSEM are then created from export CUBIT mesh files. Considering the limited choice of commercial and non-commercial codes dealing with hexahedra compared to the case of tetrahedra, we think that this tool kit offers the best alternative. However, the CUBIT mesh generation for 3-D complex structures is not totally automatic and requires many steps and user interventions. Automatic procedures have been developed by Casarotti et al. (2008b) to generate meshes according to different strategies. Embedded in a parallel Message Passing Interface (MPI) environment, they can fast create simple 3-D unstructured meshes. However, in the case of a complex basin, the outcrop is generally not honoured when this one exhibits too many variations. In this case, a robust fully 3-D unstructured algorithm for hexahedra is still not available.

The meshing technique we use in this paper has been successfully applied to the valley of Grenoble by Stupazzini et al. (2009). It consists in building a conform mesh of the model from separately meshed subdomains, to better control the size and shape of the elements. This work requires a pre-process with a CAO software which can handle NURBS (Non-Uniformal Rational B-Splines) geometries. From the digital terrain model data, this software creates NURBS curves which define the interfaces of the 3-D model, the topography, the basin basement and the numerical boundaries. The total volume is then partitioned into subdomains which are also defined by a group of curves and exported in an Initial Graphics Exchange Specification format readable by CUBIT. These subdomains are independently reconstructed by CUBIT which assembles them in a conform way to form the whole domain. Each subdomain can then be individually meshed. A more detailed description of this meshing process can be found in Delavaud (2007) and in the manual of the code. The advantage of this procedure holds in the possibility of associating to each subdomain different element sizes and types of meshing while the mesh remains conform. The number of subdomains can be substantial, depending on the complexity of the structures which have to be meshed. In the presence of a basin, this partitioning is mainly controlled by its outline in surface which needs to be isolated to correctly mesh the bend and the shape of the edge. Fig. 4 shows one of the subdomains which defines the inner outline belt of the Caracas basin. One of the triangular surface at the ends of the subdomain is first meshed with a meshing scheme called ‘triprimitive’ which applies to three side surfaces. The size of the elements is inherited from an interval size assigned to the edges. The total volume is then meshed by a projection (sweep) of the two dimension mesh along the edges towards the opposite triangular surface. As one can see, the element at the edge is particularly deformed and introduces a very small minimum distance between the GLL nodes. The variation of the free surface topography also influences the cutting into subdomains needed to ensure a homogeneous mesh in depth. The meshing strategy remains the same as the one described for the outline of the basin: the free surface is first meshed with an unstructured meshing scheme, then a sweep is applied in the vertical direction. Therefore, the mesh is unstructured only in two directions (in depth for the interior outline belt of the basin, and horizontally for the other parts). As an example, the case of a simple hemispherical basin is presented in Fig. 5.

Figure 4

Detail of a subdomain as part of the inner outline of the Caracas basin. The volume was meshed from the projection of the 2-D triangular front mesh.

Open in new tabDownload slide

Figure 5

Mesh of a quarter of a half space containing a hemispherical basin. The model is cut into four quarters independently meshed using CUBIT and then reassembled. The basin (yellow) is meshed first, then the free surface, and finally the volume. The blue elements correspond to the PML.

Open in new tabDownload slide

Meshes which honour geological discontinuities might present highly deformed elements, especially at the edges. To assess the quality of a mesh and identify such elements, CUBIT offers different metrics. Skewness, distortion or shear associated with the Jacobian of each element are possible quality measurements. In the case of the SEM, the effect of deformed elements on precision and efficiency is fortunately limited by the high degree of the method (Oliveira & Seriani 2011).

A recent study by Pelties et al. (2010) compares the method that interpolates the outline of a basin by Lagrange polynomials (e.g. Komatitsch et al. 2004; Casarotti et al. 2008b; Lee et al. 2008) with a meshing that fully honours it (e.g. Stupazzini et al. 2009). From tests performed at different frequencies in different velocity contrasts, this study provides empirical rules to ensure the reliability of Lagrange interpolation. In the case of 3-D simulations in the Grenoble valley, Casarotti et al. (2008a) observed differences in amplitude and phase of the order of 15 per cent between the two methods above 1 Hz. At lower frequencies, the detail of the basin shape has less, or even no, influence. Although it is not the scope of this paper to present a comprehensive study about the differences between honouring and not honouring interfaces, we briefly show some comparisons based on a 2-D profile of the Caracas basin. The effects of three meshing strategies on spectral ratios are presented. Spectral ratios are amplification factors with respect to the spectrum of the incident plane wave. We compute them for frequencies up to 5 Hz and for each receiver along the free surface. In Fig. 6, on the left panel, the ratios are shown for a mesh that respects the corners of the basin. These corners are critical because they generate diffracted surface waves responsible for large amplifications. Bi-dimensional effects are characterized by amplifications at frequencies higher than the fundamental frequency (0.5 Hz), especially above the thicker part of the basin at 1.5 Hz, 2.5 Hz and 3.5 Hz. We also consider a regular mesh in which the discontinuity is interpolated at the GLL nodes (Fig. 6, mid panel). In this case, we show the ratio in logarithmic scale of the spectra recorded at the free surface of the regular mesh and those recorded at the free surface of the honouring mesh. Up to 2 Hz, the spectra are almost similar. At higher frequencies, large discrepancies appear, with ratio of the order of 2. A third meshing strategy consists in respecting the interface until a depth of 60 m, about 50 per cent of the minimum propagated wavelength, and to interpolate the discontinuity in the elements at the edges of the basin (Fig. 6, right panel). The spectra obtained in this case are very similar to the fully honouring case up to 2.5 Hz. Above this frequency, large but localized discrepancies appear, also of the order of 2. This example put in evidence that the entire rock-sediment interface, especially the edge of the basin, should be respected for high frequency simulations (>2.5 Hz). For shorter frequencies, the direct discretization of the interface at the GLL nodes down to a reasonable depth seems appropriate. Advantages and drawbacks of these strategies in terms of precision and implementation complexity could be discussed in more detail based on the analysis conducted by Maday & Rønquist (1990).

Figure 6

Effects of three different meshing strategies for the Caracas basin. Left panel: the sediment-rock interface is fully honoured. The spectral ratios computed at the free surface are shown. Mid panel: the interface is interpolated by a thin GLL grid. The ratio (in logarithmic scale) of the spectra recorded at the surface and the spectra computed in the honouring mesh is plotted. Beyond 2 Hz, significant differences appear. Right panel: the whole interface is honoured except its edges. In this case, comparing the spectra with those computed in the fully honouring mesh shows large but localized discrepancies at high frequency.

Open in new tabDownload slide

2.3 Introduction of the elastic parameters

As RegSEM is able to handle a large set of meshes, it has to be versatile in introducing elastic models as well. In the continental version of the code, this is achieved thanks to a Fortran module that the user can change by himself and which is conceived to provide the elastic parameters at any location (radius, latitude and longitude) in the Earth. Both radial and azimuthal anisotropies are implemented. Moreover, the anelastic structure can be taken into account using a series of standard linear solids, as suggested by Emmerich & Korn (1987) and Komatitsch & Tromp (1999). The unstructured version of the code is limited to lossless isotropic media.

Rotation and self-gravitation, which involve non-neglectible effects at the global scale and very long periods (Komatitsch & Tromp 2002b; Chaljub & Valette 2004), are not included. Propagation in fluid has not been implemented, so the waves either from the outer core or from the oceans cannot be simulated. Nevertheless, following Komatitsch & Tromp (2002b), the mass of the oceans can be taken into account when a bathymetry is used at the surface of the Earth.

2.4 Absorbing boundary conditions

To avoid artificial reflections at the border of the chunk, it is necessary to implement efficient absorbing boundary conditions. RegSEM uses the velocity-stress formulation of the so-called Perfectly Matched Layers (hereafter PML; see Festa & Vilotte 2005). This formulation requires an unphysical splitting of the field variables along the directions of normal and parallel derivatives with respect to the interface PML volume. This means that in practice, the splitting directions have to be known at every GLL node belonging to the PML, which is not obvious when working with Cartesian coordinates in a deformed layer (such as the lowermost layer of the chunks presented in Figs 1-3). Therefore, we here make an assumption: for all the GLL nodes of a given element, the splitting directions defined at the centre of the element are used. The effect of such an assumption on the stability and accuracy of the PML is not clear. Furthermore, our PML are isotropic, so spurious reflections can be created when considering anisotropic media. Examples in the following parts of the paper will show that our implementation of the PML however provides satisfactory results.

Note that the SEM does not require to use the same polynomial order for all directions. Our code takes advantage of this flexibility: in the PML, it is possible to use a different polynomial order in the damping directions. In the following, we will always use N= 8 in these directions.

2.5 Sources

Force-vector and moment-tensor point-sources can be placed at any location in the chunk. Four different functions are provided to describe the time signal at the sources: a Gaussian, its derivative (i.e. a Ricker wavelet), its antiderivative (i.e. an error function) or the Fourier transform of a frequency band defined by a cosine taper.

More interestingly, in particular for basin response modelling, an incident wavefield can be introduced in the unstructured version of RegSEM. The method developed for that purpose is based on a decomposition technique and exploits the natural presence of the traction in the SEM formulation, that is the term in eq. 1. The wavefield is introduced on an interface in the domain, for example, the sediment-rock interface, by its action on the traction forces (Neumann condition). Similar ideas had been introduced by Bielak & Christiano (1984) in the context of finite elements for the problem of soil-structure interaction. The main interest of this method consists in avoiding the propagation of the incident field, which is known analytically or numerically, as long as it has not reached any discontinuity with which it will interact by reflection, transmission or diffraction. This type of introduction is compatible with any boundary conditions, including PML. Moreover, diffraction problems for non-vertical incidences are prevented. Finally, the computational domain does not have to be large to hold the incident wave. We refer to the manual of the code for more details about the implementation of this method.

2.6 Parallel implementation

The SEM can be easily implemented on distributed memory architectures. Given a number n of CPUs, the computational domain has to be divided into n subdomains. To do so, we use the software library METIS () that ensures an efficient partitioning which minimizes the communications between the subdomains. These communications occur at every time step of the time-marching scheme. To perform them, we use the MPI. Fig. 1 shows an example of a chunk partitioned by METIS. Fig. 7 shows the good scalability of our parallel implementation in the case of the experiment described in Section 3.1.

Figure 7

Scalability of the RegSEM parallel implementation. The simulation used to do this test is the one described in Section 3.1. The points are almost aligned, showing that the computation time goes like the inverse of the number of CPUs.

Open in new tabDownload slide

3 Validation

In this part, a series of numerical experiments are carried out to validate our code. We start with the simple case of a homogeneous medium, then we consider a layered medium (PREM) and we finally study the case of a hemispherical basin. For all these experiments, a reference solution is known.

3.1 Simulation in a homogeneous medium

We first consider a homogeneous medium in spherical geometry. In such a context, the normal mode summation technique provides a quasi-analytical solution (Capdeville 2000). The P-wave speed, S-wave speed and density of the medium are α = 8 km s^-1, β = 5 km s^-1 and ρ= 3000 kg m^-3, respectively. The chunk used in the SEM simulation is shown in Fig. 1. The elements are 1.3° large. This enables to use a 50 mHz cut-off frequency with a polynomial order N= 8. The source is an explosion located at 10 km depth. The receiver is on the free surface. The epicentral distance is 20°. On three Intel Xeon 2.5 GHz quad-core dual-processor nodes (i.e. 24 CPUs), it takes about 16 min to compute 1200 s. In Fig. 7, the computation time for other numbers of CPUs is shown. The scalability of our parallel implementation is seen to be good.

Fig. 8 shows the comparison between the SEM result and the normal mode solution. All the seismograms are normalized with respect to the amplitude of the vertical component obtained with the SEM, so the relative amplitudes are preserved. Because the source is an explosion, there is no SH energy and the waves only lie on the vertical and radial components. On these components, the traces obtained from the two methods are indistinguishable. The plot of the residual multiplied by 10 shows that the maximum error is around 2 per cent. This small error is essentially due to the finite difference time-scheme whose order is only 2. On the transverse component (whose y-axis is 100 times larger than the other components), the SEM solution is not exactly zero: around 650 s, the residual multiplied by 10 reveals a spurious reflection coming from the bottom of the chunk. This signal is extremely small and is hardly seen on the radial and vertical components, meaning that our PML are good.

Figure 8

Comparison of the SEM solution (dashed black) with the normal mode solution (red) obtained in a homogeneous medium (a= 8 km s^-1, β = 5 km s^-1, ρ = 3000 kg m^-3). The epicentral distance is 20°. The residual multiplied by 10 is plotted in green.

Open in new tabDownload slide

3.2 Simulation in PREM

We use PREM to perform a second validation test. The thin fluid layer which lies on the top of this model is replaced by the underlying solid. The chunk used in the SEM simulation is shown in Fig. 2. The elements are 0.44° large. This enables to use a 50 mHz cut-off frequency with a polynomial order N= 4. The source-receiver configuration is the same as in the previous test.

3.2.1 With no attenuation

We first do not take into account the anelastic structure of PREM. In this case, the normal mode summation gives a good reference solution. Fig. 9 shows a comparison between this solution and our SEM result. Again, the two waveforms are indistinguishable on both vertical and radial components. Nevertheless, on the radial component, the residual shows significant amplitudes after the main phases, which was not the case in the test performed in the homogeneous medium. This is explained by two reasons. First, our chunk is cut at depth, so the phases reflected at the core-mantle boundary, such as the PcS and ScS phases, are missing. Second, a spurious reflection from the PML on the vertical sides of the chunk is detected. This signal appears here because PREM is anisotropic. The magnitude of this spurious reflection is similar to the error due to the finite difference time-scheme, so eventhough they are not perfect, our PML are satisfactory. On the transverse component, a tiny reflection coming from the bottom of the chunk is observed, as it was the case in the homogeneous medium. This reflection arrives earlier here because PREM velocities increase with depth.

Figure 9

Comparison of the SEM solution (dashed black) with the normal mode solution (red) obtained in PREM with no attenuation. The epicentral distance is 20°. The residual multiplied by 10 is plotted in green.

Open in new tabDownload slide

3.2.2 With attenuation

When introducing the anelastic structure, the comparison between the SEM and the normal mode solutions (Fig. 10) does not change a lot. The main difference with the lossless case is a large residual value (5-10 per cent) for the body waves and the beginning part of surface wavetrain. This large value is probably due to the fact that the normal mode result in an attenuating medium relies on a first-order approximation. Nevertheless, we can conclude that the attenuation is well taken into account in RegSEM because the SEM and normal mode summation waveforms are very similar to each other and are both very different from the lossless case.

Figure 10

Comparison of the SEM solution (dashed black) with the normal mode solution (red) obtained in PREM with attenuation. The epicentral distance is 20°. The residual multiplied by 10 is plotted in green.

Open in new tabDownload slide

3.3 Plane wave on a hemispherical basin

To validate our code with an unstructured mesh, we now consider the half space containing the hemispherical basin of radius R (Fig. 5). For such a simple shape, the boundary element method developed by Sánchez-Sesma (1983) provides a semi-analytical solution of the diffracted wavefield in the frequency domain. The parameters of the material are summarized in Table 1. They respect the conditions on the shear modulus µ, density ρ and Poisson coefficient υ described by Sánchez-Sesma (1983): µ^R/µ^E= 0.3, ρ^R/ρ^E= 0.6, υ^R= 0.3 and υ^E= 0.25, where the exponents R and E correspond to the basin and the rest of the half-space, respectively. The model is excited by a plane P-wave with a vertical incidence. The results are presented for a normalized frequency , where Γ_P is the wavelength of the incident plane P wave. This normalized frequency determines the central frequency f_P of the Ricker pulse which defines the time function of the incident plane wave. For the normalized frequency η_P= 0.5 considered by Sánchez-Sesma (1983), f_P= 5.76 Hz.

Table 1

Open in new tabDownload slide

In the SEM simulation, a polynomial order N= 4 is chosen in the elastic medium. The signal is recorded at the free surface, along the line [Ox), where O is the centre of the basin and x a horizontal axis. We are interested in the transfer function, which is the ratio between the spectrum of the signal recorded at each receiver along this line and the spectrum of the incident plane wave. Fig. 11 shows the transfer functions, for the SEM and the semi-analytical solution, of the vertical and horizontal components of the displacement at the normalized frequency η_P= 0.5 as a function of the normalized distance x/R at t= 5 s. These transfer functions are also represented at the normalized frequency η_P= 0.7 which corresponds to the frequency f_P= 8 Hz. As a result, the two methods provides the same behaviour. Small differences of 2 per cent can be observed, which is comparable to the numerical dispersion of the SEM and the magnitude of the PML reflections.

Figure 11

Transfer functions of the vertical and horizontal components of the displacement as a function of the normalized distance x/R at t= 5 s. We compare the SEM (stars) and a semi-analytical solution (circles) based on an expansion of Bessel functions (Sánchez-Sesma 1983) in the case of a hemispherical basin, for two frequencies: 5.76 Hz (η_P= 0.5) on the left and 8 Hz (η_P= 0.7) on the right.

Open in new tabDownload slide

A variability of the behaviour according to the frequency is observed, knowing that the 1-D resonance frequency of the basin equals a^R/4R= 4.4 Hz. At 5.76 Hz, the curve is close to the 1-D case, with a maximum amplification at the centre of the basin, about 175 per cent of the amplitude obtained in the case of a homogeneous medium. The amplitude decreases when approaching the basin border and then keeps the state of a homogeneous medium. The horizontal component, on the contrary, exhibits a null amplitude at the centre of the basin, due to the symmetry, and reaches its maximum at x/R= 0.5. Then it converges to a stable amplitude. At 8 Hz, an important wave conversion appears around x/R= 0.35 where the horizontal component reaches its maximum and the vertical component its minimum. Outside the sphere, the two components tend to the same stable state as at 5.76 Hz. The maximum amplification at the centre of the basin is reached for a normalized frequency η_P= 0.61 (f_P= 7 Hz), 1.6 times the 1-D resonance frequency. These results underline the influence of the basin effects, even in a simple and symmetric case.

4 Simulations in 3-D models of Europe

In this part, we point out the influence of the shallow structures (crust and uppermost mantle) on a regional wavefield. RegSEM′s flexibility is used to generate structured meshes and simulate an event in two 3-D models of the European region. The models are Crust2.0 (Bassin et al. 2000) over the 1-D model PREM, and the crustal and upper mantle model CUB (Shapiro & Ritzwoller 2002). The corresponding synthetic seismograms are compared with the normal mode solution in PREM, and above all with the real data seismograms.

4.1 Models

Crust2.0 (Bassin et al. 2000) was obtained by compiling seismic studies and tectonic, geological settings. The model is defined as cells on a 2 × 2° geographical grid. In each cell, we average over the different layers (there are up to seven layers, including ice or water, sedimentary layers and layers in the bedrock) to get homogeneous seismic velocities and density. Then, a 2-D Gaussian filter is applied horizontally to set up the seismic parameters at every GLL nodes. Below the 3-D Moho, PREM is used to represent the deeper Earth.

CUB (Shapiro & Ritzwoller 2002) is a radially anisotropic model obtained by surface wave tomography. The model is derived from the inversion of group and phase dispersion data of the fundamental mode of both Rayleigh and Love waves in the 16-200 s period range. It is defined on a 2 × 2° geographical grid. It is characterized by a vertical block parametrization for the crust, a 3-D Moho and a spline parametrization for the upper mantle. We actually use a smoothed version of the CUB crust: the fluid and low-velocity ( β = 2.4 km s^-1) shallow layers are replaced by the underlying material, and a vertical smoothing is applied within the crust using the intrinsic interpolation law of the SEM. Moreover, because the resolution of CUB is poor deeper than 250 km, a linear transition towards PREM is imposed down to the transition zone. Again, we use a Gaussian filter in the horizontal directions. The resulting model in the Atlantic-European region is shown in Fig. 3. A map of the Moho is shown in Fig. 12.

Figure 12

Map of the Moho of model CUB (Shapiro & Ritzwoller 2002) in the Atlantic-European region. The source-receiver configuration used in Section 4.3 is also shown.

Open in new tabDownload slide

4.2 Data

The earthquake we consider has been recorded at 16 receivers on the continental Europe. All data come from broad-band stations operated by the GSN, GEOSCOPE, GEOFON and MEDNET networks. Both vertical and horizontal components on the BH channels are selected. To be compared to the data, the synthetic seismograms are convolved with the instrumental transfer function. This avoids the deconvolution of the response on the data.

The comparison between data and synthetics is performed in four period bands: 100-200 s, 50-100 s, 30-50 s and 20-30 s. To obtain a quantitative estimate of the misfit between data and synthetics, we apply a systematic cross-correlation between the two using a group velocity criterion. Two parameters of the cross-correlograms are extracted: the delay dt of the main peak and the amplitude ratio R_A between the main peak and the autocorrelation of the data. A perfect matching therefore yields dt= 0 s and R_A= 1. Our comparison is mainly based on the vertical component. The signal-to-noise ratio in the radial and transverse components is often poor, so a more severe data selection would actually be necessary to include more horizontal components in our study. Note that the seismograms in the 100-200 s period band are not presented because they show the same behaviour as those in the 50-100 s period band. Furthermore, in the shortest period band (20-30 s), waveforms are complex and some correlation parameters therefore do not make sense.

4.3 Results

The event we investigate occurred along the Mid-Atlantic Ridge on 2010 May 25, at latitude 35.41°N and longitude 35.93°W. Its depth has been estimated to 10 km and its USGS CMT solution provides a moment magnitude of 6.3. The regional chunk used for the SEM simulations of this event is shown in Fig. 3. A map view of the chunk can be found in Fig. 12, representing also the source-receiver configuration. The seismograms obtained at four relevant stations are shown in Fig. 13. The correlation parameters for these stations and the mean of the parameters over all the stations are presented in Table 2.

Figure 13

Waveforms induced by the mid-Atlantic ridge earthquake at four stations in Europe for three different period bands. We compare the real data (grey) with the waveforms obtained in three different Earth models (PREM in blue, PREM+Crust2.0 in red and CUB in green).

Open in new tabDownload slide

Table 2

Open in new tabDownload slide

In the long-period bands (100-200 s, 50-100 s), the synthetic waveforms are similar to the data with a correct amplitude. However, at every station, a systematic positive delay is observed between the data and the synthetics computed in PREM or in Crust2.0. This delay increases with the epicentral distance. CUB shows a much better fit: for example, at station KIEV, models PREM and Crust 2.0 present a large delay (dt= 20 s and dt= 34 s, respectively) whereas CUB shows a short delay (dt=-6 s). In the medium band (30-50 s), the waveforms are still similar and the amplitude is still correct. The CUB model improves the fitting to the data in comparison with what PREM and Crust2.0 do. However, when increasing the epicentral distance, the trends of the last two models differ. On stations SSB and ECH (short paths), PREM and Crust2.0 both show the same large positive delay. On stations AQU and KIEV (longer paths), the Crust2.0 model presents an increasing delay (dt= 95 s) while PREM model presents a decreasing delay (dt= 7.8 s). In the short-period band (20-30 s), the waveforms significantly differ, in particular when increasing the epicentral distance. The Crust2.0 model always shows the worst waveforms and sometimes have nonsense delays (dt= 165 s at KIEV) due to the strong coda wavetrain. At short epicentral distance (SSB for example), the CUB model presents very good fits (dt= 0.3 s) and is better than PREM (dt= 63 s). At longer epicentral distance (KIEV for example), situation is reversed: PREM is better (dt=-3.8 s vs. dt=-28 s).

In the short-period range, the effect of the crust on delay times is predominant. For oceanic paths, PREM has a too thick crust (slow velocity compared to fast velocity in the underlying lithosphere) which tends to induce positive delay times. The situation is reversed in continents, which are characterized by a thick crust and fast lithospheric velocities. Consequently, when only the effect of the crust is taken into account in the numerical simulations, the delay time tends to increase compared to what PREM does. That is exactly what we observe: for Crust2.0, the fit to the data is always poor when the path is dominantly oceanic or continental. The success of the CUB model in almost all numerical simulations is due to the fact that both crust and upper mantle 3-D structures are taken into account. For some specific paths, the right balance in the mixing of oceanic and continental paths can provide very good fits with PREM.

4.4 Conclusion

The comparison of seismograms computed in different Earth models clearly shows that the effect of the crust is large and non-linear. For pure oceanic paths or pure continental paths, the account of the crust tends to increase the residual delay times computed in PREM. This effect is well known (Montagner & Tanimoto 1991). When incorporating both crust and upper mantle 3-D structures, the time residuals are significantly improved, particularly at periods larger than 30 s. At short periods, the strong scattering effect of the crust gives rise to long coda waves which often makes our simple cross-correlation technique inefficient.

Our brief analysis puts in evidence how important the account for the shallow structures is to correctly model the seismic wave propagation at the regional scale. Previous studies (Komatitsch et al. 2002; Bozdag & Trampert 2008) used the SEM at the global scale and already noted the great influence of the shallow parts, but they did not honour the topography of the Moho. Because the computation cost is less important at the regional scale and because of the ability of RegSEM to take into account any Moho topography, this is overcome here. To exhibit the importance of honouring the Moho, we compare two seismograms computed at the same station (ANTO) in two different meshes of the CUB model (Fig. 14). In the first mesh, the Moho is fully honoured. In the second mesh, the elements at the top of the chunk are all 50 km thick, so the Moho is not honoured at all. The differences between the two traces are clear: a shift of half a period appears in the high frequency part of the signals (20-30 s), in addition with large discrepancies in amplitude (up to a factor of 2). At the global scale, similar differences are observed by Capdeville & Marigo (2008). These results illustrate the influence of the Moho on surface wavetrains and point out how crucial meshing its topography is to get accurate seismograms.

Figure 14

Rayleigh wave (vertical component) computed at station ANTO from two different simulations of a mid-Atlantic ridge earthquake in model CUB. In one simulation, the Moho is fully honoured (black line). In the other simulation, the Moho is interpolated by Lagrange polynomials of degree 4 defined in 50 km thick elements (red dashed line). The period range of the signals is 20-200 s.

Open in new tabDownload slide

Thanks to its flexibility, RegSEM is a suitable tool to assess the quality of regional tomographic models, as done at the global scale by Qin et al. (2009) and Bozdag & Trampert (2010). Nevertheless, the simulations presented here do not give precise conclusions on the quality of the CUB and Crust2.0 models. First, Crust2.0 is used with PREM although they might not be compatible as these two models have been obtained from two different data set and techniques. Second, the brutal smoothing applied to Crust2.0 and the CUB crust probably changes the effective elastic properties of the two models. A more rigorous smoothing technique, such as the one suggested by Capdeville et al. (2010) and Guillot et al. (2010), should be used to precisely assess the models.

5 Discussion

We showed that RegSEM is an efficient tool to compute seismic wavefields in geological media with possible complex geometries. Two particular cases are preferentially investigated. At the local scale, unstructured meshes of sedimentary basins (externally generated by CUBIT following the procedure suggested by Stupazzini et al. 2009) can be handled to study site effects. At longer scales (up to 90°), the influence on wave propagation of crustal and mantle 3-D structures with a possible Moho topography can be simulated using structured meshes generated by an internal routine. In both cases, the shallow structures yield important non-linear effects, even at wavelengths larger than the size of the heterogeneities. We showed that meshing these shallow structures is a crucial issue for a proper modelling of realistic wavefields.

In this work, the use of the SEM is limited to forward modelling. As mentioned in Section 1, the method can also be used to solve the inverse problem. A now classical and popular technique to achieve this goal consists in computing sensitivity kernels using the adjoint method (Tarantola 1984; Tromp et al. 2005; Fichtner et al. 2006a,b). This can be done with RegSEM, which provides kernels by simultaneously computing the adjoint wavefield and reconstructing the regular wavefield from time-reversed seismograms recorded at the boundaries (Gauthier et al. 1986). This process prevents from storing the whole regular wavefield, but it is rigorously valid in non-dissipative media only. Nevertheless, an artificial amplification can be introduced when reconstructing the regular wavefield in dissipative media (Tarantola 1988). Performing tomographic inversions using RegSEM is kept for future papers.

Another interesting application of our code is the computation of synthetic microseismic noise correlations. These data were introduced in seismology by Shapiro & Campillo (2004). Because most of the microseismic energy is contained in the 5-20 s period range and propagates as surface waves (Longuet-Higgins 1950), noise correlations are very sensitive to the shallow structure. Therefore, RegSEM is a suitable tool to perform realistic simulations (i.e. including the full complexity of wave propagation in 3-D media) of correlation waveforms. As discussed by Tromp et al. (2010), this is crucial to go beyond the ray theory classically used when inverting noise-correlation data. To mimick noise sources, the code imposes a random traction at the free surface of the chunk using the term in eq. 1. First correlations computed with RegSEM can be found in Stehly et al. (2011).

Eventhough RegSEM is a well advanced code, improvements are considered. First, there are simple options that could be easily added, such as introducing two layers of elements in thick crusts and allowing for external source time functions. Second, the creation of a realistic 3-D unstructured hexahedral mesh still remains a long and little flexible process. Moreover, the computation cost can be significantly high due to small or deformed elements resulting from the meshing of the interfaces. A de-refinement in depth could be set up to decrease the computation cost. Indeed, seismic velocities increase with depth so the wavelengths in the deeper part of the chunk are larger than in the shallow part. Therefore, keeping the same horizontal size of the elements everywhere in the medium yields an oversampling of the wavefield in the deeper part. This is particularly obvious in PREM (Fig. 2) when going into the lower mantle (i.e. crossing the 670 km discontinuity). To partially solve the problems associated with the mesh generation, we could think of a performant coupling between a tetrahedral SEM (Mercerat et al. 2006) for complex geometries and a classical hexahedral SEM for the rest of the domain. Refinement and de-refinement not only in space but also in time could also greatly help for the reduction in computation cost. More promising alternatives for the modelling of complex heterogeneities can be found in the field of mechanics. Homogeneization techniques, originally developed in material mechanics for the static case, have been recently applied to the propagation of seismic waves (Capdeville et al. 2010; Guillot et al. 2010). We plan to adapt RegSEM to this new technique in the future. Third, the PML could be improved. Figs 9 and 10 show the presence of a spurious wave in the medium. The amplitude of this wave is small but it could be even smaller if using filtering PML (Festa et al. 2005) or implementing the unsplit convolutional PML suggested by Martin & Komatitsch (2009). Finally, fluid could be incorporated in the code, using either a normal mode coupling (Capdeville et al. 2003) or the acoustic version of the SEM. All these improvements are in progress and will be available in future versions.

Acknowledgments

The authors want to thank Andreas Fichtner and Carl Tape for their insightful reviews which greatly helped in improving the manuscript. Interesting comments also emerged from discussions during different QUEST meetings. Remarks from Barbara Romanowicz, Diego Mercerat, Laurent Guillot, Ana Ferreira, Nikolaï Shapiro and Mark Panning have been appreciated. The authors thank Dimitri Komatitsch for providing the routines which compute internal forces and attenuation, as well as Marco Stupazzini for his precious help in the use of the CUBIT mesh generator. Many thanks also to Huaiyu Yuan, Laurent Stehly, Yilong Qin, Aimin Cao and Murat Erduran for using RegSEM and reporting some bugs. The experience of Geneviève Moguilny and Patrick Stoclet from the SCP (Service de Calcul Parallèle) of the IPGP has also been helpful for developing the code. The simulations presented in this paper have been performed using the Berkeley Seismological Laboratory cluster and the SCP cluster. ANR mémé partially supported this work.

References

Author notes