# A nodal discontinuous Galerkin approach to 3-D viscoelastic wave propagation in complex geological media

A nodal discontinuous Galerkin approach to 3-D viscoelastic wave propagation in complex... Summary A nodal discontinuous Galerkin (NDG) approach is developed and implemented for the computation of viscoelastic wavefields in complex geological media. The NDG approach combines unstructured tetrahedral meshes with an element-wise, high-order spatial interpolation of the wavefield based on Lagrange polynomials. Numerical fluxes are computed from an exact solution of the heterogeneous Riemann problem. Our implementation offers capabilities for modelling viscoelastic wave propagation in 1-D, 2-D and 3-D settings of very different spatial scale with little logistical overhead. It allows the import of external tetrahedral meshes provided by independent meshing software and can be run in a parallel computing environment. Computation of adjoint wavefields and an interface for the computation of waveform sensitivity kernels are offered. The method is validated in 2-D and 3-D by comparison to analytical solutions and results from a spectral element method. The capabilities of the NDG method are demonstrated through a 3-D example case taken from tunnel seismics which considers high-frequency elastic wave propagation around a curved underground tunnel cutting through inclined and faulted sedimentary strata. The NDG method was coded into the open-source software package NEXD and is available from GitHub. Numerical modelling, Computational seismology, Wave propagation 1 INTRODUCTION In geophysics, numerical simulations are a key tool for understanding physical phenomena taking place on and inside the earth. As they usually make predictions about observable quantities, they also play an essential role in inferring the current state and the physical properties of the earth’s interior. A phenomenon of particular importance in this context is seismic waves which travel through the entire globe and allow a unique look into the deep earth. Seismic waves are also useful on smaller scales, for example to find and characterize oil and gas reservoirs or to identify unknown objects in front of a newly excavated tunnel. Even on laboratory scale, elastic waves can help to increase knowledge about material behaviour, for example, to understand the mechanisms of crack propagation in rocks. In computational seismology and seismics, the need for high performance numerical simulations is continuously growing. Generally, in order to conduct seismic simulations, several different methods can be used. Overviews are given in the textbooks by Moczo et al. (2014) and Igel (2016). The probably most common method is the finite difference (FD) method (Virieux 1986; Moczo et al. 2007; Bohlen 2002) in which spatial derivatives in the governing differential equations are approximated by finite differences on an equally spaced grid (Robertsson et al. 1994). The strength of FD method is the simplicity of its numerical representation, which is however accompanied by numerical dispersion, especially for surface waves, and inaccurate representation of internal and external boundaries. Grid adaption, taking into account varying model properties, is hardly possible and a compromise between accuracy and computational cost must be carefully worked out for every simulation. Spectral and pseudospectral methods are used to simulate seismic wave propagation with high order accuracy (Tessmer & Kosloff 1994; Furumura et al. 1998). These methods rely on global basis functions such as Chebychev- or Legendre polynomials leading to simulations with little numerical dispersion. Only few points per wavelength are needed to achieve sufficient accuracy. However, the choice of global basis functions restricts the spectral methods to very smooth models, so that strong material contrasts and discontinuities are not represented adequately. Boundary Integral Methods (Sánchez-Sesma et al. 1993; Chen 2007) have been successfully used to simulate seismic wave propagation in domains with layers and interfaces, but are restricted to a limited amount of layers and are not well suited for 3-D problems as the involved matrices become very large and ill-conditioned. The finite element method (FEM) was used by Bao et al. (1998) and Marfurt (1984) to simulate seismic wave propagation. The FEM is based on the idea to subdivide the computational domain into elements, on which the equations can be solved. Albeit handling complex geometries very well, the solution of the wave equation is represented by low order polynomials resulting in numerical dispersion. For an accurate simulation, very large matrix systems have to be solved iteratively leading to high numerical costs with poor parallelization properties. Most commercial simulation software is based on FEM and is thus not well suited for seismic wave propagation problems. The idea of combining the spectral approximation properties with the flexibility of the FEM led to the development of the spectral element method (SEM; Patera 1984; Seriani & Priolo 1994) which was adapted to seismic wave propagation by Komatitsch & Vilotte (1998). The SEM uses a high order approximation of the wavefield inside the elements based on Lagrange polynomials. For the SEM, a variational formulation of the governing equations is used and a diagonal mass matrix can be constructed when taking Gauss–Lobatto–Legendre quadrature points within Cartesian elements (quadrilateral and hexahedral elements) for numerical integration. This leads to a fully explicit scheme making the SEM computationally very efficient and allowing for a simple parallelization via domain decomposition. The SEM is now widely used in the seismic community, especially for global and regional earthquake simulations (Komatitsch et al. 1999, 2003; Peter et al. 2011). Several free software implementations of the SEM exist for this purpose such as SPECFEM, SES3D and RegSEM (Geodynamics.org 2009; Fichtner et al. 2009; Cupillard et al. 2012). Unfortunately, the generation of well-structured hexahedral meshes can be very cumbersome for complex geological environments (Tautges 2001). The demand for a high order approximation of the seismic wavefield paired with the geometrical flexibility of an unstructured tetrahedral mesh led to the adaption of the discontinuous Galerkin (DG) method to seismic problems. The DG method was first proposed by Reed & Hill (1973) for solving the linear neutron transport equation. An extension to several research fields such as fluid dynamics (Bassi & Rebay 1997) and electrodynamic simulations (Cockburn et al. 2004) shows the capability of the DG method for simulating physical problems. An extensive mathematical research on the method was done by Cockburn et al. (2000) and Hesthaven & Warburton (2008). The DG method can be seen as a combination of the finite volume method (FVM; LeVeque 2002) and the SEM to solve partial differential equations with high accuracy. In contrast to the SEM and FEM, the variational formulation of the governing equations is done locally for every element and not globally. The solution obtained with the DG is allowed to be discontinuous across element boundaries. As in finite volume methods, this local formulation of the DG leads to the concept of numerical fluxes for exchanging information between adjacent elements. This results in an explicit semi-discrete numerical scheme which is very well suited for parallelization with domain decomposition techniques as only fluxes have to be communicated at element boundaries. In the DG method, one can choose a modal approach where the solution is represented by the expansion coefficients of polynomial basis functions, or a nodal approach where the solution is represented by its values at the anchor points of Lagrange polynomials used to approximate the solution within each element. For seismic wave propagation, the modal DG method was first introduced by Käser & Dumbser (2006) together with the concept of Riemann fluxes and arbitrary high order derivatives (ADER-DG). Here, time integration is coupled with the spatial resolution resulting in the same high order approximation of the solution in space and time. This approach results in a rather complicated numerical scheme due to the use of the Cauchy–Kovalewski procedure for the time integration. Intensive research on extending the ADER-DG to wave propagation problems in viscoelastic (Käser et al. 2007) and anisotropic media (de la Puente 2008) was conducted. It was also adapted to the study of rupture processes (Pelties et al. 2012). From this research sprang the community code SeisSol which became open-source in 2015 (https://github.com/SeisSol/SeisSol). The following work on SeisSol focused on increasing the numerical efficiency and adapting it to extremely large-scale seismic simulations and dynamic earthquake rupturing problems (Heinecke et al. 2014). Nodal high-order DG methods for seismic wave propagation were first developed by Delcourte et al. (2009) and Etienne et al. (2010) for the isotropic, elastic case. They use a pseudo-conservative formulation of the velocity–stress equations where the elastic material properties are assembled in a single diagonal matrix. Numerical fluxes are calculated by averaging velocity–stress values of adjacent elements. The solution is advanced in time using a leap-frog time integration. Delcourte et al. (2009) treated the 2-D case whereas Etienne et al. (2010) extended the approach to a low-order nodal DG (NDG) method for 3-D elastic wave propagation that exploits the hp-adaptivity of the DG method. The latter approach was benchmarked against other codes implementing FD, spectral element and pseudo-spectral methods in an earthquake ground motion application (Maufroy et al. 2015). Besides that, Wilcox et al. (2010) studied wave propagation for elastic-acoustic media within the DG framework based on hexahedral elements in the context of global earthquake simulation and Li (2011) used the NDG method to conduct 2-D simulations for non-linear elastic wave phenomena. One peculiarity of most DG methods (including SeisSol) is that material properties have to be constant within the elements. In case of continuously varying material properties, this assumption may lead to inaccuracies if the dominating wavelength approaches the element size because then the continuous structure is no longer properly approximated by a piecewise-constant representation. Furthermore, in extremely complicated media where discontinuities of material properties may cross an element, larger errors of the classical DG approach are to be expected. This problem was solved by Mercerat & Glinsky (2015) who extended the NDG approaches of Delcourte et al. (2009) and Etienne et al. (2010) to the case of heterogeneous elements for the 2-D, purely elastic case and demonstrated that, given equal element size, their approach yields more accurate results than the constant-element approaches. Unfortunately, it is unclear whether the pseudo-conservative form of the velocity–stress equations of Etienne et al. (2010) can be extended to the viscoelastic and anisotropic case. In this paper, we present a high-order nodal approach to the DG method based on a general mathematical treatment by Hesthaven & Warburton (2008) to model seismic wave propagation in 2-D and 3-D viscoelastic media with complex geological structures. Discretization is done using triangular (2-D) or tetrahedral (3-D) meshes. Instead of the pseudo-conservative form of Etienne et al. (2010), we use a classical, first-order in time velocity–stress formulation which allows the treatment of viscoelastic and also anisotropic media. We provide a justification for the necessity of introducing numerical fluxes into the computational scheme. We give explicit expressions for these fluxes, derived from an exact solution of the Riemann problem for the viscoelastic wave equation and expressed in terms of waves propagating away from the discontinuity. In contrast to the approach by Käser & Dumbser (2006) implemented in SeisSol, our fluxes honour differing material properties in adjacent elements. This way of treating the numerical fluxes allows for an easy extension of the method to cases where certain, potentially time-dependent discontinuities of particle velocities or stresses are prescribed at element boundaries, for example in the presence of fractures or during dynamic rupturing. The NDG approach leads to a simple, efficient and extendable numerical code with little logistical overhead. It is very intuitive and convenient as the solution vector directly represents the values of particle velocities and stresses at the nodal points. Concerning the evaluation of surface integrals on element faces required for calculating numerical fluxes, NDG is more efficient than modal DG as only information of nodes on the element faces is required. Our NDG code is capable of solving a wide range of seismic problems, from small-scale field experiments to large-scale earth models. Extensions of the code to unsaturated poroelastic (Boxberg et al. 2017) and fractured media are being worked on. Our implementation of the method has an interface to the full waveform inversion code ASKI (Schumacher et al. 2016) and offers the opportunity to compute adjoint wavefields. The paper is structured as follows: We briefly introduce the theoretical background of our approach and derive a numerical scheme suited for efficient parallel implementation. Some detail is dedicated to the derivation of the heterogeneous numerical fluxes for elastic and viscoelastic media. The general structure of the implemented algorithm is lined out and our newly developed program package NEXD is presented. The code is validated in 2-D by comparing its solution to Lamb’s problem with a semi-analytical one and by carrying out a p-convergence test. In addition, we benchmark the 2-D code against an analytical solution for a viscoelastic full space. In 3-D, we compare to results of SPECFEM (Geodynamics.org 2009) for a two-layer half-space. The potential and functionality of our NDG approach is demonstrated by means of a complex example taken from tunnel seismics where we consider 3-D propagation of elastic waves around a curved tunnel that cuts through faulted sedimentary strata. 2 A NODAL DISCONTINUOUS GALERKIN SCHEME FOR ELASTIC WAVE PROPAGATION The following section describes the theoretical basis of the numerical formulations developed and implemented in this paper. The main part deals with the adaption of the NDG method (Hesthaven & Warburton 2008) to the simulation of 3-D elastic wave propagation on unstructured tetrahedral meshes. 2.1 Velocity–stress formulation of the elastic wave equation in 3-D Following Virieux (1986) and Dumbser & Käser (2006) the elastic wave equation can be written as a system of first order hyperbolic equations. For the 3-D, elastic and isotropic case the system reads:   \begin{eqnarray} \frac{\partial }{\partial t } \sigma _{xx} -(\lambda +2\mu )\frac{\partial }{\partial x} u -\lambda \frac{\partial }{\partial y} v -\lambda \frac{\partial }{\partial z} w &=& s_{1}, \nonumber \\ \frac{\partial }{\partial t } \sigma _{yy} -\lambda \frac{\partial }{\partial x} u -(\lambda +2\mu ) \frac{\partial }{\partial y} v -\lambda \frac{\partial }{\partial z} w &=& s_{2}, \nonumber \\ \frac{\partial }{\partial t } \sigma _{zz} -\lambda \frac{\partial }{\partial x} u -\lambda \frac{\partial }{\partial y} v -(\lambda +2\mu ) \frac{\partial }{\partial z} w &=& s_{3}, \nonumber \\ \frac{\partial }{\partial t} \sigma _{xy}-\mu \left(\frac{\partial }{\partial x} v +\frac{\partial }{\partial y} u\right) &=& s_{4}, \nonumber \\ \frac{\partial }{\partial t} \sigma _{yz}-\mu \left(\frac{\partial }{\partial z} v +\frac{\partial }{\partial y} w\right) &=& s_{5}, \nonumber \\ \frac{\partial }{\partial t} \sigma _{xz}-\mu \left(\frac{\partial }{\partial x} w +\frac{\partial }{\partial z} u\right) &=& s_{6}, \nonumber \\ \frac{\partial }{\partial t} u - \frac{1}{\rho }\frac{\partial }{\partial x} \sigma _{xx} - \frac{1}{\rho }\frac{\partial }{\partial y} \sigma _{xy} - \frac{1}{\rho }\frac{\partial }{\partial z} \sigma _{xz} &=& s_{7} \nonumber ,\\ \frac{\partial }{\partial t} v - \frac{1}{\rho }\frac{\partial }{\partial x} \sigma _{xy} - \frac{1}{\rho }\frac{\partial }{\partial y} \sigma _{yy} - \frac{1}{\rho }\frac{\partial }{\partial z} \sigma _{yz} &=& s_{8}, \nonumber \\ \frac{\partial }{\partial t} w - \frac{1}{\rho }\frac{\partial }{\partial x} \sigma _{xz} - \frac{1}{\rho }\frac{\partial }{\partial y} \sigma _{yz} - \frac{1}{\rho }\frac{\partial }{\partial z} \sigma _{zz} &=& s_{9}, \end{eqnarray} (1)where σlk represent Cartesian stress components and u, v and w denote the x-, y- and z-components of the particle velocities. We use here Lamé’s elastic moduli λ and μ. Density is denoted by ρ. Note that all quantities are functions of space and time. A single force source can be defined by s = (0, 0, 0, 0, 0, 0, ax, ay, az)T, where the components of the acceleration ax, ay and az act in the last three components of the source vector. It is also possible to realize a moment tensor source which acts in the six stress components of the source vector s. In matrix-vector notation, this system of linear hyperbolic differential equations can be written in a more compact form as   \begin{eqnarray} \frac{\partial \mathbf {u}}{\partial t} + \mathbf {A} \frac{\partial \mathbf {u}}{\partial x} + \mathbf {B} \frac{\partial \mathbf {u} }{\partial y} + \mathbf {C} \frac{\partial \mathbf {u} }{\partial z} = \mathbf {s}. \end{eqnarray} (2)The vector   \begin{eqnarray*} \mathbf {u} = (\sigma _{xx}, \sigma _{yy},\sigma _{zz}, \sigma _{xy},\sigma _{yz},\sigma _{xz}, u, v, w)^{T} \nonumber \end{eqnarray*} contains the Nu = 9 velocity–stress variables, which all depend on time t and spatial coordinates x. The compact form contains the Jacobi matrices A = Aij(x), B = Bij(x) and C = Cij(x) describing the elastic properties of the medium (Dumbser & Käser 2006). In 3-D, the indices i, j range from 1 to 9. 2.2 Element-based integral forms of the velocity–stress equations The NDG method is built on an integral formulation of the velocity–stress equations. A computational domain, Ω, is defined encompassing the region of interest. It is subdivided into NE non-overlapping tetrahedral elements De, e = 1, …, NE  \begin{eqnarray} \Omega = \bigcup _{e=1}^{N_E} D^{e}. \end{eqnarray} (3)Within each element De the velocity–stress vector is represented by an interpolation formula of the form:   \begin{eqnarray} \mathbf {u}(\mathbf {x},t) = \sum _{j=1}^{N_p} \mathbf {u}(\mathbf {x}_j,t)\, l_j(\mathbf {x}), \quad \mathbf {x} \in D^e, \end{eqnarray} (4)where lj(x) is a multidimensional Lagrange polynominal anchored at Np interpolation points xi located within element De and on its faces. Since these polynomials satisfy the relation li(xj) = δij, the expansion coefficients are indeed identical to the values of the velocity–stress vector at the anchor points. Note also that no a priori assumptions are made on the continuity of the velocity–stress vector across adjacent elements. This fact is one of the name-giving properties of the DG method. Anchor points on adjacent element faces always occur as double nodes principally allowing the definition of a discontinuous velocity–stress vector across these faces. The integral form is now constructed by forming the scalar product of eq. (2) with a set of Nu × Np test vectors and integrating over element e. A system of Nu × Np equations results that can be solved for the nodal values of the velocity–stress vector. Here, we choose test vectors of the form   \begin{eqnarray} \mathbf {w}_{ni} = \mathbf {e}_n\, l_i(\mathbf {x}), \quad 1 \le i \le N_p,\quad 1 \le n \le N_u, \end{eqnarray} (5)where the nth component of en is 1 and all other components are 0. This is equivalent to multiplying the velocity–stress equations by the matrix Ili(x) where I is the identity matrix of the Nu-dimensional velocity–stress space. For each element e, we obtain Np vector equations:   \begin{eqnarray} \int _{D^e} \frac{\partial \mathbf {u}}{\partial t}l_i\, \text{d}V &=& -\int _{D^e}\left(\mathbf {A}\frac{\partial \mathbf {u}}{\partial x} +\mathbf {B}\frac{\partial \mathbf {u}}{\partial y} +\mathbf {C}\frac{\partial \mathbf {u}}{\partial z}\right)l_i \text{d}V \nonumber \\ &&+ \int _{D^e} \mathbf {s} \ l_i \, \text{d}V, \quad 1 \le i \le N_p . \end{eqnarray} (6)An integration by parts leads to:   \begin{eqnarray} \int _{D^e} \frac{\partial \mathbf {u}}{\partial t} \, l_i \text{d}V &=& -\int _{\partial D^e}(\mathbf {f}_xn_x+\mathbf {f}_yn_y+\mathbf {f}_zn_z)\, l_i\, \text{d}\Sigma \nonumber \\ &&+\int _{D^e} \left(\frac{\partial (\mathbf {A}l_i)}{\partial x} +\frac{\partial (\mathbf {B}l_i)}{\partial y} +\frac{\partial (\mathbf {C}l_i)}{\partial z}\right)\mathbf {u}\, \text{d}V \nonumber \\ &&+ \int _{D^e} \mathbf {s}\, l_i\, \text{d}V, \quad 1 \le i \le N_p . \end{eqnarray} (7)Here, we have introduced the components of the unit normal vector on the element faces, nx, ny, nz, the surface element dΣ and the fluxes fx, fy and fz defined as   \begin{eqnarray} \mathbf {f}_x = \mathbf {A}\mathbf {u}+\mathbf {r}_x,\quad \mathbf {f}_y = \mathbf {B}\mathbf {u}+\mathbf {r}_y,\quad \mathbf {f}_z = \mathbf {C}\mathbf {u}+\mathbf {r}_z . \end{eqnarray} (8)While surface integrals over expressions such as Au are the mathematically correct terms appearing after integration by parts, we have added here additional flux terms, rx, ry and rz, whose necessity is justified and discussed later. A second useful integral form of the velocity–stress equations can be obtained after another integration by parts. Expressions of the type Au appear again under surface integrals but now cancel the ones appearing in eq. (8):   \begin{eqnarray} \int _{D^e} \frac{\partial \mathbf {u}}{\partial t} \, l_i \text{d}V &=& -\int _{\partial D^e}(\mathbf {r}_xn_x+\mathbf {r}_yn_y+\mathbf {r}_zn_z)\, l_i \text{d}\Sigma \nonumber \\ &&-\int _{D^e} \left(\mathbf {A}\frac{\partial \mathbf {u}}{\partial x} +\mathbf {B}\frac{\partial \mathbf {u}}{\partial y} +\mathbf {C}\frac{\partial \mathbf {u}}{\partial z}\right)\, l_i\, \text{d}V \nonumber \\ &&+ \int _{D^e} \mathbf {s}\, l_i\, \text{d}V, \quad 1 \le i \le N_p . \end{eqnarray} (9)For both integral forms of the velocity–stress equations, we require the Np vector equations (7) or (9) to be fulfilled for each element. 2.3 Two semi-discrete schemes for the NDG method To simplify the notation in the following, we will use the abbreviations   \begin{eqnarray} \mathbf {u}^j(t) = \mathbf {u}(\mathbf {x}_j,t) \quad \mbox{and}\quad \mathbf {s}^j(t) = \mathbf {s}(\mathbf {x}_j,t),\quad 1 \le j \le N_p \end{eqnarray} (10)for the value of the velocity–stress vector and source vector, respectively, at the jth anchor point. Moreover, we define normal fluxes defined by   \begin{eqnarray} \mathbf {f}_n &=& \mathbf {f}_xn_x+\mathbf {f}_yn_y+\mathbf {f}_zn_z\quad \mbox{and} \nonumber \\ \mathbf {r}_n &=& \mathbf {r}_xn_x+\mathbf {r}_yn_y+\mathbf {r}_zn_z . \end{eqnarray} (11)By substituting eq. (4) into eq. (7) and by further assuming that the Jacobi matrices are constant within each element, we obtain   \begin{eqnarray} &&{\sum _{j=1}^{N_p} \frac{\partial \mathbf {u}^j}{\partial t} \int _{D^e} l_j(\mathbf {x}) l_i(\mathbf {x}) \text{d}V = \mathbf {A} \sum _{j=1}^{N_p} \mathbf {u}^j \int _{D^e} l_j(\mathbf {x}) \frac{\partial l_i(\mathbf {x})}{\partial x}\text{d}V }\nonumber \\ &&+\mathbf {B} \sum _{j=1}^{N_p} \mathbf {u}^j \int _{D^e} l_j(\mathbf {x}) \frac{\partial l_i(\mathbf {x})}{\partial y}\text{d}V + \mathbf {C} \sum _{j=1}^{N_p} \mathbf {u}^j \int _{D^e} l_j(\mathbf {x}) \frac{\partial l_i(\mathbf {x})}{\partial z}\text{d}V \nonumber \\ &&- \int _{\partial D^e} \mathbf {f}_nl_i(\mathbf {x})\text{d}\Sigma + \sum _{j=1}^{N_p}\mathbf {s}^j \int _{D^e} l_j(\mathbf {x}) l_i(\mathbf {x}) \text{d}V,\quad 1 \le i \le N_p . \end{eqnarray} (12)Eq. (12) motivates the introduction of the mass matrix   \begin{eqnarray} M^e_{ij} = \int _{D^e} l_i(\mathbf {x}) l_j(\mathbf {x}) \text{d}V, \quad 1 \le i,j \le N_p \end{eqnarray} (13)and the stiffness matrices   \begin{eqnarray} S^e_{ij,x} &=& \int _{D^e} l_i(\mathbf {x}) \frac{\partial l_j(\mathbf {x})}{\partial x} \text{d}V, \quad 1 \le i,j \le N_p \nonumber \\ S^e_{ij,y} &=& \int _{D^e} l_i(\mathbf {x}) \frac{\partial l_j(\mathbf {x})}{\partial y} \text{d}V, \quad 1 \le i,j \le N_p \nonumber \\ S^e_{ij,z} &=& \int _{D^e} l_i(\mathbf {x}) \frac{\partial l_j(\mathbf {x})}{\partial z} \text{d}V, \quad 1 \le i,j \le N_p , \end{eqnarray} (14)by which we can rewrite eq. (12), using Einstein summation convention, as   \begin{eqnarray} M^e_{ij}\frac{\partial \mathbf {u}^j}{\partial t} &=& \mathbf {A} S^e_{ji,x}\mathbf {u}^j + \mathbf {B} S^e_{ji,y}\mathbf {u}^j+ \mathbf {C} S^e_{ji,z}\mathbf {u}^j + M^e_{ij}\mathbf {s}^j \nonumber \\ &&- \int _{\partial D^e}\mathbf {f}_n\, l_i(\mathbf {x})\, \text{d}\Sigma , \quad 1 \le i,j \le N_p. \end{eqnarray} (15) Similarly, by inserting the nodal representation (4) into the second integral form of the velocity–stress equations (9), we obtain   \begin{eqnarray} M^e_{ij}\frac{\partial \mathbf {u}^j}{\partial t} &=& -\mathbf {A} S^e_{ij,x}\mathbf {u}^j - \mathbf {B} S^e_{ij,y}\mathbf {u}^j- \mathbf {C} S^e_{ij,z}\mathbf {u}^j + M^e_{ij}\mathbf {s}^j \nonumber \\ &&- \int _{\partial D^e}\mathbf {r}_n\, l_i(\mathbf {x})\, \text{d}\Sigma , \quad 1 \le i,j \le N_p. \end{eqnarray} (16)Subtle but notable differences between the two semi-discrete schemes, eqs (15) and (16), are different signs in front of the Jacobi matrices, the different order of indices in the stiffness matrices and the occurrence of different fluxes. Given a source vector and initial values for velocity–stress, we can use the above system of equations for each element to calculate an update of the velocity–stress vector provided we know how to compute the additional normal flux rn. 2.4 Normal Riemann fluxes The additional normal flux needs to be introduced because, by construction of the DG method, the velocity–stress vector may be discontinuous on the faces of adjacent elements. These discontinuities act as sources of seismic waves that propagate away from the discontinuity in both directions. If the flux was omitted or evaluated from the values of the velocity–stress vector of one element alone, the discontinuity would go unnotified and the additional waves would be missing in the solution. Consider a situation of two adjacent element faces with normal pointing into the x-direction. The source term is for now assumed to vanish within these two elements. Furthermore, we assume that the velocity–stress vector varies smoothly inside the elements but is discontinuous across the contact surface. We can always decompose such a velocity–stress distribution into two parts: one that is continuous at the contact surface and one that is discontinuous there and independent of x in both elements. Interpreting eq. (9) as a recipe for updating the velocity–stress vector with time, we recognize that the volume integral on the right-hand side will generate temporal changes of velocity–stress for the continuous part. For the discontinuous part, however, there will be no contribution from the term A∂u/∂x which is the only one that can produce waves propagating perpendicular to the contact surface. Hence, omitting the additional normal flux rn, implies loosing the waves generated by the discontinuity. To obtain specific expressions for the additional normal flux, we consider the so-called normal Riemann problem (LeVeque 2002). There we seek a solution of the velocity–stress equations for purely 1-D elastic wave propagation normal to an interface situated at x = 0 at both sides of which the velocity–stress vector is initially constant but discontinuous across the interface. Material properties on each side of the interface are allowed to be different. Source-free elastic wave propagation on both sides is described by   \begin{eqnarray} &&{\frac{\partial \mathbf {u}^{(1,2)}}{\partial t} + \mathbf {A}^{(1,2)} \frac{\partial \mathbf {u}^{(1,2)}}{\partial x} = \mathbf {0} \quad \mbox{with}} \nonumber \\ &&\mathbf {u}^{(1)}(x,t=0) = \mathbf {U}^{(1)}\quad \mbox{and}\quad \mathbf {u}^{(2)}(x,t=0) = \mathbf {U}^{(2)} , \end{eqnarray} (17)where the superscripts 1 and 2 refer to the left (x < 0) and right side (x > 0) of the interface, respectively. The solution to this problem is derived in detail by LeVeque (2002) in the context of the finite volume method. For positive times, we expect that waves will be generated by the discontinuity that propagate away from the interface. Waves with negative velocity will propagate to the left and waves with positive velocity to the right. In 3-D, matrix A has 9 eigenvalues, representing wave speeds c1−9, which we order in a way that c1−3 are negative, c4−6 are zero and c7−9 are positive. To calculate c1−3, we use matrix A(1), to calculate c7−9, we use matrix A(2). The amplitudes of the associated waves are obtained by decomposing the negative velocity–stress jump in terms of the corresponding right eigenvectors of matrices A(1, 2), denoted by $$\mathbf {v}^{(1,2)}_k$$, taking the ones for medium 1 when the corresponding wave speed is negative and the ones for medium 2 when the corresponding wave speed is positive:   \begin{eqnarray} \mathbf {U}^{(1)}-\mathbf {U}^{(2)} = \sum _{k=1}^3\gamma _k\mathbf {v}^{(1)}_k +\sum _{k=4}^6\gamma _k\mathbf {v}_k+\sum _{k=7}^9\gamma _k\mathbf {v}^{(2)}_k . \end{eqnarray} (18)For eigenvectors 4–6 we do not specify which side they belong to because they are independent of the material properties. The coefficients γk specify the desired amplitudes of the propagating waves. They can be calculated by multiplying the above equation by the appropriate reciprocal vectors of $$\mathbf {v}^{(1,2)}_k$$. Velocity–stress for positive times can then be written as   \begin{eqnarray} \mathbf {u}^{(1)}(x,t>0) &=& \mathbf {U}^{(1)}-\sum _{k=1}^3\gamma _k\mathbf {v}^{(1)}_kH\left(t+\frac{x}{|c_k|}\right) \nonumber \\ \mathbf {u}^{(2)}(x,t>0) &=& \mathbf {U}^{(2)}+\sum _{k=7}^9\gamma _k\mathbf {v}^{(2)}_kH\left(t-\frac{x}{c_k}\right) , \end{eqnarray} (19)where H(t) is a Heaviside function. With this solution, the discontinuity at the interface is reduced to   \begin{eqnarray} \mathbf {u}^{(2)}(0,t)-\mathbf {u}^{(1)}(0,t) &=& -\sum _{k=4}^6\gamma _k\mathbf {v}_k . \end{eqnarray} (20)Due to the structure of the eigenvectors, the remaining discontinuity occurs on the components σyy, σzz and σyz of the velocity–stress vector which do not generate waves propagating away from the interface. The remaining stress components as well as the particle velocities are continuous as required for a welded contact of two materials at the interface. Explicit expressions for the normal flux are obtained by going back to eq. (9) but staying with the purely 1-D Riemann problem. Instead of a tetrahedral element as integration volume, we choose a rectangular box attached to the right side of the interface and extending a distance L into the positive x-direction. After a small time step Δt, the change of velocity–stress in the box is given by the sum of right-propagating waves in eq. (19). As discussed above, this change must be produced by the additional normal flux rn. To avoid contributions from the right face of the box, we choose its width L greater than the distance the fastest wave can travel within time Δt. Then, instead of eq. (9), we get:   \begin{eqnarray} \int _\Sigma \text{d}\Sigma \int _0^L \Delta \mathbf {u}\, l_i(x) \text{d}x = -\Delta t\int _\Sigma \mathbf {r}_n\, l_i(0) \text{d}\Sigma , \end{eqnarray} (21)where Σ denotes the interface with outward normal pointing into the negative x-direction and li(x) is considered as a 1-D Lagrange polynomial of order N. Inserting the right-propagating waves for Δu, we find   \begin{eqnarray} &&{\int _\Sigma \text{d}\Sigma \int _0^L \sum _{k=7}^9\gamma _k\mathbf {v}^{(2)}_kH\left(\Delta t-\frac{x}{c_k}\right)l_i(x) \text{d}x }\nonumber \\ &&{\quad =\int _\Sigma \text{d}\Sigma \sum _{k=7}^9\gamma _k\mathbf {v}^{(2)}_k\int _0^{c_k\Delta t} l_i(x) \text{d}x} \nonumber \\ &&{\quad \approx \int _\Sigma \text{d}\Sigma \sum _{k=7}^9\gamma _k\mathbf {v}^{(2)}_kc_k\Delta t\, l_i(0) = -\Delta t\int _\Sigma \mathbf {r}_nl_i(0)\text{d}\Sigma,} \end{eqnarray} (22)from which we conclude for the normal flux on the right-hand side of the interface with outward normal pointing to the left:   \begin{eqnarray} \mathbf {r}^2_n = -\sum _{k=7}^9 c_k\gamma _k\mathbf {v}^{(2)}_k . \end{eqnarray} (23)The normal flux on the left-hand side of the interface with outward normal pointing to the right is obtained by analogously considering a box on the left-hand side and inserting the left-propagating waves for the change of the velocity–stress vector. Then, we find   \begin{eqnarray} \mathbf {r}^1_n = +\sum _{k=1}^3 |c_k|\gamma _k\mathbf {v}^{(1)}_k . \end{eqnarray} (24)Note that the evaluation of the normal flux is easily extended to situations where a certain, eventually time-dependent, discontinuity of the velocity–stress vector is prescribed at the interface, for example during dynamic rupturing or in the presence of fractures. To compute the normal flux in practice for a given element face, the velocity stress vector is rotated into a coordinate frame whose x-axis points into the direction of the outward normal of the face. The flux rx is computed in this frame from the velocity–stress jump and then rotated back into the original coordinate frame. This procedure is greatly alleviated by the fact that for isotropic elastic media the Jacobi matrices are invariant under rotation. It is also noted here that the expressions for the fluxes in eqs (23) and (24) account for differing material properties in adjacent elements. The normal Riemann flux given by Dumbser & Käser (2006) is calculated using only material properties of the element under consideration. The concept of fluxes in the DG method may also be exploited to realize particular boundary conditions. For example, an absorbing boundary condition can be achieved by simply cancelling the incoming fluxes at the boundaries. Similarly, by mirroring the stresses via the fluxes at the boundaries, a free boundary condition can be constructed. 2.5 NDG on tetrahedral elements One key feature of the NDG is its local character. It is not necessary to invert global matrices thus gaining flexibility for numerical implementation, especially for the parallelization of the method. For the high order approximation of the method, it is necessary to find a set of anchor points xi, i = 1, …, Np for the multidimensional Lagrange polynomials. For this purpose, the elements are mapped to a standard reference element, shown in Fig. 1 for tetrahedra, by a coordinate transformation of the form $$\mathbf {x}({\boldsymbol\xi })$$ where $${\boldsymbol\xi }=(\xi ,\eta ,\zeta )$$ denotes the position vector in the reference frame. The anchor points within a tetrahedral element can be constructed by a technique called ‘warp and blending’ introduced by Warburton (2006) (Fig. 1). The number of points is dependent on the order N of the interpolation and is, for the 3-D case on a tetrahedral element (Hesthaven & Warburton 2008),   \begin{eqnarray} N_p = \frac{(N+1)(N+2)(N+3)}{6}. \end{eqnarray} (25)Using the standard reference element, the mass matrix may be calculated once in advance according to   \begin{eqnarray} M_{ij} = \int _{I} l_i({\boldsymbol\xi }) l_j({\boldsymbol\xi }) \text{d}V_{\mbox{ref}} = \frac{1}{J^e}\int _{D_k} l_i(\mathbf {x}) l_j(\mathbf {x}) \text{d}V\, \end{eqnarray} (26)for the physical element De and reference element I, where Je is the constant Jacobian of the element De describing the volume change during the transformation from the physical to the reference element. The stiffness matrices Sij,ξ, Sij,η and Sij,ζ in the reference frame are given by   \begin{eqnarray} S_{ij,\xi _l} = \int _\mathbf {I} l_i({\boldsymbol\xi }) \frac{\partial l_j({\boldsymbol\xi })}{\partial \xi _l} \text{d}V_{\mbox{ref}}, \end{eqnarray} (27)and are related to the ones defined in physical space by (using Einstein summation convention)   \begin{eqnarray} S^e_{ij,x_m} = J^e\frac{\partial \xi _l}{\partial x_m}\, S_{ij,\xi _l}. \end{eqnarray} (28)Evaluation of the stiffness matrices in the reference frame is conveniently done by expanding the derivatives of the Lagrange polynomials in terms of Lagrange polynomials:   \begin{eqnarray} \frac{\partial l_j({\boldsymbol\xi })}{\partial \xi _l} = \sum _{r=1}^{N_p}\frac{\partial l_j}{\partial \xi _l}({\boldsymbol\xi }_r)l_r({\boldsymbol\xi }), \end{eqnarray} (29)leading to   \begin{eqnarray} S_{ij,\xi _l} = \sum _{r=1}^{N_p}\frac{\partial l_j}{\partial \xi _l}({\boldsymbol\xi }_r) \int _\mathbf {I} l_i({\boldsymbol\xi }) l_r({\boldsymbol\xi })\text{d}V_{{\rm ref}} . \end{eqnarray} (30)Defining the derivative matrix   \begin{eqnarray} D_{rj,\xi _l} = \frac{\partial l_j}{\partial \xi _l}({\boldsymbol\xi }_r) \end{eqnarray} (31)the stiffness matrices can be expressed as (using summation convention)   \begin{eqnarray} S_{ij,\xi _l} = M_{ir}D_{rj,\xi _l}\quad \mbox{and}\quad S^e_{ij,x_m} = J^e\frac{\partial \xi _l}{\partial x_m}M_{ir}D_{rj,\xi _l} . \end{eqnarray} (32)The mass matrix and the derivative matrix on the reference element can be conveniently evaluated by expanding the Lagrange polynomials and their derivatives in terms of orthogonal Jacobi polynomials (Hesthaven & Warburton 2008). Figure 1. View largeDownload slide Standard reference element for a tetrahedral element. Coordinates in this reference frame vary between −1 and +1. Small grey squares show position of anchor points for an interpolation order of N = 4. Figure 1. View largeDownload slide Standard reference element for a tetrahedral element. Coordinates in this reference frame vary between −1 and +1. Small grey squares show position of anchor points for an interpolation order of N = 4. For quantities defined on an element face f such as the fluxes, we perform an expansion in terms of 2-D Lagrange polynomials $$l_s^{2D}(\mathbf {x})$$ anchored at the face nodes. These are identical to those 3-D Lagrange polynomials whose anchor points reside on the element face, whereas all other 3-D Lagrange polynomials identically vanish on the element face. The number of 2-D polynomials is given by Ns = (N + 1)(N + 2)/2. We write for the normal flux:   \begin{eqnarray} \mathbf {f}_{n,f}(\mathbf {x},t) = \sum _{s=1}^{N_s} \mathbf {f}_{n,f}^s(t)\, l_s^{2D}(\mathbf {x}), \, \, \mathbf {x} \in \partial D^e, \end{eqnarray} (33)leading to the appearance of a face mass matrix $$M^{f}_{is}$$ defined as   \begin{eqnarray} M^{f}_{is} = \int _{I_{F}} l_i({\boldsymbol\xi }) l^{2D}_s({\boldsymbol\xi }) \text{d}\Sigma _{\mbox{ref}} = \frac{1}{J^f}\int _{\partial D^e} l_i(\mathbf {x}) l^{2D}_s(\mathbf {x}) \text{d}\Sigma \end{eqnarray} (34)for every face f of the element De with associated reference face IF and J f denoting the Jacobian for the face. To obtain the final semi-discrete numerical scheme for the integral form of eq. (15), we multiply it by the inverse mass matrix of the reference element (eq. 26) leading to (using summation convention)   \begin{eqnarray} \frac{\partial \mathbf {u}^m(t)}{\partial t} &=& \left(\mathbf {A}\frac{\partial \xi _l}{\partial x}+\mathbf {B} \frac{\partial \xi _l}{\partial y}+ \mathbf {C}\frac{\partial \xi _l}{\partial z}\right)(M^{-1}D^T_{\xi _l}M^T)_{mj}\mathbf {u}^j(t) \nonumber \\ &&-\frac{J^f}{J^e}\sum _{f=1}^{4}(M^{-1}M^{f})_{ms}\mathbf {f}^s_{n,f}(t) +\mathbf {s}^m(t) . \end{eqnarray} (35)For the integral form eq. (16) we obtain the following semi-discrete numerical scheme (using summation convention):   \begin{eqnarray} \frac{\partial \mathbf {u}^m(t)}{\partial t} &=& -\left(\mathbf {A}\frac{\partial \xi _l}{\partial x}+\mathbf {B}\frac{\partial \xi _l}{\partial y}+ \mathbf {C}\frac{\partial \xi _l}{\partial z}\right)(D_{\xi _l})_{mj}\mathbf {u}^j(t) \nonumber \\ &&-\frac{J^f}{J^e}\sum _{f=1}^{4}(M^{-1}M^{f})_{ms}\mathbf {r}^s_{n,f}(t)+\mathbf {s}^m(t) , \end{eqnarray} (36)where the $$\mathbf {r}^s_{n,f}$$ are now nodal coefficients of the normal Riemann flux with respect to 2-D Lagrange polynomials defined on the element faces. The two equations above represent the numerical schemes of the NDG method. They can be efficiently embedded into a numerical implementation and solved for the velocity–stress vector u. A similar equation can be derived for the 2-D case, where the system reduces to a system of five equations, to be solved for a velocity–stress vector with components u = (σxx, σyy, σxy, u, v)T. 2.6 Rules for spatial discretization The choice of element size, that is, the edge length of an element, is governed by the polynomial order of the Lagrange expansion, N, and the smallest significant wavelength, λmin, to be propagated stably across the model or parts of the model. According to our experience, we need 5 anchor points per minimum wavelength. Since the edge length of an element, ℓ, is about N times the spacing between anchor points, we get   $$\ell \le N\frac{\lambda _\mathrm{min}}{5} .$$ (37)The smallest wavelength is calculated from the highest significant frequency fmax and the smallest propagation velocity, vmin as   $$\lambda _\mathrm{min} = \frac{v_\mathrm{min}}{f_\mathrm{max}},$$ (38)leading finally to   $$\ell \le \frac{N}{5}\frac{v_\mathrm{min}}{f_\mathrm{max}}.$$ (39)The general procedure is to select one or several (for model parts) guiding values for element size following eq. (39) according to which the meshing software constructs a mesh. This mesh is then checked for stability again using eq. (39) but now solved for the highest frequency for which a stable simulation is possible:   $$f_\mathrm{stable} =\frac{N}{5}\frac{v_\mathrm{min}}{\ell }.$$ (40)The mesh passes the check if fstable ≤ fmax for all elements. If some elements do not satisfy this condition, a new mesh with a smaller guiding value for element size is constructed and the check is repeated. 2.7 Extension to anelasticity To incorporate anelasticity into our scheme, we follow the approach of Käser et al. (2007) with an adaptation to our way of flux computation. Based on original work by Emmerich & Korn (1987), it is assumed that the viscoelastic rheology can be described by a generalized Maxwell body mechanically represented by a spring in parallel to a series of Maxwell bodies each of which is made up of a spring and a dashpot. The stress–strain relation then generalizes to the form   \begin{eqnarray} \sigma _i(t) = M^u_{ij}\left(\epsilon _j(t)-\sum _{r=1}^{N_r}Y_r\omega _r\int _{-\infty }^t\epsilon _j(\tau )e^{-\omega _r(t-\tau )}{\rm d}\tau \right), \end{eqnarray} (41)where the σi stand for the six stresses, the εj for the corresponding strains, and the $$M^u_{ij}$$ for unrelaxed elastic moduli. Each of the Nr Maxwell bodies is characterized by a relaxation frequency ωr and an anelasticity coefficient Yr. For an isotropic medium, all moduli can be expressed through the Lamé constants λ and μ. Integration of this relation into the velocity–stress formulation is done by defining new material independent memory variables related to the anelastic stresses (Moczo et al. 2007):   \begin{eqnarray} \phi ^r_i(t) = \omega _r\frac{\partial }{\partial t}\int _{-\infty }^t\epsilon _i(\tau )e^{-\omega _r(t-\tau )}\text{d}\tau , \quad 1\le r\le N_r . \end{eqnarray} (42)According to eq. (41), all equations of the elastic velocity–stress formulation (eq. 1) containing time derivatives of stress will be extended by terms containing these new variables. The memory variables satisfy the differential equations   \begin{eqnarray} \frac{\partial \phi ^r_i(t)}{\partial t}+\omega _r\phi _i^r(t) = \omega _r\frac{\partial \epsilon _i}{\partial t} , \quad 1\le r\le N_r , \end{eqnarray} (43)which can be appended to the original velocity–stress eq. (1) because the strain-derivative on the right-hand side can be expressed by spatial derivatives of particle velocities. For each Maxwell body, six equations of this type are added to the system. Owing to the occurrence of spatial derivatives of particle velocities via the strain rates, the Jacobi matrices need to be extended to a dimension of 9 + 6Nr by adding 6Nr zero columns and 6Nr rows with nine entries each which are only non-zero if associated with particle velocity components. A modification of the Jacobi matrices affects the computation of the Riemann normal fluxes which were based (after rotation) on eigenvectors of the Jacobi matrix A. Here, we show that the coefficients of eqs. (23) and (24) describing amplitudes of waves propagating away from an element face remain unchanged and can be computed in exactly the same way as before. However, the eigenvectors change because these waves now also carry contributions in the new anelastic stress variables. The extended Jacobi matrix takes the following block structure:   \begin{eqnarray} \mathbf {A}^{\prime } = \left( \begin{array}{c{@}{\quad}c}\mathbf {A} & \mathbf {Z}_1 \\ \mathbf {P} & \mathbf {Z}_2 \end{array} \right) , \end{eqnarray} (44)where P is a 6Nr × 9 matrix, Z1 is a 9 × 6Nr zero matrix and Z2 is a 6Nr × 6Nr square zero matrix. Owing to the zero matrices, the extended Jacobi matrix has the same rank as the original one. Hence, all additional 6Nr eigenvalues must vanish. Now consider the case where the vector vk is an eigenvector of A with non-zero eigenvalue. This happens according to our numbering for index ranges 1−3 and 7−9. Then, there is an eigenvector of the extended Jacobi matrix, $$\mathbf {v}^{\prime }_k$$, satisfying $$\mathbf {A}^{\prime }\mathbf {v}^{\prime }_k = c_k\mathbf {v}^{\prime }_k$$, given by   \begin{eqnarray} \mathbf {v}^{\prime }_k = \left(\mathbf {v}_k,\mathbf {q}_k\right)^T\quad \mbox{with}\quad \mathbf {q}_k = \frac{1}{c_k}\mathbf {P}\mathbf {v}_k . \end{eqnarray} (45)If vk is an eigenvector of A but with a zero eigenvalue, then we can choose an eigenvector of A΄ of the form $$\mathbf {v}^{\prime }_k = (\mathbf {v_k},\mathbf {0})^T$$. This requires Pvk = 0, which is satisfied because the eigenvectors 4–6 of A vanish on the particle velocity components where P is non-zero. Finally, the remaining 6Nr eigenvectors with zero eigenvalue can be chosen as $$\mathbf {v}^{\prime }_k=(\mathbf {0},\mathbf {e}_k)^T$$ where ek is non-zero on the kth component only. An extended velocity–stress discontinuity is now decomposed as   \begin{eqnarray} \mathbf {U}^{\prime (1)}-\mathbf {U}^{\prime (2)} &=& \sum _{k=1}^3\gamma _k\left(\begin{array}{c}\mathbf {v}^{(1)}_k \\ \mathbf {q}^{(1)}_k\end{array}\right) +\sum _{k=4}^6\gamma _k\left(\begin{array}{c}\mathbf {v}_k \\ \mathbf {0}\end{array}\right) \nonumber \\ &&+\sum _{k=7}^9\gamma _k\left(\begin{array}{c}\mathbf {v}^{(2)}_k \\ \mathbf {q}^{(2)}_k\end{array}\right) +\sum _{k=10}^{9+6N_r}\gamma _k\left(\begin{array}{c}\mathbf {0} \\ \mathbf {e}_k\end{array}\right). \end{eqnarray} (46)However, considering only the first nine components, we arrive at exactly the same equation as eq. (18) from which the desired coefficients γk, 1 ≤ k ≤ 9 can be computed. Still, the associated eigenvectors changed. Hence, the extended normal flux on the right side of an interface now reads   \begin{eqnarray} \mathbf {r}^{2^{\prime }}_n = -\sum _{k=7}^9 c_k\gamma _k\left(\begin{array}{c}\mathbf {v}^{(2)}_k \\ \mathbf {q}^{(2)}_k\end{array}\right) . \end{eqnarray} (47) 2.8 Time discretization The time discretization is, in general, independent of the treatment of the spatial derivatives. Hence, established and well-proven methods for the time integrations can be used. For example, Komatitsch & Tromp (1999) use the classical second-order Newmark scheme to advance the SEM in time. For higher orders, the Runge–Kutta methods (Jameson et al. 1981) are available. Käser & Dumbser (2006) introduced the ADER-approach to the DG method for seismic problems by using the Cauchy–Kovalevsky theorem to solve the problem in the same order of approximation with respect to spatial and temporal variations. To advance the numerical schemes eqs (35) and (36) in time we apply a Total Variation Diminishing (TVD) Runge–Kutta method (Gottlieb & Shu 1998) which is third order accurate. It has three stages but needs only one additional array to store the velocity–stress–field for one physical integration step. The scheme has the form:   \begin{eqnarray} y^{n+1} & =& \frac{1}{3}y^n + \frac{2}{3}(z^{(2)} + \Delta t f(z^{(2)}))\quad \text{with} \nonumber \\ z^{(2)} &= &\frac{3}{4}y^n + \frac{1}{4}(z^{(1)} + \Delta t f(z^{(1)})),\quad \text{and} \nonumber \\ z^{(1)} &= &y^n + \Delta t f(y^n) . \end{eqnarray} (48)It computes the solution for the time step n + 1 going through two intermediate steps with the help of the auxiliary fields z(1) and z(2). Δt is the time step which has to satisfy the Courant–Friedrichs–Lewy criterion (Courant et al. 1928) according to which the velocity of the propagating waves must be smaller than the ratio of minimum point distance Δx divided by the time step Δt. The maximum time step is thus given by   \begin{eqnarray} \Delta t \le C_{\text{CFL}} \frac{\Delta x}{c_{\text{max}}}. \end{eqnarray} (49)Stability is ensured for CCFL ≤ 1. Typically CCFL is chosen smaller than 0.4 to run stable simulations on deformed meshes. The constant highly depends on the quality of the mesh and can be different for 2-D and 3-D simulations. Also, boundary conditions influence the value of the Courant constant. 3 NUMERICAL IMPLEMENTATION The NDG method described above was implemented into a software package called NEXD (Nodal Discontinuous Galerkin Finite Element in X Dimensions). As programming language, we use modern FORTRAN (FORTRAN95 and later) with object-oriented features. Parallelization is done by means of Open MPI (www.open-mpi.org). The code is organized in a way to allow a consistent work flow (Fig. 2). The model and the mesh can be generated with external meshing software. It is only required to provide plain text files containing the basic information about the mesh such as node numbering and connectivity, coordinates of vertices, boundary conditions and material parameters in the elements. A Python interface is available to use the Trelis/CUBIT software (www.csimsoft.com) for model and mesh generation. Three different program versions were developed, a 1-D, a 2-D and a 3-D one. Each program consists of three different main programs, the mesher, the solver and a post-processor to generate movie files. The mesher reads the mesh files and calculates all necessary information for the solver. It calculates the anchor points of the Lagrange polynomials and all necessary transformations, finds source and receiver positions and prepares the domain decomposition of the mesh with the help of the partitionizer ‘METIS’ (Karypis & Kumar 1999) for parallel simulations. All information is stored in a database available for each processor. Additional information for debugging and visualization purposes can be output on demand. Figure 2. View largeDownload slide The general procedure of a numerical simulation with NEXD. First, a mesh is created via a Python script and a meshing program such as Trelis/CUBIT. The ‘mesher’ reads in the mesh and prepares the forward calculation. If the simulation runs in parallel, the mesh is partitioned using the external library ‘Metis’. Database files with all required information about the mesh are created. These are read by the ‘solver’ which calculates the seismic wavefield. Parallelization of the solver is done with the help of MPI. Seismograms for displacement, velocity and acceleration are stored in plain ASCII files and databases are generated containing wavefield information. These files can be processed with the ‘movie’ program to create ‘VTK’ files to display videos of the propagating wavefield, for example, using the program ‘Paraview’. Figure 2. View largeDownload slide The general procedure of a numerical simulation with NEXD. First, a mesh is created via a Python script and a meshing program such as Trelis/CUBIT. The ‘mesher’ reads in the mesh and prepares the forward calculation. If the simulation runs in parallel, the mesh is partitioned using the external library ‘Metis’. Database files with all required information about the mesh are created. These are read by the ‘solver’ which calculates the seismic wavefield. Parallelization of the solver is done with the help of MPI. Seismograms for displacement, velocity and acceleration are stored in plain ASCII files and databases are generated containing wavefield information. These files can be processed with the ‘movie’ program to create ‘VTK’ files to display videos of the propagating wavefield, for example, using the program ‘Paraview’. The solver reads the database provided by the mesher and executes the forward simulation, generally in parallel mode using MPI. It controls the main time loop to advance the solution to the velocity–stress equations. For each time step, it iterates over all elements, deals with boundary conditions and calculates fluxes on element faces. Then, the velocity–stress vector is updated in time. At the end of the time loop, seismograms for displacement, velocity and acceleration are stored for every receiver position and component. On demand, snapshots of the wavefield are saved into databases which can be converted into VTK-files (Schroeder et al. 2006) for visualization with, for example ‘Paraview’ (Ayachit 2015). The code can handle different source time functions such as a Ricker or a Gaussian wavelet as well as externally provided source wavelets defined as time series. Both single force and moment tensor excitation are implemented. External velocity models defined on a regular grid can be mapped to an existing mesh by interpolation allowing a separation of meshing and assignment of material properties. Absorbing boundary conditions can be either realized by cancellation of incoming fluxes at the boundaries, or, much more efficiently, by using perfectly matching layers (PML). We have implemented a special variant of the PMLs in 2-D and 3-D, termed nearly perfectly matching layers (NPML) as described by Cummer (2003). Both 2-D and 3-D version of the code support the calculation of adjoint wavefields that may be used in the context of full waveform inversion. In addition, there is an interface to the full-waveform inversion code ASKI (Schumacher et al. 2016) for the computation of waveform sensitivity kernels. The code will be made available for download under GitHub (https://github.com/seismology-RUB, last accessed 1 December 2017). 4 NUMERICAL VALIDATION We validate our NDG approach for some test cases by comparing with analytical solutions, if available, and results from high-order SPECFEM simulations. 4.1 Lamb’s problem – 2-D elastic case The first test case is Lamb’s problem, the elastic response of a 2-D homogeneous half space to a single normal force at the free surface. The solution contains three main phases: a direct P-wave, a direct S-wave and a non-dispersive Rayleigh wave at the free surface. Lamb’s problem is a simple but challenging test to validate numerical codes especially for their accuracy with regard to dispersion. We compare solutions obtained with the NDG method to a semi-analytical reference solution, calculated with the program EX2DDIR (Berg et al. 1994). This program uses the Carniard-de Hoop (De Hoop 1960; Aki & Richards 1980) method to calculate the 2-D seismic response of a homogeneous half space. The model configuration and values for P-wave speed, S-wave speed and density of the homogeneous half-space are depicted in Fig. 3. Geometrical setting and material parameters exactly follow the configuration chosen by Käser & Dumbser (2006) and Komatitsch & Vilotte (1998). The free surface is inclined by an angle of α = 10° to avoid the Rayleigh wave propagating parallel to a coordinate axis. Due to the inclination, the solution is expected to be more sensitive to the numerical scheme and the implementation. The single force acting normal to the free surface resides at $$\mathbf {s} = (1720 \, \text{m},2303\, \text{m})$$ and two receivers are placed at the free surface at a distance of 990 m and 1706 m from the source. Receiver components are oriented normally and tangentially to the tilted free surface. As a source time function, a Ricker wavelet with a central frequency of fc = 14.5 Hz is used beginning at $$t_0 = -\frac{1.2}{fc}\, \text{s}$$ to ensure that the main peak of the wavelet is located at zero time. Figure 3. View largeDownload slide Model setup for computing the response of a homogeneous elastic half-space to a single normal force at the free surface, also known as Lamb’s problem. The free surface is inclined to avoid surface wave propagation parallel to a coordinate axis. Values of elastic wave speeds as well as position of source and receivers are shown together with the geometry of the model area. Figure 3. View largeDownload slide Model setup for computing the response of a homogeneous elastic half-space to a single normal force at the free surface, also known as Lamb’s problem. The free surface is inclined to avoid surface wave propagation parallel to a coordinate axis. Values of elastic wave speeds as well as position of source and receivers are shown together with the geometry of the model area. For the NDG method, the volume is discretized by triangular elements with maximum edge length of 45 m leading to a total of 17 700 elements for the mesh. A polynomial degree of five is used for the spatial resolution combined with a TVD Runge–Kutta scheme of third-order accuracy for time integration. The time step is Δt = 1.5 × 10−4 s and Nt = 10 000 time steps are calculated. On the left, right and bottom boundaries of the model, absorbing boundary conditions are used to mimic an infinitely wide half-space. A snapshot of the displacement field at t = 0.575 s is shown in Fig. 4. P- and S-waves propagate away from the source position followed by a Rayleigh wave with high amplitudes near the surface. A comparison of the resulting synthetic seismograms at the two receiver positions with the semi-analytical solution computed with EX2DDIR is depicted in Figs 5(a) and (b) for the normal component and Figs 5(c) and (d) for the tangential component. The seismograms fit very well. Maximum amplitude differences of less than 2 percent are visible in the high-amplitude Rayleigh wave train. This test validates the NDG method for accurate simulation of body waves and, in particular, surface waves. Figure 4. View largeDownload slide Snapshot of the absolute amplitude of the displacement field for Lamb’s problem after t = 0.575 s. The single force source is indicated as red dot, the two receivers as green dots. Three main phases are visible: a direct P-wave followed by an S-wave of smaller wave length due to the lower wave speed and a high amplitude surface wave near the free surface. Figure 4. View largeDownload slide Snapshot of the absolute amplitude of the displacement field for Lamb’s problem after t = 0.575 s. The single force source is indicated as red dot, the two receivers as green dots. Three main phases are visible: a direct P-wave followed by an S-wave of smaller wave length due to the lower wave speed and a high amplitude surface wave near the free surface. Figure 5. View largeDownload slide Displacement seismograms for Lamb’s problem for the normal and tangential component at receiver r1 (a, c) and receiver r2 (b, d). The analytical solution obtained with EX2DDIR is plotted in black. Results of the NDG simulation are plotted in grey. The difference is plotted as dashed line. A spatial order of N = 5 is used to obtain the NDG solution. The seismograms agree very well, validating the NDG for the 2-D case. The maximum relative difference occurring on the tangential component of receiver 1 is 2 per cent. No phase shift is notable for the two distant receivers, indicating that the NDG scheme calculates surface waves accurately. Figure 5. View largeDownload slide Displacement seismograms for Lamb’s problem for the normal and tangential component at receiver r1 (a, c) and receiver r2 (b, d). The analytical solution obtained with EX2DDIR is plotted in black. Results of the NDG simulation are plotted in grey. The difference is plotted as dashed line. A spatial order of N = 5 is used to obtain the NDG solution. The seismograms agree very well, validating the NDG for the 2-D case. The maximum relative difference occurring on the tangential component of receiver 1 is 2 per cent. No phase shift is notable for the two distant receivers, indicating that the NDG scheme calculates surface waves accurately. 4.2 P-convergence Test A different 2-D simulation in a homogeneous, elastic half-space was carried out to test the convergence of the numerical solution with increasing polynomial order of interpolation. Model setup, material properties and positions of a single force and a receiver are depicted in Fig. 6. At the top of the model, a free surface is assumed whereas on the other sides absorbing boundary conditions are applied. The force acts in horizontal (x)-direction and the source time function is a Ricker wavelet with centre frequency of fc = 600 Hz. The setup allows to observe the direct P- and S-waves on both receiver components and additionally converted and reflected waves coming from the free surface (Fig. 7). Given the seismic wave speeds and frequency spectrum of the source, the wave length of the S-wave at the centre frequency is λc ≃ 2.5 m while the minimum wavelength at the upper end of the frequency spectrum is λmin ≃ 1 m. Thus, the S-wave, in this setting, propagates a distance of about 20 wavelengths between source and receiver. The mesh consists of 21 606 elements leading to 648 180 degrees of freedom for the second-order test case (N = 2) and up to 7 129 980 degrees of freedom for the 10th-order calculation (N = 10) (Table 1) For comparison, reference seismograms are computed with a 12th-order SPECFEM simulation. Synthetic seismograms for both components displayed in Fig. 8 show small differences between the low-order (N = 3) NDG and the SPECFEM reference simulation which disappear for the 10th-order NDG solution. Figure 6. View largeDownload slide Model setup, material properties and positions of horizontal single force and receiver used for the p-convergence test. Figure 6. View largeDownload slide Model setup, material properties and positions of horizontal single force and receiver used for the p-convergence test. Figure 7. View largeDownload slide Snapshot of the absolute amplitude of the displacement field for the p-convergence test model. Time is T = 0.037 s and a polynomial order of N = 5 was used. Source position is marked by a red dot, receiver position by a green one. The mesh consists of 21 606 triangular elements. The direct P- and S-phases as well as the PP, PS and SP conversions from the free surface have already passed the receiver while the SS and PPS-phases are about to reach the receiver. Figure 7. View largeDownload slide Snapshot of the absolute amplitude of the displacement field for the p-convergence test model. Time is T = 0.037 s and a polynomial order of N = 5 was used. Source position is marked by a red dot, receiver position by a green one. The mesh consists of 21 606 triangular elements. The direct P- and S-phases as well as the PP, PS and SP conversions from the free surface have already passed the receiver while the SS and PPS-phases are about to reach the receiver. Figure 8. View largeDownload slide Seismograms of the displacement field for the (a) x- and (b) z-component for the p-convergence test. In both figures, the SPECFEM 12th-order reference solution is shown in black. The NDG solutions are green for the 3rd–order calculation and red for the 10th-order result. Whereas the third-order solution still has visible misfit, compared to the SPECFEM reference seismograms, no visible misfit can be detected with the N = 10 order NDG seismograms. Figure 8. View largeDownload slide Seismograms of the displacement field for the (a) x- and (b) z-component for the p-convergence test. In both figures, the SPECFEM 12th-order reference solution is shown in black. The NDG solutions are green for the 3rd–order calculation and red for the 10th-order result. Whereas the third-order solution still has visible misfit, compared to the SPECFEM reference seismograms, no visible misfit can be detected with the N = 10 order NDG seismograms. Table 1. Number of degrees of freedom for different polynomial orders of interpolation for a mesh with 21 606 triangular elements Order  Points per element  Degrees of freedom  2  6  648 180  3  10  1 080 300  4  15  1 620 450  5  21  2 268 630  6  28  3 024 840  8  45  4 861 350  10  66  7 129 980  Order  Points per element  Degrees of freedom  2  6  648 180  3  10  1 080 300  4  15  1 620 450  5  21  2 268 630  6  28  3 024 840  8  45  4 861 350  10  66  7 129 980  View Large To quantify the accuracy of NDG simulation with different polynomial orders, a misfit between the NDG and SPECFEM seismograms, expressed as L2 norm is calculated (Fig. 9). The misfit quickly decreases with polynomial order and converges to a minimum residual. Differences become very small from 5th order onwards. With regard to the numerical effort, it is sufficient to use fourth- or fifth-order polynomial interpolation for acceptably accurate simulations. Figure 9. View largeDownload slide L2 misfit between a 12-order SPECFEM simulation and NDG simulations of polynomial orders ranging from 2 to 10 for the p-convergence test case. As expected, the misfit decreases with increasing order and the numerical solution converges to the reference solution. Figure 9. View largeDownload slide L2 misfit between a 12-order SPECFEM simulation and NDG simulations of polynomial orders ranging from 2 to 10 for the p-convergence test case. As expected, the misfit decreases with increasing order and the numerical solution converges to the reference solution. 4.3 A 3-D simulation in a two-layer half-space This test simulation is performed with our 3-D implementation of the NDG method. The reference solution is, as before, computed with SPECFEM. The test shows the capability of the 3-D implementation and validates the NDG method with respect to 3-D simulations in discontinuous elastic media. Besides that, the potential of NDG to employ velocity adapted meshes is demonstrated. The considered model is a 3-D box of dimensions 100 m × 100 m× 50 m. Absorbing boundary conditions are used to simulate an infinitely wide model at all boundaries except at the top boundary, where a free surface boundary condition is assumed. The model is composed of a layer over a faster half-space. The top layer has a thickness of z = 12.5 m (Fig. 10). The material properties for this model are listed in Table 2. Figure 10. View largeDownload slide Two different meshes for a two layer simulation. Model dimensions are x = 100 m, y = 100 m and z = 50 m. The top layer has a thickness of z = 12.5 m and velocities vp = 3000 m s−1, vs = 1700 m s−1 and density ρ = 2000 kg m−3. The bottom layer is faster with material properties vp = 5500 m s−1, vs = 3400 m s−1 and ρ = 2000 kg m−3. The top surface is a free boundary. Absorbing boundary conditions are used for the other sides to mimic an infinite half-space. Model (a) has an element size of h = 2.5 m leading to 299 425 elements in total. The coarsened mesh in (b) has an element size of 5 m and 104 414 elements in total. Figure 10. View largeDownload slide Two different meshes for a two layer simulation. Model dimensions are x = 100 m, y = 100 m and z = 50 m. The top layer has a thickness of z = 12.5 m and velocities vp = 3000 m s−1, vs = 1700 m s−1 and density ρ = 2000 kg m−3. The bottom layer is faster with material properties vp = 5500 m s−1, vs = 3400 m s−1 and ρ = 2000 kg m−3. The top surface is a free boundary. Absorbing boundary conditions are used for the other sides to mimic an infinite half-space. Model (a) has an element size of h = 2.5 m leading to 299 425 elements in total. The coarsened mesh in (b) has an element size of 5 m and 104 414 elements in total. Table 2. Elastic properties of the two-layer half-space example in Section 4.3. Layer  vp (m s−1)  vs (m s−1)  ρ (kg m−3)  Yellow  3000  1700  2000  Green  5500  3400  2000  Layer  vp (m s−1)  vs (m s−1)  ρ (kg m−3)  Yellow  3000  1700  2000  Green  5500  3400  2000  View Large A vertical single force source is located in the half-space, at coordinates $$\mathbf {s} = (30\, \text{m}, 30\, \text{m}, 25\, \text{m})^T$$, 12.5 m below the bottom of the top layer. Note that in this example, the vertical coordinate z points bottom up. As source time function, a Ricker wavelet with dominant frequency of fc = 160 Hz is used. We perform two simulations with the NDG code: one on a mesh with approximately uniformly sized tetrahedral elements of an edge length of around h = 2.5 m leading to 299 425 elements and 94.3 × 106 degrees of freedom in total, and one on a mesh that is coarsened in the bottom layer to an edge length of around 5 m leading to 104 414 elements and 32.9 × 106 degrees of freedom in total. In the first simulation, we perform 16 000 time steps of length Δt = 6.415 × 10−6 s whereas in the second one, we run through 12 600 time steps of length Δt = 8.38 × 10−6 s. Polynomial order of interpolation is N = 4 in both simulations. Computation time of the first simulation is by a factor of 3.5 greater than that of the second one. For comparison, we run the same model with SPECFEM using hexahedral elements and polynomial order of 4. Here, 4500 time steps of length Δt = 2.34 × 10−5 s are performed. Element edge length is 2 m leading to 148 137 elements with 30 × 106 degrees of freedom in total. To record the wavefield, nine receivers are placed at the surface along a line parallel to the x-axis from $$\mathbf {r}_1=(10\, \text{m},50\, \text{m},50\, \text{m})^T$$ to $$\mathbf {r}_{9}=(90\, \text{m},50\, \text{m},50\, \text{m})^T$$. Synthetic seismograms obtained by all three simulations at the receivers are depicted in Fig. 11. No differences are visible on the scale of the figure. A quantitative calculation yields normalized rms-misfits for all seismograms of 0.2  per cent between SPECFEM and NDG with the uniform mesh, of 0.2  per cent between SPECFEM and NDG with the coarsened mesh and of 0.04  per cent between NDG with the uniform and NDG with the coarsened mesh. This test validates the 3-D-version of our DG method. Adapting the mesh to velocity changes in the model can significantly help to reduce computation times without loosing accuracy, and is easily realized with the DG method and tetrahedral elements. Figure 11. View largeDownload slide Seismic section of the z-component of particle velocity for the 3-D two-layer test. The line of nine receivers is located at the top of the model and has endpoint coordinates $$\mathbf {r}_1=(10\, \text{m},50\, \text{m},50\, \text{m})^T$$ and $$\mathbf {r}_{9}=(90\, \text{m},50\, \text{m},50\, \text{m})^T$$. The source is in the fast layer at $$\mathbf {s} = (30\, \text{m}, 30\, \text{m}, 25\, \text{m})^T$$. Black line: SEM reference simulation, dark grey line: NDG simulation with the coarsened mesh in the bottom layer. Light grey line: NDG simulation with about equal element size in both layers. For each receiver, the seismograms are normalized to the maximum of the SEM seismogram. Seismograms are shifted relative to each other. Figure 11. View largeDownload slide Seismic section of the z-component of particle velocity for the 3-D two-layer test. The line of nine receivers is located at the top of the model and has endpoint coordinates $$\mathbf {r}_1=(10\, \text{m},50\, \text{m},50\, \text{m})^T$$ and $$\mathbf {r}_{9}=(90\, \text{m},50\, \text{m},50\, \text{m})^T$$. The source is in the fast layer at $$\mathbf {s} = (30\, \text{m}, 30\, \text{m}, 25\, \text{m})^T$$. Black line: SEM reference simulation, dark grey line: NDG simulation with the coarsened mesh in the bottom layer. Light grey line: NDG simulation with about equal element size in both layers. For each receiver, the seismograms are normalized to the maximum of the SEM seismogram. Seismograms are shifted relative to each other. 4.4 Homogeneous full space—2-D anelastic case To validate our anelastic simulations we compare them to an analytical solution given by Carcione et al. (1988). They derived Green functions for a homogeneous 2-D full space in the frequency domain. Anelasticity is introduced by using complex-valued, frequency-dependent P- and S-wave velocities calculated from complex elastic moduli. A back transformation to the time domain results in the displacement field at the receiver position. For comparison with our DG method, we choose a homogeneous, quadratic model domain of 2000 m × 2000 m with a density of ρ = 2000 kg m−3. At a reference frequency of f = 50 Hz, the real part of P wave velocity is set to vp = 3000 m s−1 and that of S wave velocity to vs = 2000 m s−1. Anelasticity coefficients and relaxation frequencies are chosen in a way to achieve nearly constant values of the quality factors around Qκ = 30 and Qμ = 20. A vertical point force is located at the centre of the domain at x = y = 1000 m and a receiver resides at x = 1500 m, y = 500 m. Anelasticity is modelled, as described by Käser et al. (2007), using 3 Maxwell bodies leading to frequency dependent wave velocities and quality factors shown in Fig. 12. To evaluate the analytical expressions for the Green functions, we use exactly the same dispersion relation. The analytical Green function is convolved with the source time function used by Carcione et al. (1988) who chose the product of a Gaussian and a cosine function. Comparisons of the anelastic and elastic DG solutions with their analytical counterparts are shown in Fig. 13. For the elastic case, the values of P and S wave velocities at the reference frequency were used. In both cases, the numerical and analytical solution match very well. The relative difference between numerical and analytic solution at the extremal points is less than 10−3. Figure 12. View largeDownload slide Velocities for P and S waves (left) and quality factors, Qκ and Qμ, (right) versus frequency of the anelastic DG model. Figure 12. View largeDownload slide Velocities for P and S waves (left) and quality factors, Qκ and Qμ, (right) versus frequency of the anelastic DG model. Figure 13. View largeDownload slide Comparison of nodal DG simulation with an analytical solution after Carcione et al. (1988) for the anelastic and elastic cases. Traces are shifted relative to each other for better visibility. Figure 13. View largeDownload slide Comparison of nodal DG simulation with an analytical solution after Carcione et al. (1988) for the anelastic and elastic cases. Traces are shifted relative to each other for better visibility. 5 A COMPLEX NUMERICAL EXAMPLE This example case is designed to demonstrate the capability of our DG approach to model wave propagation in media exhibiting complicated geological structures and, in addition, containing bodies of non-trivial geometric forms. It is taken from related work in the field of tunnel seismics where seismic waves are used to image disturbing structures in front of the tunnel boring machine. The model is composed of a curved cylindrical tunnel embedded into three sedimentary layers which are offset relative to each other by two thrust faults (Fig. 14). The tunnel axis lies in the xy-plane. At the entry, the tunnel axis is parallel to the y-axis but later turns away to the right, that is, towards increasing x-values. The vertical coordinate, z, of the tunnel axis is constant. The interfaces between the sedimentary layers are slightly inclined relative to the plane of the tunnel axis and intersect the tunnel axis at some distance from the tunnel entry implying different material properties above and below the tunnel from there on. In addition, the thrust faults, inclined relative to the vertical, intersect the tunnel at some distance away from the entry. The elastic material parameters of the sedimentary layers are rather different and are listed in Table 3. Two high-velocity layers at the top and bottom of the model enclose a low-velocity layer in between. Figure 14. View largeDownload slide Model motivated by tunnel seismics containing a curved cylindrical tunnel drilled into three inclined sedimentary layers offset relative to each other by thrust faults. The tunnel is l = 50 m long and ends in the fast third layer (red) with material properties vp = 4650 m s−1, vs = 2700 m s−1 and ρ = 2600 kg m−3. This layer is covered by a slow layer (grey) with vp = 2000 m s−1, vs = 1200 m s−1 and ρ = 2350 kg m−3. The top layer (blue) is, again, a fast layer having a P-wave velocity of vp = 4330 m s−1 a S-wave velocity of vs = 2500 m s−1 and a density of ρ = 2500 kg m−3. The overall dimensions of the model are wx = 64 m, wy = 100 m and wz = 64 m. (a) Entire model box, (b) a horizontal section to visualize the track of the tunnel. The tetrahedral mesh contains Ne = 1 541 398 elements resulting in a simulation with 4.856 × 108 degrees of freedom. Figure 14. View largeDownload slide Model motivated by tunnel seismics containing a curved cylindrical tunnel drilled into three inclined sedimentary layers offset relative to each other by thrust faults. The tunnel is l = 50 m long and ends in the fast third layer (red) with material properties vp = 4650 m s−1, vs = 2700 m s−1 and ρ = 2600 kg m−3. This layer is covered by a slow layer (grey) with vp = 2000 m s−1, vs = 1200 m s−1 and ρ = 2350 kg m−3. The top layer (blue) is, again, a fast layer having a P-wave velocity of vp = 4330 m s−1 a S-wave velocity of vs = 2500 m s−1 and a density of ρ = 2500 kg m−3. The overall dimensions of the model are wx = 64 m, wy = 100 m and wz = 64 m. (a) Entire model box, (b) a horizontal section to visualize the track of the tunnel. The tetrahedral mesh contains Ne = 1 541 398 elements resulting in a simulation with 4.856 × 108 degrees of freedom. Table 3. Elastic properties of the sedimentary layers in the tunnel seismics model Layer  vp (m s−1)  vs (m s−1)  ρ (kg m−3)  Blue  4330  2500  2500  Grey  2000  1200  2350  Red  4650  2700  2600  Layer  vp (m s−1)  vs (m s−1)  ρ (kg m−3)  Blue  4330  2500  2500  Grey  2000  1200  2350  Red  4650  2700  2600  View Large The dimensions of the entire model box are 64 m × 100 m × 64 m. The tunnel is approximately 50 m long, ending in the fast bottom layer. Its diameter is 12 m. The model is meshed with tetrahedral elements for the NDG simulation with an edge size of h = 2 m for both upper and lower layer and with a size of h = 1 m for the middle slow layer. This choice allows stable simulations up to fc = 500 Hz. To save computational effort, some mesh coarsening was applied in the outer regions of the mesh. In total, Ne = 1 541 398 tetrahedral elements are required to discretize the model. Polynomial order is 4 and number of Lagrange anchor points Np is 35. With Nu = 9 field variables we end up with NpNuNe = 4.8 × 108 degrees of freedom for this simulation. Absorbing boundary conditions based on flux cancellation are applied on all boundaries except for the top boundary where a free surface is assumed. We also attempted to mesh this model with hexahedral elements, but failed to produce a mesh of acceptable quality. For this reason, a computation with SPECFEM could not be performed. Seismic waves in the model are excited by a single force residing at the sidewall of the tunnel and acting normally to it. Its precise location is $$\mathbf {x}_s = (34.53 \, \text{m}, 37.73 \, \text{m}, 32.00\, \text{m})$$. The source time function is a Ricker wavelet with centre frequency of 500 Hz. The time step is Δt = 3.414 × 10−6 s and Nt = 20 000 time steps were calculated. A line of 46 receivers is placed on the right tunnel sidewall (viewed from the tunnel entry). Receiver spacing is 1 m (Fig. 14). The simulation was carried out on 192 CPU cores. Different stages of the wavefield are displayed as snapshots in the xy-plane in Fig. 15. The source emits P- and S-waves and a high-amplitude surface wave attached to the tunnel sidewall. The wave speed discontinuity due to the first thrust fault is well recognized in Fig. 15(c). When the tunnel Rayleigh wave detaches from the tunnel front face, it converts into an S-wave which propagates away from the tunnel (Figs 15c–f). When hitting the second thrust fault (Fig. 15e), the curvature of the transmitted S-wave front changes and a reflection is produced (Fig. 15f). On arrival at the tunnel front face, this reflection converts into a tunnel Rayleigh wave (Fig. 15g) and runs along the entire tunnel side wall (Figs 15h–l). This scenario of a tunnel Rayleigh wave converting into an S-wave and back into a tunnel Rayleigh wave after being reflected has been described earlier by Bohlen et al. (2007). Figure 15. View largeDownload slide Snapshots of the wavefield for the curved tunnel model. The absolute value of the displacement field is shown in the xy−plane (see also Fig. 14b). The source is marked as red dot while the receivers are highlighted by green points. The time between the individual frames is Δt = 3.41 × 10−3 s. A single force normal to the tunnel sidewall (a) excites a high amplitude surface wave propagating along the tunnel (b) which is converted into an S-wave (c–e) propagating ahead of the tunnel. A reflection at the second thrust fault occurs (f) and a part of the signal propagates back to the tunnel where it is converted to a tunnel surface wave (g). This signal can be observed in the seismograms and propagates all along the free surface of the tunnel (h)–(l). Figure 15. View largeDownload slide Snapshots of the wavefield for the curved tunnel model. The absolute value of the displacement field is shown in the xy−plane (see also Fig. 14b). The source is marked as red dot while the receivers are highlighted by green points. The time between the individual frames is Δt = 3.41 × 10−3 s. A single force normal to the tunnel sidewall (a) excites a high amplitude surface wave propagating along the tunnel (b) which is converted into an S-wave (c–e) propagating ahead of the tunnel. A reflection at the second thrust fault occurs (f) and a part of the signal propagates back to the tunnel where it is converted to a tunnel surface wave (g). This signal can be observed in the seismograms and propagates all along the free surface of the tunnel (h)–(l). Synthetic seismograms of the y-component at the receivers are shown in Fig. 16. The traces are ordered from bottom to top according to their numbering which starts at 1 on the tunnel front face. The source resides close to receiver 11 and the first thrust fault lies between receivers 19 and 20. Close to the source (receivers 1–19), we recognize a superposition of the S-wave and the Rayleigh wave (R1). P-waves are very weak owing to the fact that the force acts normally to the tunnel sidewall. When arriving at the first thrust fault, a converted P-wave is excited (P) and the tunnel Rayleigh wave changes speed (R2). In addition, the tunnel Rayleigh wave is reflected from this discontinuity and propagates back to the tunnel front face (R4 and R3). The Rayleigh wave generated by the S-wave coming back from the second thrust fault is denoted by R5. It arrives first at the tunnel front face and then propagates all along the receiver line. Phase L is a reflection from the interface above the tunnel whereas ‘surf’ is a reflection from the free surface. Phase ‘AB’ is a reflection from the boundary that contains the tunnel entry owing to a not perfectly working absorbing boundary condition. Figure 16. View largeDownload slide Seismic section of the y-component (along axis at tunnel entry) of the displacement field. Receivers are distributed at the right-hand side of the tunnel sidewall (see Fig. 15, green dots, and Fig. 14b, white dots). Receiver 1 (bottom of section) resides at the tunnel front face. Receiver 46 (top of section) sits at the tunnel entry. The source is very close to receiver 11. ‘R1’ denotes the direct surface wave excited by a single force acting normal to the tunnel sidewall. ‘R2’ is its continuation after having crossed the first thrust fault which intersects the tunnel. ‘R3’ is a surface wave reflected from the thrust fault. ‘R4’ is a surface wave reflected from the tunnel front face. R1 also excites a P wave (‘P’) when hitting the thrust fault. ‘R5’ is the so-called RSSR-wave, a surface wave generated through conversion from an S-wave which is reflected from the second thrust fault ahead of the tunnel front face. ‘L’ denotes a reflection from the layer interface above the tunnel while ‘surf’ is the signal from the top free surface. At the tunnel entry, a nonphysical reflection owing to imperfect absorbing boundary conditions evolves and is marked as ‘AB’. Figure 16. View largeDownload slide Seismic section of the y-component (along axis at tunnel entry) of the displacement field. Receivers are distributed at the right-hand side of the tunnel sidewall (see Fig. 15, green dots, and Fig. 14b, white dots). Receiver 1 (bottom of section) resides at the tunnel front face. Receiver 46 (top of section) sits at the tunnel entry. The source is very close to receiver 11. ‘R1’ denotes the direct surface wave excited by a single force acting normal to the tunnel sidewall. ‘R2’ is its continuation after having crossed the first thrust fault which intersects the tunnel. ‘R3’ is a surface wave reflected from the thrust fault. ‘R4’ is a surface wave reflected from the tunnel front face. R1 also excites a P wave (‘P’) when hitting the thrust fault. ‘R5’ is the so-called RSSR-wave, a surface wave generated through conversion from an S-wave which is reflected from the second thrust fault ahead of the tunnel front face. ‘L’ denotes a reflection from the layer interface above the tunnel while ‘surf’ is the signal from the top free surface. At the tunnel entry, a nonphysical reflection owing to imperfect absorbing boundary conditions evolves and is marked as ‘AB’. 6 DISCUSSION The DG method is extraordinarily suitable for wave propagation simulations in complicated geological environments due to its ability to work with tetrahedral elements. With tetrahedra, a high-quality mesh is much easier to achieve than with the hexahedral elements used by the SEM. Creating hexahedral meshes for complicated models may take significantly more time and experience than a comparable tetrahedral mesh. The curved tunnel example demonstrates that it is extremely difficult to produce a high-quality hexahedral mesh, while for the same model a tetrahedral mesh can be constructed semi-automatically. Another big advantage of tetrahedral meshes is the fact that element sizes can be easily adapted to the physical properties of the medium without affecting the accuracy of the simulation. In this way, adaptive mesh construction can partially compensate the higher numerical costs of the DG method. This advantage pays off particularly well in complex geological situations, which often occur in unconsolidated sediments and are thus relevant for near surface seismic modelling. Nevertheless, it is worthwhile to carefully consider which method is used for a given simulation problem. In many cases, the SEM may be the method of choice. The realization of absorbing or radiation boundary conditions by just omitting incoming fluxes at boundary elements is not fully satisfactory especially for waves with grazing incidence. For this reason, NPML boundary conditions, a variation of PML conditions, are implemented. Spurious reflections from the absorbing boundary nearly vanish for the simulations using NPML. However, PML conditions are only weakly stable for the seismic wave equation (Fichtner et al. 2011) leading to cases where the solution grows exponentially with time. Especially large material contrast close to the absorbing boundary may favour instable behaviour. This problem is overcome by a switching mechanism. When using NPML boundary conditions, the total energy of the medium is monitored and if the energy, after a certain time, begins to grow, the NPML is switched to the classical, flux-controlled boundary conditions thereby saving the simulation. Exponential growth of energy in the solution is detected with an STA/LTA trigger, originally developed for earthquake detection. Viscoelastic attenuation is modelled using a series of parallel Maxwell bodies. The comparison of a viscoelastic simulations in 2-D with an analytical solution exhibits very good agreement. The numerical example also shows that three parallel Maxwell bodies are sufficient to obtain nearly constant quality factors and a log-linear dependency of wave velocities across a wide frequency range. Introducing more mechanisms is normally not required and increases the numerical costs significantly. The current implementation of the NDG method successfully passes several accuracy tests. For example, synthetic seismograms for Lamb’s problem very well match seismograms computed from an analytical solution. Further comparisons with SEM also show very good agreement. The success of these tests also validates the parallel implementation of the NDG code using MPI. A convergence test, where simulations are done with increasing spatial order, underline the successive reduction of misfit to a reference solution with increasing order. It turns out that a spatial order of N = 4 is a viable compromise between numerical effort and accuracy for the NDG method. Very high order simulations require a correspondingly small time step making the computations inefficient. The scaling of the parallel implementation with increasing numbers of processors (not described in this paper) could only be tested on a single machine with 48 CPUs owing to the heterogeneity of the available computer cluster. At least for this machine, the computation time decreases quasi-linearly with CPU number. We expect that this behaviour is also preserved for a larger number of CPUs. 7 CONCLUSIONS The NDG method is adapted to 1-D, 2-D and 3-D elastic and anelastic seismic wave propagation. Formulations of numerical fluxes for the elastic and anelastic case are given which are derived from an exact solution of the heterogeneous Riemann problem, and honour possibly differing material properties in adjacent elements. Absorbing boundary conditions relying on the Nearly Perfectly Matched Layer approach are adapted to the NDG as an alternative to flux based boundary conditions. All features are implemented in the software package NEXD allowing efficiently parallelized numerical simulations. The accuracy and implementation is successfully validated through different test cases such as Lamb’s problem by either comparing with analytical solutions or spectral element simulations. The capabilities of the method are demonstrated by a challenging example case derived from tunnel reconnaissance with complicated geological structures and internal bodies of non-trivial geometric form. The NDG method is particularly attractive for such scenarios due to its high order approximation of the seismic wavefield as well as its use of triangular and tetrahedral meshes for the model discretization. The code may be used as forward solver for full waveform inversion as it offers the calculation of adjoint wavefields and provides an interface to the computation of waveform sensitivity kernels. The code is available for download under GitHub (https://github.com/seismology-RUB, last accessed 1 December 2017). ACKNOWLEDGEMENTS This work was supported by the Collaborative Research Center 837 ‘Interaction models in mechanized tunneling’ funded by the Deutsche Forschungsgemeinschaft (DFG). REFERENCES Aki K., Richards P.G., 1980. Quantivative Seismology, Theorie and Methods , University Science Books. Ayachit U., 2015. The ParaView Guide: A Parallel Visualization Application , Kitware. Bao H., Bielak J., Ghattas O., Kallivokas L., O’Hallaron D., Shewchuk J., Xu J., 1998. Large-scale simulation of elastic wave propagation in heterogeneous media on parallel computers, Comput. Methods Appl. Mech. Eng. , 152(1–2), 85– 102. https://doi.org/10.1016/S0045-7825(97)00183-7 Google Scholar CrossRef Search ADS   Bassi F., Rebay S., 1997. A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier–Stokes equations, J. Comput. Phys. , 131( 2), 267– 279. https://doi.org/10.1006/jcph.1996.5572 Google Scholar CrossRef Search ADS   Berg P., If F., Nielsen P., Skovgaard O., 1994. Analytical reference solutions, Model. Earth Oil Explor. , 1, 421– 427. Bohlen T., 2002. Parallel 3-D viscoelastic finite difference seismic modelling, Comput. Geosci. , 28, 887– 899. https://doi.org/10.1016/S0098-3004(02)00006-7 Google Scholar CrossRef Search ADS   Bohlen T., Lorang U., Rabbel W., Mueller C., Giese R., Lueth S., Jetschny S., 2007. Rayleigh-to-shear wave conversion at the tunnel face—from 3D-FD modeling to ahead-of-drill exploration, Geophysics , 72( 6), T67– T79. https://doi.org/10.1190/1.2785978 Google Scholar CrossRef Search ADS   Boxberg M. S., Heuel J., Friederich W., 2017. A nodal discontinuous Galerkin solver for modeling seismic wave propagation in porous media, in Poromechanics VI , eds Vandamme M., Dangla P., Pereira J-M., Ghabezloo S., American Society of Civil Engineers. Carcione J. M., Kosloff D., Kosloff R., 1988. Wave propagation simulation in a linear viscoelastic medium, Geophys. J. Int. , 95( 3), 597– 611. https://doi.org/10.1111/j.1365-246X.1988.tb06706.x Google Scholar CrossRef Search ADS   Chen X.-F., 2007. Generation and propagation of seismic SH waves in multi-layered media with irregular interfaces, Adv. Geophys. , 48, 191– 264. Google Scholar CrossRef Search ADS   Cockburn B., Karniadakis G.E., Shu C.-W., 2000. The Development of Discontinuous Galerkin Methods , Springer. Google Scholar CrossRef Search ADS   Cockburn B., Li F., Shu C.-W., 2004. Locally divergence-free discontinuous Galerkin methods for the Maxwell equations, J. Comput. Phys. , 194( 2), 588– 610. https://doi.org/10.1016/j.jcp.2003.09.007 Google Scholar CrossRef Search ADS   Courant R., Friedrichs K., Lewy H., 1928. Über die partiellen Differenzengleichungen der mathematischen Physik, Math. Ann. , 100( 1), 32– 74. https://doi.org/10.1007/BF01448839 Google Scholar CrossRef Search ADS   Cummer S.A., 2003. A simple, nearly perfectly matched layer for general electromagnetic media, IEEE Microw. Wirel. Compon. Lett. , 13( 3), 128– 130. https://doi.org/10.1109/LMWC.2003.810124 Google Scholar CrossRef Search ADS   Cupillard P., Delavaud E., Burgos G., Festa G., Vilotte J.-P., Capdeville Y., Montagner J.-P., 2012. Regsem: a versatile code based on the spectral element method to compute seismic wave propagation at the regional scale, Geophys. J. Int. , 188( 3), 1203– 1220. https://doi.org/10.1111/j.1365-246X.2011.05311.x Google Scholar CrossRef Search ADS   De Hoop A., 1960. A modification of Cagniard’s method for solving seismic pulse problems, Appl. Sci. Res., Sec. B , 8( 1), 349– 356. https://doi.org/10.1007/BF02920068 Google Scholar CrossRef Search ADS   de la Puente J., 2008. Seismic Wave Propagation for Complex Rheologies , VDM Verlag Dr. Müller. Delcourte S., Fezoui L., Glinsky-Olivier N., 2009. A high-order discontinuous Galerkin method for the seismic wave propagation, in ESAIM: Proceedings , vol. 27, pp. 70– 89, EDP Sciences. Google Scholar CrossRef Search ADS   Dumbser M., Käser M., 2006. An arbitrary high order discontinuous Galerkin method for elastic waves on unstructured meshes II: the three-dimensional isotropic case, Geophys. J. Int. , 167( 1), 319– 336. https://doi.org/10.1111/j.1365-246X.2006.03120.x Google Scholar CrossRef Search ADS   Emmerich H., Korn M., 1987. Incorporation of attenuation into time-domain computations of seismic wave fields, Geophysics , 52( 9), 1252– 1264. https://doi.org/10.1190/1.1442386 Google Scholar CrossRef Search ADS   Etienne V., Chaljub E., Virieux J., Glinsky N., 2010. An hp-adaptive discontinuous Galerkin finite-element method for 3-D elastic wave modelling, Geophys. J. Int. , 183( 2), 941– 962. https://doi.org/10.1111/j.1365-246X.2010.04764.x Google Scholar CrossRef Search ADS   Fichtner A., Kennett B.L., Igel H., Bunge H.-P., 2009. Full seismic waveform tomography for upper-mantle structure in the Australasian region using adjoint methods, Geophys. J. Int. , 179( 3), 1703– 1725. https://doi.org/10.1111/j.1365-246X.2009.04368.x Google Scholar CrossRef Search ADS   Fichtner A., Bleibinhaus F., Capdeville Y., 2011. Full Seismic Waveform Modelling and Inversion , Springer. Google Scholar CrossRef Search ADS   Furumura T., Kennett B., Furumura M., 1998. Seismic wavefield calculation for laterally heterogeneous whole earth models using the pseudospectral method, Geophys. J. Int. , 135( 3), 845– 860. https://doi.org/10.1046/j.1365-246X.1998.00682.x Google Scholar CrossRef Search ADS   Geodynamics.org, 2009. Computational infrastructure for geodynamics, http://www.geodynamics.org, last accessed 1 December 2017. Gottlieb S., Shu C.-W., 1998. Total variation diminishing Runge-Kutta schemes, Math. Comput. Am. Math. Soc. , 67( 221), 73– 85. https://doi.org/10.1090/S0025-5718-98-00913-2 Google Scholar CrossRef Search ADS   Heinecke A. et al.  , 2014. Petascale High Order Dynamic Rupture Earthquake Simulations on Heterogeneous Supercomputers, in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis SC14 , pp. 3– 14, IEEE, New Orleans, LA, USA, Gordon Bell Finalist. Hesthaven J.S., Warburton T., 2008. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications , vol. 54, Springer. Igel H., 2016. Computational Seismology—A Practical Introduction , Oxford Univ. Press. Google Scholar CrossRef Search ADS   Jameson A., Schmidt W., Turkel E., 1981. Numerical solution of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes, AIAA 14th Fluid and Plasma Dynamic Conference , June 23–25, 1981, Palo Alto, California, paper 1981-1259. Karypis G., Kumar V., 1999. A fast and highly quality multilevel scheme for partitioning irregular graphs, SIAM J. Sci. Comput. , 20, 359– 392. https://doi.org/10.1137/S1064827595287997 Google Scholar CrossRef Search ADS   Käser M., Dumbser M., 2006. An arbitrary high order discontinuous galerkin method for elastic waves on unstructured meshes. I: The two-dimensional isotropic case with external source terms, Geophys. J. Int. , 166( 2), 855– 877. https://doi.org/10.1111/j.1365-246X.2006.03051.x Google Scholar CrossRef Search ADS   Käser M., Dumbser M., de La Puente J., Igel H., 2007. An arbitrary high-order Discontinuous Galerkin method for elastic waves on unstructured meshes—III. Viscoelastic attenuation, Geophys. J. Int. , 168( 1), 224– 242. https://doi.org/10.1111/j.1365-246X.2006.03193.x Google Scholar CrossRef Search ADS   Komatitsch D., Tromp J., 1999. Introduction to the spectral-element method for 3-D seismic wave propagation, Geophys. J. Int. , 139( 3), 806– 822. https://doi.org/10.1046/j.1365-246x.1999.00967.x Google Scholar CrossRef Search ADS   Komatitsch D., Vilotte J.P., 1998. The spectral-element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures, Bull. seism. Soc. Am. , 88( 2), 368– 392. Komatitsch D., Vilotte J.P., Vai R., Castillo-Covarrubias J.M., Sánchez-Sesma F.J., 1999. The spectral-element method for elastic wave equations: application to 2D and 3D seismic problems, Int. J. Numer. Methods Eng. , 45( 9), 1139– 1164. https://doi.org/10.1002/(SICI)1097-0207(19990730)45:9<1139::AID-NME617>3.0.CO;2-T Google Scholar CrossRef Search ADS   Komatitsch D., Tsuboi S., Ji C., Tromp J., 2003. A 14.6 billion degrees of freedom, 5 teraflops, 2.5 terabyte earthquake simulation on the Earth Simulator, Proceedings of the ACM/IEEE Supercomputing SC’2003 conference , pp. 4– 11, Gordon Bell Prize winner article. LeVeque R.J., 2002. Finite Volume Methods for Hyperbolic Problems , vol. 31, Cambridge Univ. Press. Google Scholar CrossRef Search ADS   Li Y., 2011. Development of numerical simulation method for nonlinear elastodynamic: application to acoustic imaging of defect with the help of cavity chaotic transducer, PhD thesis , Ecole Centrale de Lille. Marfurt K., 1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave-equations, Geophysics , 49( 5), 533– 549. https://doi.org/10.1190/1.1441689 Google Scholar CrossRef Search ADS   Maufroy E. et al.  , 2015. Earthquake ground motion in the Mygdonian Basin, Greece: the E2VP verification and validation of 3D numerical simulation up to 4 Hz, Bull. seism. Soc. Am. , 105, 1398– 1418. https://doi.org/10.1785/0120140228 Google Scholar CrossRef Search ADS   Mercerat E.D., Glinsky N., 2015. A nodal high-order discontinuous Galerkin method for elastic wave propagation in arbitrary heterogeneous media, Bull. seism. Soc. Am. , 201, 1101– 1118. Moczo P., Robertsson J.O., Eisner L., 2007. The finite-difference time-domain method for modeling of seismic wave propagation, Adv. Geophys. , 48, 421– 516. Google Scholar CrossRef Search ADS   Moczo P., Kristek J., Galis M., 2014. The Finite-Difference Modelling of Earthquake Motions: Waves and Ruptures , Cambridge Univ. Press. Google Scholar CrossRef Search ADS   Patera A.T., 1984. A spectral selement method for fluid dynamics: laminar flow in a channel expansion, J. Comput. Phys. , 54, 468– 488. https://doi.org/10.1016/0021-9991(84)90128-1 Google Scholar CrossRef Search ADS   Pelties C., de la Puente J., Ampuero J.-P., Brietzke G.B., Käser M., 2012. Three-dimensional dynamic rupture simulation with a high-order discontinuous Galerkin method on unstructured tetrahedral meshes, J. geophys. Res. , 117, B02309, doi:10.1029/2011JB008857. https://doi.org/10.1029/2011JB008857 Peter D. et al.  , 2011. Forward and adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes, Geophys. J. Int. , 186( 2), 721– 739. https://doi.org/10.1111/j.1365-246X.2011.05044.x Google Scholar CrossRef Search ADS   Reed W.H., Hill T.R., 1973. Triangular mesh methods for the neutron transport equation, Los Alamos Report LA-UR-73-479. Robertsson J. O.A., Blanch J.O., Symes W.W., 1994. Viscoelastic finite-difference modeling, Geophysics , 59( 9), 1444– 1456. https://doi.org/10.1190/1.1443701 Google Scholar CrossRef Search ADS   Sánchez-Sesma F., Ramos-Martinez J., Campillo M., 1993. An indirect boundary element method applied to simulate the seismic response of alluvial valleys for incident P, S and Rayleigh waves, Earthq. Eng. Struct. Dyn. , 22( 4), 279– 295. https://doi.org/10.1002/eqe.4290220402 Google Scholar CrossRef Search ADS   Schroeder W., Martin K., Lorensen B., 2006. The Visualization Toolkit , Kitware. Schumacher F., Friederich W., Lamara S., 2016. A flexible, extendable, modular and computationally efficient approach to scattering-integral-based seismic full waveform inversion, Geophys. J. Int. , 204( 2), 1100– 1119. https://doi.org/10.1093/gji/ggv505 Google Scholar CrossRef Search ADS   Seriani G., Priolo E., 1994. Spectral element method for acoustic wave simulation in heterogeneous media, Finite Elem. Anal. Des. , 16(3–4), 337– 348. https://doi.org/10.1016/0168-874X(94)90076-0 Google Scholar CrossRef Search ADS   Tautges T.J., 2001. The generation of hexahedral meshes for assembly geometry: survey and progress, Int. J. Numer. Methods Eng. , 50( 12), 2617– 2642. https://doi.org/10.1002/nme.139 Google Scholar CrossRef Search ADS   Tessmer E., Kosloff D., 1994. 3-D elastic modeling with surface-topography by a Chebyshev spectral method, Geophysics , 59( 3), 464– 473. https://doi.org/10.1190/1.1443608 Google Scholar CrossRef Search ADS   Virieux J., 1986. P–SV wave propagation in heterogeneous media: velocity-stress finite-difference method, Geophysics , 51( 4), 889– 901. https://doi.org/10.1190/1.1442147 Google Scholar CrossRef Search ADS   Warburton T., 2006. An explicit construction of interpolation nodes on the simplex, J. Eng. Math. , 56( 3), 247– 262. https://doi.org/10.1007/s10665-006-9086-6 Google Scholar CrossRef Search ADS   Wilcox L.C., Stadler G., Burstedde C., Ghattas O., 2010. A high-order discontinuous Galerkin method for wave propagation through coupled elastic–acoustic media, J. Comput. Phys. , 229( 24), 9373– 9396. https://doi.org/10.1016/j.jcp.2010.09.008 Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# A nodal discontinuous Galerkin approach to 3-D viscoelastic wave propagation in complex geological media

, Volume 212 (3) – Mar 1, 2018
18 pages

Loading next page...

/lp/ou_press/a-nodal-discontinuous-galerkin-approach-to-3-d-viscoelastic-wave-hwOV7nZyUp
Publisher
The Royal Astronomical Society
Copyright
© The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx494
Publisher site
See Article on Publisher Site

### Abstract

Summary A nodal discontinuous Galerkin (NDG) approach is developed and implemented for the computation of viscoelastic wavefields in complex geological media. The NDG approach combines unstructured tetrahedral meshes with an element-wise, high-order spatial interpolation of the wavefield based on Lagrange polynomials. Numerical fluxes are computed from an exact solution of the heterogeneous Riemann problem. Our implementation offers capabilities for modelling viscoelastic wave propagation in 1-D, 2-D and 3-D settings of very different spatial scale with little logistical overhead. It allows the import of external tetrahedral meshes provided by independent meshing software and can be run in a parallel computing environment. Computation of adjoint wavefields and an interface for the computation of waveform sensitivity kernels are offered. The method is validated in 2-D and 3-D by comparison to analytical solutions and results from a spectral element method. The capabilities of the NDG method are demonstrated through a 3-D example case taken from tunnel seismics which considers high-frequency elastic wave propagation around a curved underground tunnel cutting through inclined and faulted sedimentary strata. The NDG method was coded into the open-source software package NEXD and is available from GitHub. Numerical modelling, Computational seismology, Wave propagation 1 INTRODUCTION In geophysics, numerical simulations are a key tool for understanding physical phenomena taking place on and inside the earth. As they usually make predictions about observable quantities, they also play an essential role in inferring the current state and the physical properties of the earth’s interior. A phenomenon of particular importance in this context is seismic waves which travel through the entire globe and allow a unique look into the deep earth. Seismic waves are also useful on smaller scales, for example to find and characterize oil and gas reservoirs or to identify unknown objects in front of a newly excavated tunnel. Even on laboratory scale, elastic waves can help to increase knowledge about material behaviour, for example, to understand the mechanisms of crack propagation in rocks. In computational seismology and seismics, the need for high performance numerical simulations is continuously growing. Generally, in order to conduct seismic simulations, several different methods can be used. Overviews are given in the textbooks by Moczo et al. (2014) and Igel (2016). The probably most common method is the finite difference (FD) method (Virieux 1986; Moczo et al. 2007; Bohlen 2002) in which spatial derivatives in the governing differential equations are approximated by finite differences on an equally spaced grid (Robertsson et al. 1994). The strength of FD method is the simplicity of its numerical representation, which is however accompanied by numerical dispersion, especially for surface waves, and inaccurate representation of internal and external boundaries. Grid adaption, taking into account varying model properties, is hardly possible and a compromise between accuracy and computational cost must be carefully worked out for every simulation. Spectral and pseudospectral methods are used to simulate seismic wave propagation with high order accuracy (Tessmer & Kosloff 1994; Furumura et al. 1998). These methods rely on global basis functions such as Chebychev- or Legendre polynomials leading to simulations with little numerical dispersion. Only few points per wavelength are needed to achieve sufficient accuracy. However, the choice of global basis functions restricts the spectral methods to very smooth models, so that strong material contrasts and discontinuities are not represented adequately. Boundary Integral Methods (Sánchez-Sesma et al. 1993; Chen 2007) have been successfully used to simulate seismic wave propagation in domains with layers and interfaces, but are restricted to a limited amount of layers and are not well suited for 3-D problems as the involved matrices become very large and ill-conditioned. The finite element method (FEM) was used by Bao et al. (1998) and Marfurt (1984) to simulate seismic wave propagation. The FEM is based on the idea to subdivide the computational domain into elements, on which the equations can be solved. Albeit handling complex geometries very well, the solution of the wave equation is represented by low order polynomials resulting in numerical dispersion. For an accurate simulation, very large matrix systems have to be solved iteratively leading to high numerical costs with poor parallelization properties. Most commercial simulation software is based on FEM and is thus not well suited for seismic wave propagation problems. The idea of combining the spectral approximation properties with the flexibility of the FEM led to the development of the spectral element method (SEM; Patera 1984; Seriani & Priolo 1994) which was adapted to seismic wave propagation by Komatitsch & Vilotte (1998). The SEM uses a high order approximation of the wavefield inside the elements based on Lagrange polynomials. For the SEM, a variational formulation of the governing equations is used and a diagonal mass matrix can be constructed when taking Gauss–Lobatto–Legendre quadrature points within Cartesian elements (quadrilateral and hexahedral elements) for numerical integration. This leads to a fully explicit scheme making the SEM computationally very efficient and allowing for a simple parallelization via domain decomposition. The SEM is now widely used in the seismic community, especially for global and regional earthquake simulations (Komatitsch et al. 1999, 2003; Peter et al. 2011). Several free software implementations of the SEM exist for this purpose such as SPECFEM, SES3D and RegSEM (Geodynamics.org 2009; Fichtner et al. 2009; Cupillard et al. 2012). Unfortunately, the generation of well-structured hexahedral meshes can be very cumbersome for complex geological environments (Tautges 2001). The demand for a high order approximation of the seismic wavefield paired with the geometrical flexibility of an unstructured tetrahedral mesh led to the adaption of the discontinuous Galerkin (DG) method to seismic problems. The DG method was first proposed by Reed & Hill (1973) for solving the linear neutron transport equation. An extension to several research fields such as fluid dynamics (Bassi & Rebay 1997) and electrodynamic simulations (Cockburn et al. 2004) shows the capability of the DG method for simulating physical problems. An extensive mathematical research on the method was done by Cockburn et al. (2000) and Hesthaven & Warburton (2008). The DG method can be seen as a combination of the finite volume method (FVM; LeVeque 2002) and the SEM to solve partial differential equations with high accuracy. In contrast to the SEM and FEM, the variational formulation of the governing equations is done locally for every element and not globally. The solution obtained with the DG is allowed to be discontinuous across element boundaries. As in finite volume methods, this local formulation of the DG leads to the concept of numerical fluxes for exchanging information between adjacent elements. This results in an explicit semi-discrete numerical scheme which is very well suited for parallelization with domain decomposition techniques as only fluxes have to be communicated at element boundaries. In the DG method, one can choose a modal approach where the solution is represented by the expansion coefficients of polynomial basis functions, or a nodal approach where the solution is represented by its values at the anchor points of Lagrange polynomials used to approximate the solution within each element. For seismic wave propagation, the modal DG method was first introduced by Käser & Dumbser (2006) together with the concept of Riemann fluxes and arbitrary high order derivatives (ADER-DG). Here, time integration is coupled with the spatial resolution resulting in the same high order approximation of the solution in space and time. This approach results in a rather complicated numerical scheme due to the use of the Cauchy–Kovalewski procedure for the time integration. Intensive research on extending the ADER-DG to wave propagation problems in viscoelastic (Käser et al. 2007) and anisotropic media (de la Puente 2008) was conducted. It was also adapted to the study of rupture processes (Pelties et al. 2012). From this research sprang the community code SeisSol which became open-source in 2015 (https://github.com/SeisSol/SeisSol). The following work on SeisSol focused on increasing the numerical efficiency and adapting it to extremely large-scale seismic simulations and dynamic earthquake rupturing problems (Heinecke et al. 2014). Nodal high-order DG methods for seismic wave propagation were first developed by Delcourte et al. (2009) and Etienne et al. (2010) for the isotropic, elastic case. They use a pseudo-conservative formulation of the velocity–stress equations where the elastic material properties are assembled in a single diagonal matrix. Numerical fluxes are calculated by averaging velocity–stress values of adjacent elements. The solution is advanced in time using a leap-frog time integration. Delcourte et al. (2009) treated the 2-D case whereas Etienne et al. (2010) extended the approach to a low-order nodal DG (NDG) method for 3-D elastic wave propagation that exploits the hp-adaptivity of the DG method. The latter approach was benchmarked against other codes implementing FD, spectral element and pseudo-spectral methods in an earthquake ground motion application (Maufroy et al. 2015). Besides that, Wilcox et al. (2010) studied wave propagation for elastic-acoustic media within the DG framework based on hexahedral elements in the context of global earthquake simulation and Li (2011) used the NDG method to conduct 2-D simulations for non-linear elastic wave phenomena. One peculiarity of most DG methods (including SeisSol) is that material properties have to be constant within the elements. In case of continuously varying material properties, this assumption may lead to inaccuracies if the dominating wavelength approaches the element size because then the continuous structure is no longer properly approximated by a piecewise-constant representation. Furthermore, in extremely complicated media where discontinuities of material properties may cross an element, larger errors of the classical DG approach are to be expected. This problem was solved by Mercerat & Glinsky (2015) who extended the NDG approaches of Delcourte et al. (2009) and Etienne et al. (2010) to the case of heterogeneous elements for the 2-D, purely elastic case and demonstrated that, given equal element size, their approach yields more accurate results than the constant-element approaches. Unfortunately, it is unclear whether the pseudo-conservative form of the velocity–stress equations of Etienne et al. (2010) can be extended to the viscoelastic and anisotropic case. In this paper, we present a high-order nodal approach to the DG method based on a general mathematical treatment by Hesthaven & Warburton (2008) to model seismic wave propagation in 2-D and 3-D viscoelastic media with complex geological structures. Discretization is done using triangular (2-D) or tetrahedral (3-D) meshes. Instead of the pseudo-conservative form of Etienne et al. (2010), we use a classical, first-order in time velocity–stress formulation which allows the treatment of viscoelastic and also anisotropic media. We provide a justification for the necessity of introducing numerical fluxes into the computational scheme. We give explicit expressions for these fluxes, derived from an exact solution of the Riemann problem for the viscoelastic wave equation and expressed in terms of waves propagating away from the discontinuity. In contrast to the approach by Käser & Dumbser (2006) implemented in SeisSol, our fluxes honour differing material properties in adjacent elements. This way of treating the numerical fluxes allows for an easy extension of the method to cases where certain, potentially time-dependent discontinuities of particle velocities or stresses are prescribed at element boundaries, for example in the presence of fractures or during dynamic rupturing. The NDG approach leads to a simple, efficient and extendable numerical code with little logistical overhead. It is very intuitive and convenient as the solution vector directly represents the values of particle velocities and stresses at the nodal points. Concerning the evaluation of surface integrals on element faces required for calculating numerical fluxes, NDG is more efficient than modal DG as only information of nodes on the element faces is required. Our NDG code is capable of solving a wide range of seismic problems, from small-scale field experiments to large-scale earth models. Extensions of the code to unsaturated poroelastic (Boxberg et al. 2017) and fractured media are being worked on. Our implementation of the method has an interface to the full waveform inversion code ASKI (Schumacher et al. 2016) and offers the opportunity to compute adjoint wavefields. The paper is structured as follows: We briefly introduce the theoretical background of our approach and derive a numerical scheme suited for efficient parallel implementation. Some detail is dedicated to the derivation of the heterogeneous numerical fluxes for elastic and viscoelastic media. The general structure of the implemented algorithm is lined out and our newly developed program package NEXD is presented. The code is validated in 2-D by comparing its solution to Lamb’s problem with a semi-analytical one and by carrying out a p-convergence test. In addition, we benchmark the 2-D code against an analytical solution for a viscoelastic full space. In 3-D, we compare to results of SPECFEM (Geodynamics.org 2009) for a two-layer half-space. The potential and functionality of our NDG approach is demonstrated by means of a complex example taken from tunnel seismics where we consider 3-D propagation of elastic waves around a curved tunnel that cuts through faulted sedimentary strata. 2 A NODAL DISCONTINUOUS GALERKIN SCHEME FOR ELASTIC WAVE PROPAGATION The following section describes the theoretical basis of the numerical formulations developed and implemented in this paper. The main part deals with the adaption of the NDG method (Hesthaven & Warburton 2008) to the simulation of 3-D elastic wave propagation on unstructured tetrahedral meshes. 2.1 Velocity–stress formulation of the elastic wave equation in 3-D Following Virieux (1986) and Dumbser & Käser (2006) the elastic wave equation can be written as a system of first order hyperbolic equations. For the 3-D, elastic and isotropic case the system reads:   \begin{eqnarray} \frac{\partial }{\partial t } \sigma _{xx} -(\lambda +2\mu )\frac{\partial }{\partial x} u -\lambda \frac{\partial }{\partial y} v -\lambda \frac{\partial }{\partial z} w &=& s_{1}, \nonumber \\ \frac{\partial }{\partial t } \sigma _{yy} -\lambda \frac{\partial }{\partial x} u -(\lambda +2\mu ) \frac{\partial }{\partial y} v -\lambda \frac{\partial }{\partial z} w &=& s_{2}, \nonumber \\ \frac{\partial }{\partial t } \sigma _{zz} -\lambda \frac{\partial }{\partial x} u -\lambda \frac{\partial }{\partial y} v -(\lambda +2\mu ) \frac{\partial }{\partial z} w &=& s_{3}, \nonumber \\ \frac{\partial }{\partial t} \sigma _{xy}-\mu \left(\frac{\partial }{\partial x} v +\frac{\partial }{\partial y} u\right) &=& s_{4}, \nonumber \\ \frac{\partial }{\partial t} \sigma _{yz}-\mu \left(\frac{\partial }{\partial z} v +\frac{\partial }{\partial y} w\right) &=& s_{5}, \nonumber \\ \frac{\partial }{\partial t} \sigma _{xz}-\mu \left(\frac{\partial }{\partial x} w +\frac{\partial }{\partial z} u\right) &=& s_{6}, \nonumber \\ \frac{\partial }{\partial t} u - \frac{1}{\rho }\frac{\partial }{\partial x} \sigma _{xx} - \frac{1}{\rho }\frac{\partial }{\partial y} \sigma _{xy} - \frac{1}{\rho }\frac{\partial }{\partial z} \sigma _{xz} &=& s_{7} \nonumber ,\\ \frac{\partial }{\partial t} v - \frac{1}{\rho }\frac{\partial }{\partial x} \sigma _{xy} - \frac{1}{\rho }\frac{\partial }{\partial y} \sigma _{yy} - \frac{1}{\rho }\frac{\partial }{\partial z} \sigma _{yz} &=& s_{8}, \nonumber \\ \frac{\partial }{\partial t} w - \frac{1}{\rho }\frac{\partial }{\partial x} \sigma _{xz} - \frac{1}{\rho }\frac{\partial }{\partial y} \sigma _{yz} - \frac{1}{\rho }\frac{\partial }{\partial z} \sigma _{zz} &=& s_{9}, \end{eqnarray} (1)where σlk represent Cartesian stress components and u, v and w denote the x-, y- and z-components of the particle velocities. We use here Lamé’s elastic moduli λ and μ. Density is denoted by ρ. Note that all quantities are functions of space and time. A single force source can be defined by s = (0, 0, 0, 0, 0, 0, ax, ay, az)T, where the components of the acceleration ax, ay and az act in the last three components of the source vector. It is also possible to realize a moment tensor source which acts in the six stress components of the source vector s. In matrix-vector notation, this system of linear hyperbolic differential equations can be written in a more compact form as   \begin{eqnarray} \frac{\partial \mathbf {u}}{\partial t} + \mathbf {A} \frac{\partial \mathbf {u}}{\partial x} + \mathbf {B} \frac{\partial \mathbf {u} }{\partial y} + \mathbf {C} \frac{\partial \mathbf {u} }{\partial z} = \mathbf {s}. \end{eqnarray} (2)The vector   \begin{eqnarray*} \mathbf {u} = (\sigma _{xx}, \sigma _{yy},\sigma _{zz}, \sigma _{xy},\sigma _{yz},\sigma _{xz}, u, v, w)^{T} \nonumber \end{eqnarray*} contains the Nu = 9 velocity–stress variables, which all depend on time t and spatial coordinates x. The compact form contains the Jacobi matrices A = Aij(x), B = Bij(x) and C = Cij(x) describing the elastic properties of the medium (Dumbser & Käser 2006). In 3-D, the indices i, j range from 1 to 9. 2.2 Element-based integral forms of the velocity–stress equations The NDG method is built on an integral formulation of the velocity–stress equations. A computational domain, Ω, is defined encompassing the region of interest. It is subdivided into NE non-overlapping tetrahedral elements De, e = 1, …, NE  \begin{eqnarray} \Omega = \bigcup _{e=1}^{N_E} D^{e}. \end{eqnarray} (3)Within each element De the velocity–stress vector is represented by an interpolation formula of the form:   \begin{eqnarray} \mathbf {u}(\mathbf {x},t) = \sum _{j=1}^{N_p} \mathbf {u}(\mathbf {x}_j,t)\, l_j(\mathbf {x}), \quad \mathbf {x} \in D^e, \end{eqnarray} (4)where lj(x) is a multidimensional Lagrange polynominal anchored at Np interpolation points xi located within element De and on its faces. Since these polynomials satisfy the relation li(xj) = δij, the expansion coefficients are indeed identical to the values of the velocity–stress vector at the anchor points. Note also that no a priori assumptions are made on the continuity of the velocity–stress vector across adjacent elements. This fact is one of the name-giving properties of the DG method. Anchor points on adjacent element faces always occur as double nodes principally allowing the definition of a discontinuous velocity–stress vector across these faces. The integral form is now constructed by forming the scalar product of eq. (2) with a set of Nu × Np test vectors and integrating over element e. A system of Nu × Np equations results that can be solved for the nodal values of the velocity–stress vector. Here, we choose test vectors of the form   \begin{eqnarray} \mathbf {w}_{ni} = \mathbf {e}_n\, l_i(\mathbf {x}), \quad 1 \le i \le N_p,\quad 1 \le n \le N_u, \end{eqnarray} (5)where the nth component of en is 1 and all other components are 0. This is equivalent to multiplying the velocity–stress equations by the matrix Ili(x) where I is the identity matrix of the Nu-dimensional velocity–stress space. For each element e, we obtain Np vector equations:   \begin{eqnarray} \int _{D^e} \frac{\partial \mathbf {u}}{\partial t}l_i\, \text{d}V &=& -\int _{D^e}\left(\mathbf {A}\frac{\partial \mathbf {u}}{\partial x} +\mathbf {B}\frac{\partial \mathbf {u}}{\partial y} +\mathbf {C}\frac{\partial \mathbf {u}}{\partial z}\right)l_i \text{d}V \nonumber \\ &&+ \int _{D^e} \mathbf {s} \ l_i \, \text{d}V, \quad 1 \le i \le N_p . \end{eqnarray} (6)An integration by parts leads to:   \begin{eqnarray} \int _{D^e} \frac{\partial \mathbf {u}}{\partial t} \, l_i \text{d}V &=& -\int _{\partial D^e}(\mathbf {f}_xn_x+\mathbf {f}_yn_y+\mathbf {f}_zn_z)\, l_i\, \text{d}\Sigma \nonumber \\ &&+\int _{D^e} \left(\frac{\partial (\mathbf {A}l_i)}{\partial x} +\frac{\partial (\mathbf {B}l_i)}{\partial y} +\frac{\partial (\mathbf {C}l_i)}{\partial z}\right)\mathbf {u}\, \text{d}V \nonumber \\ &&+ \int _{D^e} \mathbf {s}\, l_i\, \text{d}V, \quad 1 \le i \le N_p . \end{eqnarray} (7)Here, we have introduced the components of the unit normal vector on the element faces, nx, ny, nz, the surface element dΣ and the fluxes fx, fy and fz defined as   \begin{eqnarray} \mathbf {f}_x = \mathbf {A}\mathbf {u}+\mathbf {r}_x,\quad \mathbf {f}_y = \mathbf {B}\mathbf {u}+\mathbf {r}_y,\quad \mathbf {f}_z = \mathbf {C}\mathbf {u}+\mathbf {r}_z . \end{eqnarray} (8)While surface integrals over expressions such as Au are the mathematically correct terms appearing after integration by parts, we have added here additional flux terms, rx, ry and rz, whose necessity is justified and discussed later. A second useful integral form of the velocity–stress equations can be obtained after another integration by parts. Expressions of the type Au appear again under surface integrals but now cancel the ones appearing in eq. (8):   \begin{eqnarray} \int _{D^e} \frac{\partial \mathbf {u}}{\partial t} \, l_i \text{d}V &=& -\int _{\partial D^e}(\mathbf {r}_xn_x+\mathbf {r}_yn_y+\mathbf {r}_zn_z)\, l_i \text{d}\Sigma \nonumber \\ &&-\int _{D^e} \left(\mathbf {A}\frac{\partial \mathbf {u}}{\partial x} +\mathbf {B}\frac{\partial \mathbf {u}}{\partial y} +\mathbf {C}\frac{\partial \mathbf {u}}{\partial z}\right)\, l_i\, \text{d}V \nonumber \\ &&+ \int _{D^e} \mathbf {s}\, l_i\, \text{d}V, \quad 1 \le i \le N_p . \end{eqnarray} (9)For both integral forms of the velocity–stress equations, we require the Np vector equations (7) or (9) to be fulfilled for each element. 2.3 Two semi-discrete schemes for the NDG method To simplify the notation in the following, we will use the abbreviations   \begin{eqnarray} \mathbf {u}^j(t) = \mathbf {u}(\mathbf {x}_j,t) \quad \mbox{and}\quad \mathbf {s}^j(t) = \mathbf {s}(\mathbf {x}_j,t),\quad 1 \le j \le N_p \end{eqnarray} (10)for the value of the velocity–stress vector and source vector, respectively, at the jth anchor point. Moreover, we define normal fluxes defined by   \begin{eqnarray} \mathbf {f}_n &=& \mathbf {f}_xn_x+\mathbf {f}_yn_y+\mathbf {f}_zn_z\quad \mbox{and} \nonumber \\ \mathbf {r}_n &=& \mathbf {r}_xn_x+\mathbf {r}_yn_y+\mathbf {r}_zn_z . \end{eqnarray} (11)By substituting eq. (4) into eq. (7) and by further assuming that the Jacobi matrices are constant within each element, we obtain   \begin{eqnarray} &&{\sum _{j=1}^{N_p} \frac{\partial \mathbf {u}^j}{\partial t} \int _{D^e} l_j(\mathbf {x}) l_i(\mathbf {x}) \text{d}V = \mathbf {A} \sum _{j=1}^{N_p} \mathbf {u}^j \int _{D^e} l_j(\mathbf {x}) \frac{\partial l_i(\mathbf {x})}{\partial x}\text{d}V }\nonumber \\ &&+\mathbf {B} \sum _{j=1}^{N_p} \mathbf {u}^j \int _{D^e} l_j(\mathbf {x}) \frac{\partial l_i(\mathbf {x})}{\partial y}\text{d}V + \mathbf {C} \sum _{j=1}^{N_p} \mathbf {u}^j \int _{D^e} l_j(\mathbf {x}) \frac{\partial l_i(\mathbf {x})}{\partial z}\text{d}V \nonumber \\ &&- \int _{\partial D^e} \mathbf {f}_nl_i(\mathbf {x})\text{d}\Sigma + \sum _{j=1}^{N_p}\mathbf {s}^j \int _{D^e} l_j(\mathbf {x}) l_i(\mathbf {x}) \text{d}V,\quad 1 \le i \le N_p . \end{eqnarray} (12)Eq. (12) motivates the introduction of the mass matrix   \begin{eqnarray} M^e_{ij} = \int _{D^e} l_i(\mathbf {x}) l_j(\mathbf {x}) \text{d}V, \quad 1 \le i,j \le N_p \end{eqnarray} (13)and the stiffness matrices   \begin{eqnarray} S^e_{ij,x} &=& \int _{D^e} l_i(\mathbf {x}) \frac{\partial l_j(\mathbf {x})}{\partial x} \text{d}V, \quad 1 \le i,j \le N_p \nonumber \\ S^e_{ij,y} &=& \int _{D^e} l_i(\mathbf {x}) \frac{\partial l_j(\mathbf {x})}{\partial y} \text{d}V, \quad 1 \le i,j \le N_p \nonumber \\ S^e_{ij,z} &=& \int _{D^e} l_i(\mathbf {x}) \frac{\partial l_j(\mathbf {x})}{\partial z} \text{d}V, \quad 1 \le i,j \le N_p , \end{eqnarray} (14)by which we can rewrite eq. (12), using Einstein summation convention, as   \begin{eqnarray} M^e_{ij}\frac{\partial \mathbf {u}^j}{\partial t} &=& \mathbf {A} S^e_{ji,x}\mathbf {u}^j + \mathbf {B} S^e_{ji,y}\mathbf {u}^j+ \mathbf {C} S^e_{ji,z}\mathbf {u}^j + M^e_{ij}\mathbf {s}^j \nonumber \\ &&- \int _{\partial D^e}\mathbf {f}_n\, l_i(\mathbf {x})\, \text{d}\Sigma , \quad 1 \le i,j \le N_p. \end{eqnarray} (15) Similarly, by inserting the nodal representation (4) into the second integral form of the velocity–stress equations (9), we obtain   \begin{eqnarray} M^e_{ij}\frac{\partial \mathbf {u}^j}{\partial t} &=& -\mathbf {A} S^e_{ij,x}\mathbf {u}^j - \mathbf {B} S^e_{ij,y}\mathbf {u}^j- \mathbf {C} S^e_{ij,z}\mathbf {u}^j + M^e_{ij}\mathbf {s}^j \nonumber \\ &&- \int _{\partial D^e}\mathbf {r}_n\, l_i(\mathbf {x})\, \text{d}\Sigma , \quad 1 \le i,j \le N_p. \end{eqnarray} (16)Subtle but notable differences between the two semi-discrete schemes, eqs (15) and (16), are different signs in front of the Jacobi matrices, the different order of indices in the stiffness matrices and the occurrence of different fluxes. Given a source vector and initial values for velocity–stress, we can use the above system of equations for each element to calculate an update of the velocity–stress vector provided we know how to compute the additional normal flux rn. 2.4 Normal Riemann fluxes The additional normal flux needs to be introduced because, by construction of the DG method, the velocity–stress vector may be discontinuous on the faces of adjacent elements. These discontinuities act as sources of seismic waves that propagate away from the discontinuity in both directions. If the flux was omitted or evaluated from the values of the velocity–stress vector of one element alone, the discontinuity would go unnotified and the additional waves would be missing in the solution. Consider a situation of two adjacent element faces with normal pointing into the x-direction. The source term is for now assumed to vanish within these two elements. Furthermore, we assume that the velocity–stress vector varies smoothly inside the elements but is discontinuous across the contact surface. We can always decompose such a velocity–stress distribution into two parts: one that is continuous at the contact surface and one that is discontinuous there and independent of x in both elements. Interpreting eq. (9) as a recipe for updating the velocity–stress vector with time, we recognize that the volume integral on the right-hand side will generate temporal changes of velocity–stress for the continuous part. For the discontinuous part, however, there will be no contribution from the term A∂u/∂x which is the only one that can produce waves propagating perpendicular to the contact surface. Hence, omitting the additional normal flux rn, implies loosing the waves generated by the discontinuity. To obtain specific expressions for the additional normal flux, we consider the so-called normal Riemann problem (LeVeque 2002). There we seek a solution of the velocity–stress equations for purely 1-D elastic wave propagation normal to an interface situated at x = 0 at both sides of which the velocity–stress vector is initially constant but discontinuous across the interface. Material properties on each side of the interface are allowed to be different. Source-free elastic wave propagation on both sides is described by   \begin{eqnarray} &&{\frac{\partial \mathbf {u}^{(1,2)}}{\partial t} + \mathbf {A}^{(1,2)} \frac{\partial \mathbf {u}^{(1,2)}}{\partial x} = \mathbf {0} \quad \mbox{with}} \nonumber \\ &&\mathbf {u}^{(1)}(x,t=0) = \mathbf {U}^{(1)}\quad \mbox{and}\quad \mathbf {u}^{(2)}(x,t=0) = \mathbf {U}^{(2)} , \end{eqnarray} (17)where the superscripts 1 and 2 refer to the left (x < 0) and right side (x > 0) of the interface, respectively. The solution to this problem is derived in detail by LeVeque (2002) in the context of the finite volume method. For positive times, we expect that waves will be generated by the discontinuity that propagate away from the interface. Waves with negative velocity will propagate to the left and waves with positive velocity to the right. In 3-D, matrix A has 9 eigenvalues, representing wave speeds c1−9, which we order in a way that c1−3 are negative, c4−6 are zero and c7−9 are positive. To calculate c1−3, we use matrix A(1), to calculate c7−9, we use matrix A(2). The amplitudes of the associated waves are obtained by decomposing the negative velocity–stress jump in terms of the corresponding right eigenvectors of matrices A(1, 2), denoted by $$\mathbf {v}^{(1,2)}_k$$, taking the ones for medium 1 when the corresponding wave speed is negative and the ones for medium 2 when the corresponding wave speed is positive:   \begin{eqnarray} \mathbf {U}^{(1)}-\mathbf {U}^{(2)} = \sum _{k=1}^3\gamma _k\mathbf {v}^{(1)}_k +\sum _{k=4}^6\gamma _k\mathbf {v}_k+\sum _{k=7}^9\gamma _k\mathbf {v}^{(2)}_k . \end{eqnarray} (18)For eigenvectors 4–6 we do not specify which side they belong to because they are independent of the material properties. The coefficients γk specify the desired amplitudes of the propagating waves. They can be calculated by multiplying the above equation by the appropriate reciprocal vectors of $$\mathbf {v}^{(1,2)}_k$$. Velocity–stress for positive times can then be written as   \begin{eqnarray} \mathbf {u}^{(1)}(x,t>0) &=& \mathbf {U}^{(1)}-\sum _{k=1}^3\gamma _k\mathbf {v}^{(1)}_kH\left(t+\frac{x}{|c_k|}\right) \nonumber \\ \mathbf {u}^{(2)}(x,t>0) &=& \mathbf {U}^{(2)}+\sum _{k=7}^9\gamma _k\mathbf {v}^{(2)}_kH\left(t-\frac{x}{c_k}\right) , \end{eqnarray} (19)where H(t) is a Heaviside function. With this solution, the discontinuity at the interface is reduced to   \begin{eqnarray} \mathbf {u}^{(2)}(0,t)-\mathbf {u}^{(1)}(0,t) &=& -\sum _{k=4}^6\gamma _k\mathbf {v}_k . \end{eqnarray} (20)Due to the structure of the eigenvectors, the remaining discontinuity occurs on the components σyy, σzz and σyz of the velocity–stress vector which do not generate waves propagating away from the interface. The remaining stress components as well as the particle velocities are continuous as required for a welded contact of two materials at the interface. Explicit expressions for the normal flux are obtained by going back to eq. (9) but staying with the purely 1-D Riemann problem. Instead of a tetrahedral element as integration volume, we choose a rectangular box attached to the right side of the interface and extending a distance L into the positive x-direction. After a small time step Δt, the change of velocity–stress in the box is given by the sum of right-propagating waves in eq. (19). As discussed above, this change must be produced by the additional normal flux rn. To avoid contributions from the right face of the box, we choose its width L greater than the distance the fastest wave can travel within time Δt. Then, instead of eq. (9), we get:   \begin{eqnarray} \int _\Sigma \text{d}\Sigma \int _0^L \Delta \mathbf {u}\, l_i(x) \text{d}x = -\Delta t\int _\Sigma \mathbf {r}_n\, l_i(0) \text{d}\Sigma , \end{eqnarray} (21)where Σ denotes the interface with outward normal pointing into the negative x-direction and li(x) is considered as a 1-D Lagrange polynomial of order N. Inserting the right-propagating waves for Δu, we find   \begin{eqnarray} &&{\int _\Sigma \text{d}\Sigma \int _0^L \sum _{k=7}^9\gamma _k\mathbf {v}^{(2)}_kH\left(\Delta t-\frac{x}{c_k}\right)l_i(x) \text{d}x }\nonumber \\ &&{\quad =\int _\Sigma \text{d}\Sigma \sum _{k=7}^9\gamma _k\mathbf {v}^{(2)}_k\int _0^{c_k\Delta t} l_i(x) \text{d}x} \nonumber \\ &&{\quad \approx \int _\Sigma \text{d}\Sigma \sum _{k=7}^9\gamma _k\mathbf {v}^{(2)}_kc_k\Delta t\, l_i(0) = -\Delta t\int _\Sigma \mathbf {r}_nl_i(0)\text{d}\Sigma,} \end{eqnarray} (22)from which we conclude for the normal flux on the right-hand side of the interface with outward normal pointing to the left:   \begin{eqnarray} \mathbf {r}^2_n = -\sum _{k=7}^9 c_k\gamma _k\mathbf {v}^{(2)}_k . \end{eqnarray} (23)The normal flux on the left-hand side of the interface with outward normal pointing to the right is obtained by analogously considering a box on the left-hand side and inserting the left-propagating waves for the change of the velocity–stress vector. Then, we find   \begin{eqnarray} \mathbf {r}^1_n = +\sum _{k=1}^3 |c_k|\gamma _k\mathbf {v}^{(1)}_k . \end{eqnarray} (24)Note that the evaluation of the normal flux is easily extended to situations where a certain, eventually time-dependent, discontinuity of the velocity–stress vector is prescribed at the interface, for example during dynamic rupturing or in the presence of fractures. To compute the normal flux in practice for a given element face, the velocity stress vector is rotated into a coordinate frame whose x-axis points into the direction of the outward normal of the face. The flux rx is computed in this frame from the velocity–stress jump and then rotated back into the original coordinate frame. This procedure is greatly alleviated by the fact that for isotropic elastic media the Jacobi matrices are invariant under rotation. It is also noted here that the expressions for the fluxes in eqs (23) and (24) account for differing material properties in adjacent elements. The normal Riemann flux given by Dumbser & Käser (2006) is calculated using only material properties of the element under consideration. The concept of fluxes in the DG method may also be exploited to realize particular boundary conditions. For example, an absorbing boundary condition can be achieved by simply cancelling the incoming fluxes at the boundaries. Similarly, by mirroring the stresses via the fluxes at the boundaries, a free boundary condition can be constructed. 2.5 NDG on tetrahedral elements One key feature of the NDG is its local character. It is not necessary to invert global matrices thus gaining flexibility for numerical implementation, especially for the parallelization of the method. For the high order approximation of the method, it is necessary to find a set of anchor points xi, i = 1, …, Np for the multidimensional Lagrange polynomials. For this purpose, the elements are mapped to a standard reference element, shown in Fig. 1 for tetrahedra, by a coordinate transformation of the form $$\mathbf {x}({\boldsymbol\xi })$$ where $${\boldsymbol\xi }=(\xi ,\eta ,\zeta )$$ denotes the position vector in the reference frame. The anchor points within a tetrahedral element can be constructed by a technique called ‘warp and blending’ introduced by Warburton (2006) (Fig. 1). The number of points is dependent on the order N of the interpolation and is, for the 3-D case on a tetrahedral element (Hesthaven & Warburton 2008),   \begin{eqnarray} N_p = \frac{(N+1)(N+2)(N+3)}{6}. \end{eqnarray} (25)Using the standard reference element, the mass matrix may be calculated once in advance according to   \begin{eqnarray} M_{ij} = \int _{I} l_i({\boldsymbol\xi }) l_j({\boldsymbol\xi }) \text{d}V_{\mbox{ref}} = \frac{1}{J^e}\int _{D_k} l_i(\mathbf {x}) l_j(\mathbf {x}) \text{d}V\, \end{eqnarray} (26)for the physical element De and reference element I, where Je is the constant Jacobian of the element De describing the volume change during the transformation from the physical to the reference element. The stiffness matrices Sij,ξ, Sij,η and Sij,ζ in the reference frame are given by   \begin{eqnarray} S_{ij,\xi _l} = \int _\mathbf {I} l_i({\boldsymbol\xi }) \frac{\partial l_j({\boldsymbol\xi })}{\partial \xi _l} \text{d}V_{\mbox{ref}}, \end{eqnarray} (27)and are related to the ones defined in physical space by (using Einstein summation convention)   \begin{eqnarray} S^e_{ij,x_m} = J^e\frac{\partial \xi _l}{\partial x_m}\, S_{ij,\xi _l}. \end{eqnarray} (28)Evaluation of the stiffness matrices in the reference frame is conveniently done by expanding the derivatives of the Lagrange polynomials in terms of Lagrange polynomials:   \begin{eqnarray} \frac{\partial l_j({\boldsymbol\xi })}{\partial \xi _l} = \sum _{r=1}^{N_p}\frac{\partial l_j}{\partial \xi _l}({\boldsymbol\xi }_r)l_r({\boldsymbol\xi }), \end{eqnarray} (29)leading to   \begin{eqnarray} S_{ij,\xi _l} = \sum _{r=1}^{N_p}\frac{\partial l_j}{\partial \xi _l}({\boldsymbol\xi }_r) \int _\mathbf {I} l_i({\boldsymbol\xi }) l_r({\boldsymbol\xi })\text{d}V_{{\rm ref}} . \end{eqnarray} (30)Defining the derivative matrix   \begin{eqnarray} D_{rj,\xi _l} = \frac{\partial l_j}{\partial \xi _l}({\boldsymbol\xi }_r) \end{eqnarray} (31)the stiffness matrices can be expressed as (using summation convention)   \begin{eqnarray} S_{ij,\xi _l} = M_{ir}D_{rj,\xi _l}\quad \mbox{and}\quad S^e_{ij,x_m} = J^e\frac{\partial \xi _l}{\partial x_m}M_{ir}D_{rj,\xi _l} . \end{eqnarray} (32)The mass matrix and the derivative matrix on the reference element can be conveniently evaluated by expanding the Lagrange polynomials and their derivatives in terms of orthogonal Jacobi polynomials (Hesthaven & Warburton 2008). Figure 1. View largeDownload slide Standard reference element for a tetrahedral element. Coordinates in this reference frame vary between −1 and +1. Small grey squares show position of anchor points for an interpolation order of N = 4. Figure 1. View largeDownload slide Standard reference element for a tetrahedral element. Coordinates in this reference frame vary between −1 and +1. Small grey squares show position of anchor points for an interpolation order of N = 4. For quantities defined on an element face f such as the fluxes, we perform an expansion in terms of 2-D Lagrange polynomials $$l_s^{2D}(\mathbf {x})$$ anchored at the face nodes. These are identical to those 3-D Lagrange polynomials whose anchor points reside on the element face, whereas all other 3-D Lagrange polynomials identically vanish on the element face. The number of 2-D polynomials is given by Ns = (N + 1)(N + 2)/2. We write for the normal flux:   \begin{eqnarray} \mathbf {f}_{n,f}(\mathbf {x},t) = \sum _{s=1}^{N_s} \mathbf {f}_{n,f}^s(t)\, l_s^{2D}(\mathbf {x}), \, \, \mathbf {x} \in \partial D^e, \end{eqnarray} (33)leading to the appearance of a face mass matrix $$M^{f}_{is}$$ defined as   \begin{eqnarray} M^{f}_{is} = \int _{I_{F}} l_i({\boldsymbol\xi }) l^{2D}_s({\boldsymbol\xi }) \text{d}\Sigma _{\mbox{ref}} = \frac{1}{J^f}\int _{\partial D^e} l_i(\mathbf {x}) l^{2D}_s(\mathbf {x}) \text{d}\Sigma \end{eqnarray} (34)for every face f of the element De with associated reference face IF and J f denoting the Jacobian for the face. To obtain the final semi-discrete numerical scheme for the integral form of eq. (15), we multiply it by the inverse mass matrix of the reference element (eq. 26) leading to (using summation convention)   \begin{eqnarray} \frac{\partial \mathbf {u}^m(t)}{\partial t} &=& \left(\mathbf {A}\frac{\partial \xi _l}{\partial x}+\mathbf {B} \frac{\partial \xi _l}{\partial y}+ \mathbf {C}\frac{\partial \xi _l}{\partial z}\right)(M^{-1}D^T_{\xi _l}M^T)_{mj}\mathbf {u}^j(t) \nonumber \\ &&-\frac{J^f}{J^e}\sum _{f=1}^{4}(M^{-1}M^{f})_{ms}\mathbf {f}^s_{n,f}(t) +\mathbf {s}^m(t) . \end{eqnarray} (35)For the integral form eq. (16) we obtain the following semi-discrete numerical scheme (using summation convention):   \begin{eqnarray} \frac{\partial \mathbf {u}^m(t)}{\partial t} &=& -\left(\mathbf {A}\frac{\partial \xi _l}{\partial x}+\mathbf {B}\frac{\partial \xi _l}{\partial y}+ \mathbf {C}\frac{\partial \xi _l}{\partial z}\right)(D_{\xi _l})_{mj}\mathbf {u}^j(t) \nonumber \\ &&-\frac{J^f}{J^e}\sum _{f=1}^{4}(M^{-1}M^{f})_{ms}\mathbf {r}^s_{n,f}(t)+\mathbf {s}^m(t) , \end{eqnarray} (36)where the $$\mathbf {r}^s_{n,f}$$ are now nodal coefficients of the normal Riemann flux with respect to 2-D Lagrange polynomials defined on the element faces. The two equations above represent the numerical schemes of the NDG method. They can be efficiently embedded into a numerical implementation and solved for the velocity–stress vector u. A similar equation can be derived for the 2-D case, where the system reduces to a system of five equations, to be solved for a velocity–stress vector with components u = (σxx, σyy, σxy, u, v)T. 2.6 Rules for spatial discretization The choice of element size, that is, the edge length of an element, is governed by the polynomial order of the Lagrange expansion, N, and the smallest significant wavelength, λmin, to be propagated stably across the model or parts of the model. According to our experience, we need 5 anchor points per minimum wavelength. Since the edge length of an element, ℓ, is about N times the spacing between anchor points, we get   $$\ell \le N\frac{\lambda _\mathrm{min}}{5} .$$ (37)The smallest wavelength is calculated from the highest significant frequency fmax and the smallest propagation velocity, vmin as   $$\lambda _\mathrm{min} = \frac{v_\mathrm{min}}{f_\mathrm{max}},$$ (38)leading finally to   $$\ell \le \frac{N}{5}\frac{v_\mathrm{min}}{f_\mathrm{max}}.$$ (39)The general procedure is to select one or several (for model parts) guiding values for element size following eq. (39) according to which the meshing software constructs a mesh. This mesh is then checked for stability again using eq. (39) but now solved for the highest frequency for which a stable simulation is possible:   $$f_\mathrm{stable} =\frac{N}{5}\frac{v_\mathrm{min}}{\ell }.$$ (40)The mesh passes the check if fstable ≤ fmax for all elements. If some elements do not satisfy this condition, a new mesh with a smaller guiding value for element size is constructed and the check is repeated. 2.7 Extension to anelasticity To incorporate anelasticity into our scheme, we follow the approach of Käser et al. (2007) with an adaptation to our way of flux computation. Based on original work by Emmerich & Korn (1987), it is assumed that the viscoelastic rheology can be described by a generalized Maxwell body mechanically represented by a spring in parallel to a series of Maxwell bodies each of which is made up of a spring and a dashpot. The stress–strain relation then generalizes to the form   \begin{eqnarray} \sigma _i(t) = M^u_{ij}\left(\epsilon _j(t)-\sum _{r=1}^{N_r}Y_r\omega _r\int _{-\infty }^t\epsilon _j(\tau )e^{-\omega _r(t-\tau )}{\rm d}\tau \right), \end{eqnarray} (41)where the σi stand for the six stresses, the εj for the corresponding strains, and the $$M^u_{ij}$$ for unrelaxed elastic moduli. Each of the Nr Maxwell bodies is characterized by a relaxation frequency ωr and an anelasticity coefficient Yr. For an isotropic medium, all moduli can be expressed through the Lamé constants λ and μ. Integration of this relation into the velocity–stress formulation is done by defining new material independent memory variables related to the anelastic stresses (Moczo et al. 2007):   \begin{eqnarray} \phi ^r_i(t) = \omega _r\frac{\partial }{\partial t}\int _{-\infty }^t\epsilon _i(\tau )e^{-\omega _r(t-\tau )}\text{d}\tau , \quad 1\le r\le N_r . \end{eqnarray} (42)According to eq. (41), all equations of the elastic velocity–stress formulation (eq. 1) containing time derivatives of stress will be extended by terms containing these new variables. The memory variables satisfy the differential equations   \begin{eqnarray} \frac{\partial \phi ^r_i(t)}{\partial t}+\omega _r\phi _i^r(t) = \omega _r\frac{\partial \epsilon _i}{\partial t} , \quad 1\le r\le N_r , \end{eqnarray} (43)which can be appended to the original velocity–stress eq. (1) because the strain-derivative on the right-hand side can be expressed by spatial derivatives of particle velocities. For each Maxwell body, six equations of this type are added to the system. Owing to the occurrence of spatial derivatives of particle velocities via the strain rates, the Jacobi matrices need to be extended to a dimension of 9 + 6Nr by adding 6Nr zero columns and 6Nr rows with nine entries each which are only non-zero if associated with particle velocity components. A modification of the Jacobi matrices affects the computation of the Riemann normal fluxes which were based (after rotation) on eigenvectors of the Jacobi matrix A. Here, we show that the coefficients of eqs. (23) and (24) describing amplitudes of waves propagating away from an element face remain unchanged and can be computed in exactly the same way as before. However, the eigenvectors change because these waves now also carry contributions in the new anelastic stress variables. The extended Jacobi matrix takes the following block structure:   \begin{eqnarray} \mathbf {A}^{\prime } = \left( \begin{array}{c{@}{\quad}c}\mathbf {A} & \mathbf {Z}_1 \\ \mathbf {P} & \mathbf {Z}_2 \end{array} \right) , \end{eqnarray} (44)where P is a 6Nr × 9 matrix, Z1 is a 9 × 6Nr zero matrix and Z2 is a 6Nr × 6Nr square zero matrix. Owing to the zero matrices, the extended Jacobi matrix has the same rank as the original one. Hence, all additional 6Nr eigenvalues must vanish. Now consider the case where the vector vk is an eigenvector of A with non-zero eigenvalue. This happens according to our numbering for index ranges 1−3 and 7−9. Then, there is an eigenvector of the extended Jacobi matrix, $$\mathbf {v}^{\prime }_k$$, satisfying $$\mathbf {A}^{\prime }\mathbf {v}^{\prime }_k = c_k\mathbf {v}^{\prime }_k$$, given by   \begin{eqnarray} \mathbf {v}^{\prime }_k = \left(\mathbf {v}_k,\mathbf {q}_k\right)^T\quad \mbox{with}\quad \mathbf {q}_k = \frac{1}{c_k}\mathbf {P}\mathbf {v}_k . \end{eqnarray} (45)If vk is an eigenvector of A but with a zero eigenvalue, then we can choose an eigenvector of A΄ of the form $$\mathbf {v}^{\prime }_k = (\mathbf {v_k},\mathbf {0})^T$$. This requires Pvk = 0, which is satisfied because the eigenvectors 4–6 of A vanish on the particle velocity components where P is non-zero. Finally, the remaining 6Nr eigenvectors with zero eigenvalue can be chosen as $$\mathbf {v}^{\prime }_k=(\mathbf {0},\mathbf {e}_k)^T$$ where ek is non-zero on the kth component only. An extended velocity–stress discontinuity is now decomposed as   \begin{eqnarray} \mathbf {U}^{\prime (1)}-\mathbf {U}^{\prime (2)} &=& \sum _{k=1}^3\gamma _k\left(\begin{array}{c}\mathbf {v}^{(1)}_k \\ \mathbf {q}^{(1)}_k\end{array}\right) +\sum _{k=4}^6\gamma _k\left(\begin{array}{c}\mathbf {v}_k \\ \mathbf {0}\end{array}\right) \nonumber \\ &&+\sum _{k=7}^9\gamma _k\left(\begin{array}{c}\mathbf {v}^{(2)}_k \\ \mathbf {q}^{(2)}_k\end{array}\right) +\sum _{k=10}^{9+6N_r}\gamma _k\left(\begin{array}{c}\mathbf {0} \\ \mathbf {e}_k\end{array}\right). \end{eqnarray} (46)However, considering only the first nine components, we arrive at exactly the same equation as eq. (18) from which the desired coefficients γk, 1 ≤ k ≤ 9 can be computed. Still, the associated eigenvectors changed. Hence, the extended normal flux on the right side of an interface now reads   \begin{eqnarray} \mathbf {r}^{2^{\prime }}_n = -\sum _{k=7}^9 c_k\gamma _k\left(\begin{array}{c}\mathbf {v}^{(2)}_k \\ \mathbf {q}^{(2)}_k\end{array}\right) . \end{eqnarray} (47) 2.8 Time discretization The time discretization is, in general, independent of the treatment of the spatial derivatives. Hence, established and well-proven methods for the time integrations can be used. For example, Komatitsch & Tromp (1999) use the classical second-order Newmark scheme to advance the SEM in time. For higher orders, the Runge–Kutta methods (Jameson et al. 1981) are available. Käser & Dumbser (2006) introduced the ADER-approach to the DG method for seismic problems by using the Cauchy–Kovalevsky theorem to solve the problem in the same order of approximation with respect to spatial and temporal variations. To advance the numerical schemes eqs (35) and (36) in time we apply a Total Variation Diminishing (TVD) Runge–Kutta method (Gottlieb & Shu 1998) which is third order accurate. It has three stages but needs only one additional array to store the velocity–stress–field for one physical integration step. The scheme has the form:   \begin{eqnarray} y^{n+1} & =& \frac{1}{3}y^n + \frac{2}{3}(z^{(2)} + \Delta t f(z^{(2)}))\quad \text{with} \nonumber \\ z^{(2)} &= &\frac{3}{4}y^n + \frac{1}{4}(z^{(1)} + \Delta t f(z^{(1)})),\quad \text{and} \nonumber \\ z^{(1)} &= &y^n + \Delta t f(y^n) . \end{eqnarray} (48)It computes the solution for the time step n + 1 going through two intermediate steps with the help of the auxiliary fields z(1) and z(2). Δt is the time step which has to satisfy the Courant–Friedrichs–Lewy criterion (Courant et al. 1928) according to which the velocity of the propagating waves must be smaller than the ratio of minimum point distance Δx divided by the time step Δt. The maximum time step is thus given by   \begin{eqnarray} \Delta t \le C_{\text{CFL}} \frac{\Delta x}{c_{\text{max}}}. \end{eqnarray} (49)Stability is ensured for CCFL ≤ 1. Typically CCFL is chosen smaller than 0.4 to run stable simulations on deformed meshes. The constant highly depends on the quality of the mesh and can be different for 2-D and 3-D simulations. Also, boundary conditions influence the value of the Courant constant. 3 NUMERICAL IMPLEMENTATION The NDG method described above was implemented into a software package called NEXD (Nodal Discontinuous Galerkin Finite Element in X Dimensions). As programming language, we use modern FORTRAN (FORTRAN95 and later) with object-oriented features. Parallelization is done by means of Open MPI (www.open-mpi.org). The code is organized in a way to allow a consistent work flow (Fig. 2). The model and the mesh can be generated with external meshing software. It is only required to provide plain text files containing the basic information about the mesh such as node numbering and connectivity, coordinates of vertices, boundary conditions and material parameters in the elements. A Python interface is available to use the Trelis/CUBIT software (www.csimsoft.com) for model and mesh generation. Three different program versions were developed, a 1-D, a 2-D and a 3-D one. Each program consists of three different main programs, the mesher, the solver and a post-processor to generate movie files. The mesher reads the mesh files and calculates all necessary information for the solver. It calculates the anchor points of the Lagrange polynomials and all necessary transformations, finds source and receiver positions and prepares the domain decomposition of the mesh with the help of the partitionizer ‘METIS’ (Karypis & Kumar 1999) for parallel simulations. All information is stored in a database available for each processor. Additional information for debugging and visualization purposes can be output on demand. Figure 2. View largeDownload slide The general procedure of a numerical simulation with NEXD. First, a mesh is created via a Python script and a meshing program such as Trelis/CUBIT. The ‘mesher’ reads in the mesh and prepares the forward calculation. If the simulation runs in parallel, the mesh is partitioned using the external library ‘Metis’. Database files with all required information about the mesh are created. These are read by the ‘solver’ which calculates the seismic wavefield. Parallelization of the solver is done with the help of MPI. Seismograms for displacement, velocity and acceleration are stored in plain ASCII files and databases are generated containing wavefield information. These files can be processed with the ‘movie’ program to create ‘VTK’ files to display videos of the propagating wavefield, for example, using the program ‘Paraview’. Figure 2. View largeDownload slide The general procedure of a numerical simulation with NEXD. First, a mesh is created via a Python script and a meshing program such as Trelis/CUBIT. The ‘mesher’ reads in the mesh and prepares the forward calculation. If the simulation runs in parallel, the mesh is partitioned using the external library ‘Metis’. Database files with all required information about the mesh are created. These are read by the ‘solver’ which calculates the seismic wavefield. Parallelization of the solver is done with the help of MPI. Seismograms for displacement, velocity and acceleration are stored in plain ASCII files and databases are generated containing wavefield information. These files can be processed with the ‘movie’ program to create ‘VTK’ files to display videos of the propagating wavefield, for example, using the program ‘Paraview’. The solver reads the database provided by the mesher and executes the forward simulation, generally in parallel mode using MPI. It controls the main time loop to advance the solution to the velocity–stress equations. For each time step, it iterates over all elements, deals with boundary conditions and calculates fluxes on element faces. Then, the velocity–stress vector is updated in time. At the end of the time loop, seismograms for displacement, velocity and acceleration are stored for every receiver position and component. On demand, snapshots of the wavefield are saved into databases which can be converted into VTK-files (Schroeder et al. 2006) for visualization with, for example ‘Paraview’ (Ayachit 2015). The code can handle different source time functions such as a Ricker or a Gaussian wavelet as well as externally provided source wavelets defined as time series. Both single force and moment tensor excitation are implemented. External velocity models defined on a regular grid can be mapped to an existing mesh by interpolation allowing a separation of meshing and assignment of material properties. Absorbing boundary conditions can be either realized by cancellation of incoming fluxes at the boundaries, or, much more efficiently, by using perfectly matching layers (PML). We have implemented a special variant of the PMLs in 2-D and 3-D, termed nearly perfectly matching layers (NPML) as described by Cummer (2003). Both 2-D and 3-D version of the code support the calculation of adjoint wavefields that may be used in the context of full waveform inversion. In addition, there is an interface to the full-waveform inversion code ASKI (Schumacher et al. 2016) for the computation of waveform sensitivity kernels. The code will be made available for download under GitHub (https://github.com/seismology-RUB, last accessed 1 December 2017). 4 NUMERICAL VALIDATION We validate our NDG approach for some test cases by comparing with analytical solutions, if available, and results from high-order SPECFEM simulations. 4.1 Lamb’s problem – 2-D elastic case The first test case is Lamb’s problem, the elastic response of a 2-D homogeneous half space to a single normal force at the free surface. The solution contains three main phases: a direct P-wave, a direct S-wave and a non-dispersive Rayleigh wave at the free surface. Lamb’s problem is a simple but challenging test to validate numerical codes especially for their accuracy with regard to dispersion. We compare solutions obtained with the NDG method to a semi-analytical reference solution, calculated with the program EX2DDIR (Berg et al. 1994). This program uses the Carniard-de Hoop (De Hoop 1960; Aki & Richards 1980) method to calculate the 2-D seismic response of a homogeneous half space. The model configuration and values for P-wave speed, S-wave speed and density of the homogeneous half-space are depicted in Fig. 3. Geometrical setting and material parameters exactly follow the configuration chosen by Käser & Dumbser (2006) and Komatitsch & Vilotte (1998). The free surface is inclined by an angle of α = 10° to avoid the Rayleigh wave propagating parallel to a coordinate axis. Due to the inclination, the solution is expected to be more sensitive to the numerical scheme and the implementation. The single force acting normal to the free surface resides at $$\mathbf {s} = (1720 \, \text{m},2303\, \text{m})$$ and two receivers are placed at the free surface at a distance of 990 m and 1706 m from the source. Receiver components are oriented normally and tangentially to the tilted free surface. As a source time function, a Ricker wavelet with a central frequency of fc = 14.5 Hz is used beginning at $$t_0 = -\frac{1.2}{fc}\, \text{s}$$ to ensure that the main peak of the wavelet is located at zero time. Figure 3. View largeDownload slide Model setup for computing the response of a homogeneous elastic half-space to a single normal force at the free surface, also known as Lamb’s problem. The free surface is inclined to avoid surface wave propagation parallel to a coordinate axis. Values of elastic wave speeds as well as position of source and receivers are shown together with the geometry of the model area. Figure 3. View largeDownload slide Model setup for computing the response of a homogeneous elastic half-space to a single normal force at the free surface, also known as Lamb’s problem. The free surface is inclined to avoid surface wave propagation parallel to a coordinate axis. Values of elastic wave speeds as well as position of source and receivers are shown together with the geometry of the model area. For the NDG method, the volume is discretized by triangular elements with maximum edge length of 45 m leading to a total of 17 700 elements for the mesh. A polynomial degree of five is used for the spatial resolution combined with a TVD Runge–Kutta scheme of third-order accuracy for time integration. The time step is Δt = 1.5 × 10−4 s and Nt = 10 000 time steps are calculated. On the left, right and bottom boundaries of the model, absorbing boundary conditions are used to mimic an infinitely wide half-space. A snapshot of the displacement field at t = 0.575 s is shown in Fig. 4. P- and S-waves propagate away from the source position followed by a Rayleigh wave with high amplitudes near the surface. A comparison of the resulting synthetic seismograms at the two receiver positions with the semi-analytical solution computed with EX2DDIR is depicted in Figs 5(a) and (b) for the normal component and Figs 5(c) and (d) for the tangential component. The seismograms fit very well. Maximum amplitude differences of less than 2 percent are visible in the high-amplitude Rayleigh wave train. This test validates the NDG method for accurate simulation of body waves and, in particular, surface waves. Figure 4. View largeDownload slide Snapshot of the absolute amplitude of the displacement field for Lamb’s problem after t = 0.575 s. The single force source is indicated as red dot, the two receivers as green dots. Three main phases are visible: a direct P-wave followed by an S-wave of smaller wave length due to the lower wave speed and a high amplitude surface wave near the free surface. Figure 4. View largeDownload slide Snapshot of the absolute amplitude of the displacement field for Lamb’s problem after t = 0.575 s. The single force source is indicated as red dot, the two receivers as green dots. Three main phases are visible: a direct P-wave followed by an S-wave of smaller wave length due to the lower wave speed and a high amplitude surface wave near the free surface. Figure 5. View largeDownload slide Displacement seismograms for Lamb’s problem for the normal and tangential component at receiver r1 (a, c) and receiver r2 (b, d). The analytical solution obtained with EX2DDIR is plotted in black. Results of the NDG simulation are plotted in grey. The difference is plotted as dashed line. A spatial order of N = 5 is used to obtain the NDG solution. The seismograms agree very well, validating the NDG for the 2-D case. The maximum relative difference occurring on the tangential component of receiver 1 is 2 per cent. No phase shift is notable for the two distant receivers, indicating that the NDG scheme calculates surface waves accurately. Figure 5. View largeDownload slide Displacement seismograms for Lamb’s problem for the normal and tangential component at receiver r1 (a, c) and receiver r2 (b, d). The analytical solution obtained with EX2DDIR is plotted in black. Results of the NDG simulation are plotted in grey. The difference is plotted as dashed line. A spatial order of N = 5 is used to obtain the NDG solution. The seismograms agree very well, validating the NDG for the 2-D case. The maximum relative difference occurring on the tangential component of receiver 1 is 2 per cent. No phase shift is notable for the two distant receivers, indicating that the NDG scheme calculates surface waves accurately. 4.2 P-convergence Test A different 2-D simulation in a homogeneous, elastic half-space was carried out to test the convergence of the numerical solution with increasing polynomial order of interpolation. Model setup, material properties and positions of a single force and a receiver are depicted in Fig. 6. At the top of the model, a free surface is assumed whereas on the other sides absorbing boundary conditions are applied. The force acts in horizontal (x)-direction and the source time function is a Ricker wavelet with centre frequency of fc = 600 Hz. The setup allows to observe the direct P- and S-waves on both receiver components and additionally converted and reflected waves coming from the free surface (Fig. 7). Given the seismic wave speeds and frequency spectrum of the source, the wave length of the S-wave at the centre frequency is λc ≃ 2.5 m while the minimum wavelength at the upper end of the frequency spectrum is λmin ≃ 1 m. Thus, the S-wave, in this setting, propagates a distance of about 20 wavelengths between source and receiver. The mesh consists of 21 606 elements leading to 648 180 degrees of freedom for the second-order test case (N = 2) and up to 7 129 980 degrees of freedom for the 10th-order calculation (N = 10) (Table 1) For comparison, reference seismograms are computed with a 12th-order SPECFEM simulation. Synthetic seismograms for both components displayed in Fig. 8 show small differences between the low-order (N = 3) NDG and the SPECFEM reference simulation which disappear for the 10th-order NDG solution. Figure 6. View largeDownload slide Model setup, material properties and positions of horizontal single force and receiver used for the p-convergence test. Figure 6. View largeDownload slide Model setup, material properties and positions of horizontal single force and receiver used for the p-convergence test. Figure 7. View largeDownload slide Snapshot of the absolute amplitude of the displacement field for the p-convergence test model. Time is T = 0.037 s and a polynomial order of N = 5 was used. Source position is marked by a red dot, receiver position by a green one. The mesh consists of 21 606 triangular elements. The direct P- and S-phases as well as the PP, PS and SP conversions from the free surface have already passed the receiver while the SS and PPS-phases are about to reach the receiver. Figure 7. View largeDownload slide Snapshot of the absolute amplitude of the displacement field for the p-convergence test model. Time is T = 0.037 s and a polynomial order of N = 5 was used. Source position is marked by a red dot, receiver position by a green one. The mesh consists of 21 606 triangular elements. The direct P- and S-phases as well as the PP, PS and SP conversions from the free surface have already passed the receiver while the SS and PPS-phases are about to reach the receiver. Figure 8. View largeDownload slide Seismograms of the displacement field for the (a) x- and (b) z-component for the p-convergence test. In both figures, the SPECFEM 12th-order reference solution is shown in black. The NDG solutions are green for the 3rd–order calculation and red for the 10th-order result. Whereas the third-order solution still has visible misfit, compared to the SPECFEM reference seismograms, no visible misfit can be detected with the N = 10 order NDG seismograms. Figure 8. View largeDownload slide Seismograms of the displacement field for the (a) x- and (b) z-component for the p-convergence test. In both figures, the SPECFEM 12th-order reference solution is shown in black. The NDG solutions are green for the 3rd–order calculation and red for the 10th-order result. Whereas the third-order solution still has visible misfit, compared to the SPECFEM reference seismograms, no visible misfit can be detected with the N = 10 order NDG seismograms. Table 1. Number of degrees of freedom for different polynomial orders of interpolation for a mesh with 21 606 triangular elements Order  Points per element  Degrees of freedom  2  6  648 180  3  10  1 080 300  4  15  1 620 450  5  21  2 268 630  6  28  3 024 840  8  45  4 861 350  10  66  7 129 980  Order  Points per element  Degrees of freedom  2  6  648 180  3  10  1 080 300  4  15  1 620 450  5  21  2 268 630  6  28  3 024 840  8  45  4 861 350  10  66  7 129 980  View Large To quantify the accuracy of NDG simulation with different polynomial orders, a misfit between the NDG and SPECFEM seismograms, expressed as L2 norm is calculated (Fig. 9). The misfit quickly decreases with polynomial order and converges to a minimum residual. Differences become very small from 5th order onwards. With regard to the numerical effort, it is sufficient to use fourth- or fifth-order polynomial interpolation for acceptably accurate simulations. Figure 9. View largeDownload slide L2 misfit between a 12-order SPECFEM simulation and NDG simulations of polynomial orders ranging from 2 to 10 for the p-convergence test case. As expected, the misfit decreases with increasing order and the numerical solution converges to the reference solution. Figure 9. View largeDownload slide L2 misfit between a 12-order SPECFEM simulation and NDG simulations of polynomial orders ranging from 2 to 10 for the p-convergence test case. As expected, the misfit decreases with increasing order and the numerical solution converges to the reference solution. 4.3 A 3-D simulation in a two-layer half-space This test simulation is performed with our 3-D implementation of the NDG method. The reference solution is, as before, computed with SPECFEM. The test shows the capability of the 3-D implementation and validates the NDG method with respect to 3-D simulations in discontinuous elastic media. Besides that, the potential of NDG to employ velocity adapted meshes is demonstrated. The considered model is a 3-D box of dimensions 100 m × 100 m× 50 m. Absorbing boundary conditions are used to simulate an infinitely wide model at all boundaries except at the top boundary, where a free surface boundary condition is assumed. The model is composed of a layer over a faster half-space. The top layer has a thickness of z = 12.5 m (Fig. 10). The material properties for this model are listed in Table 2. Figure 10. View largeDownload slide Two different meshes for a two layer simulation. Model dimensions are x = 100 m, y = 100 m and z = 50 m. The top layer has a thickness of z = 12.5 m and velocities vp = 3000 m s−1, vs = 1700 m s−1 and density ρ = 2000 kg m−3. The bottom layer is faster with material properties vp = 5500 m s−1, vs = 3400 m s−1 and ρ = 2000 kg m−3. The top surface is a free boundary. Absorbing boundary conditions are used for the other sides to mimic an infinite half-space. Model (a) has an element size of h = 2.5 m leading to 299 425 elements in total. The coarsened mesh in (b) has an element size of 5 m and 104 414 elements in total. Figure 10. View largeDownload slide Two different meshes for a two layer simulation. Model dimensions are x = 100 m, y = 100 m and z = 50 m. The top layer has a thickness of z = 12.5 m and velocities vp = 3000 m s−1, vs = 1700 m s−1 and density ρ = 2000 kg m−3. The bottom layer is faster with material properties vp = 5500 m s−1, vs = 3400 m s−1 and ρ = 2000 kg m−3. The top surface is a free boundary. Absorbing boundary conditions are used for the other sides to mimic an infinite half-space. Model (a) has an element size of h = 2.5 m leading to 299 425 elements in total. The coarsened mesh in (b) has an element size of 5 m and 104 414 elements in total. Table 2. Elastic properties of the two-layer half-space example in Section 4.3. Layer  vp (m s−1)  vs (m s−1)  ρ (kg m−3)  Yellow  3000  1700  2000  Green  5500  3400  2000  Layer  vp (m s−1)  vs (m s−1)  ρ (kg m−3)  Yellow  3000  1700  2000  Green  5500  3400  2000  View Large A vertical single force source is located in the half-space, at coordinates $$\mathbf {s} = (30\, \text{m}, 30\, \text{m}, 25\, \text{m})^T$$, 12.5 m below the bottom of the top layer. Note that in this example, the vertical coordinate z points bottom up. As source time function, a Ricker wavelet with dominant frequency of fc = 160 Hz is used. We perform two simulations with the NDG code: one on a mesh with approximately uniformly sized tetrahedral elements of an edge length of around h = 2.5 m leading to 299 425 elements and 94.3 × 106 degrees of freedom in total, and one on a mesh that is coarsened in the bottom layer to an edge length of around 5 m leading to 104 414 elements and 32.9 × 106 degrees of freedom in total. In the first simulation, we perform 16 000 time steps of length Δt = 6.415 × 10−6 s whereas in the second one, we run through 12 600 time steps of length Δt = 8.38 × 10−6 s. Polynomial order of interpolation is N = 4 in both simulations. Computation time of the first simulation is by a factor of 3.5 greater than that of the second one. For comparison, we run the same model with SPECFEM using hexahedral elements and polynomial order of 4. Here, 4500 time steps of length Δt = 2.34 × 10−5 s are performed. Element edge length is 2 m leading to 148 137 elements with 30 × 106 degrees of freedom in total. To record the wavefield, nine receivers are placed at the surface along a line parallel to the x-axis from $$\mathbf {r}_1=(10\, \text{m},50\, \text{m},50\, \text{m})^T$$ to $$\mathbf {r}_{9}=(90\, \text{m},50\, \text{m},50\, \text{m})^T$$. Synthetic seismograms obtained by all three simulations at the receivers are depicted in Fig. 11. No differences are visible on the scale of the figure. A quantitative calculation yields normalized rms-misfits for all seismograms of 0.2  per cent between SPECFEM and NDG with the uniform mesh, of 0.2  per cent between SPECFEM and NDG with the coarsened mesh and of 0.04  per cent between NDG with the uniform and NDG with the coarsened mesh. This test validates the 3-D-version of our DG method. Adapting the mesh to velocity changes in the model can significantly help to reduce computation times without loosing accuracy, and is easily realized with the DG method and tetrahedral elements. Figure 11. View largeDownload slide Seismic section of the z-component of particle velocity for the 3-D two-layer test. The line of nine receivers is located at the top of the model and has endpoint coordinates $$\mathbf {r}_1=(10\, \text{m},50\, \text{m},50\, \text{m})^T$$ and $$\mathbf {r}_{9}=(90\, \text{m},50\, \text{m},50\, \text{m})^T$$. The source is in the fast layer at $$\mathbf {s} = (30\, \text{m}, 30\, \text{m}, 25\, \text{m})^T$$. Black line: SEM reference simulation, dark grey line: NDG simulation with the coarsened mesh in the bottom layer. Light grey line: NDG simulation with about equal element size in both layers. For each receiver, the seismograms are normalized to the maximum of the SEM seismogram. Seismograms are shifted relative to each other. Figure 11. View largeDownload slide Seismic section of the z-component of particle velocity for the 3-D two-layer test. The line of nine receivers is located at the top of the model and has endpoint coordinates $$\mathbf {r}_1=(10\, \text{m},50\, \text{m},50\, \text{m})^T$$ and $$\mathbf {r}_{9}=(90\, \text{m},50\, \text{m},50\, \text{m})^T$$. The source is in the fast layer at $$\mathbf {s} = (30\, \text{m}, 30\, \text{m}, 25\, \text{m})^T$$. Black line: SEM reference simulation, dark grey line: NDG simulation with the coarsened mesh in the bottom layer. Light grey line: NDG simulation with about equal element size in both layers. For each receiver, the seismograms are normalized to the maximum of the SEM seismogram. Seismograms are shifted relative to each other. 4.4 Homogeneous full space—2-D anelastic case To validate our anelastic simulations we compare them to an analytical solution given by Carcione et al. (1988). They derived Green functions for a homogeneous 2-D full space in the frequency domain. Anelasticity is introduced by using complex-valued, frequency-dependent P- and S-wave velocities calculated from complex elastic moduli. A back transformation to the time domain results in the displacement field at the receiver position. For comparison with our DG method, we choose a homogeneous, quadratic model domain of 2000 m × 2000 m with a density of ρ = 2000 kg m−3. At a reference frequency of f = 50 Hz, the real part of P wave velocity is set to vp = 3000 m s−1 and that of S wave velocity to vs = 2000 m s−1. Anelasticity coefficients and relaxation frequencies are chosen in a way to achieve nearly constant values of the quality factors around Qκ = 30 and Qμ = 20. A vertical point force is located at the centre of the domain at x = y = 1000 m and a receiver resides at x = 1500 m, y = 500 m. Anelasticity is modelled, as described by Käser et al. (2007), using 3 Maxwell bodies leading to frequency dependent wave velocities and quality factors shown in Fig. 12. To evaluate the analytical expressions for the Green functions, we use exactly the same dispersion relation. The analytical Green function is convolved with the source time function used by Carcione et al. (1988) who chose the product of a Gaussian and a cosine function. Comparisons of the anelastic and elastic DG solutions with their analytical counterparts are shown in Fig. 13. For the elastic case, the values of P and S wave velocities at the reference frequency were used. In both cases, the numerical and analytical solution match very well. The relative difference between numerical and analytic solution at the extremal points is less than 10−3. Figure 12. View largeDownload slide Velocities for P and S waves (left) and quality factors, Qκ and Qμ, (right) versus frequency of the anelastic DG model. Figure 12. View largeDownload slide Velocities for P and S waves (left) and quality factors, Qκ and Qμ, (right) versus frequency of the anelastic DG model. Figure 13. View largeDownload slide Comparison of nodal DG simulation with an analytical solution after Carcione et al. (1988) for the anelastic and elastic cases. Traces are shifted relative to each other for better visibility. Figure 13. View largeDownload slide Comparison of nodal DG simulation with an analytical solution after Carcione et al. (1988) for the anelastic and elastic cases. Traces are shifted relative to each other for better visibility. 5 A COMPLEX NUMERICAL EXAMPLE This example case is designed to demonstrate the capability of our DG approach to model wave propagation in media exhibiting complicated geological structures and, in addition, containing bodies of non-trivial geometric forms. It is taken from related work in the field of tunnel seismics where seismic waves are used to image disturbing structures in front of the tunnel boring machine. The model is composed of a curved cylindrical tunnel embedded into three sedimentary layers which are offset relative to each other by two thrust faults (Fig. 14). The tunnel axis lies in the xy-plane. At the entry, the tunnel axis is parallel to the y-axis but later turns away to the right, that is, towards increasing x-values. The vertical coordinate, z, of the tunnel axis is constant. The interfaces between the sedimentary layers are slightly inclined relative to the plane of the tunnel axis and intersect the tunnel axis at some distance from the tunnel entry implying different material properties above and below the tunnel from there on. In addition, the thrust faults, inclined relative to the vertical, intersect the tunnel at some distance away from the entry. The elastic material parameters of the sedimentary layers are rather different and are listed in Table 3. Two high-velocity layers at the top and bottom of the model enclose a low-velocity layer in between. Figure 14. View largeDownload slide Model motivated by tunnel seismics containing a curved cylindrical tunnel drilled into three inclined sedimentary layers offset relative to each other by thrust faults. The tunnel is l = 50 m long and ends in the fast third layer (red) with material properties vp = 4650 m s−1, vs = 2700 m s−1 and ρ = 2600 kg m−3. This layer is covered by a slow layer (grey) with vp = 2000 m s−1, vs = 1200 m s−1 and ρ = 2350 kg m−3. The top layer (blue) is, again, a fast layer having a P-wave velocity of vp = 4330 m s−1 a S-wave velocity of vs = 2500 m s−1 and a density of ρ = 2500 kg m−3. The overall dimensions of the model are wx = 64 m, wy = 100 m and wz = 64 m. (a) Entire model box, (b) a horizontal section to visualize the track of the tunnel. The tetrahedral mesh contains Ne = 1 541 398 elements resulting in a simulation with 4.856 × 108 degrees of freedom. Figure 14. View largeDownload slide Model motivated by tunnel seismics containing a curved cylindrical tunnel drilled into three inclined sedimentary layers offset relative to each other by thrust faults. The tunnel is l = 50 m long and ends in the fast third layer (red) with material properties vp = 4650 m s−1, vs = 2700 m s−1 and ρ = 2600 kg m−3. This layer is covered by a slow layer (grey) with vp = 2000 m s−1, vs = 1200 m s−1 and ρ = 2350 kg m−3. The top layer (blue) is, again, a fast layer having a P-wave velocity of vp = 4330 m s−1 a S-wave velocity of vs = 2500 m s−1 and a density of ρ = 2500 kg m−3. The overall dimensions of the model are wx = 64 m, wy = 100 m and wz = 64 m. (a) Entire model box, (b) a horizontal section to visualize the track of the tunnel. The tetrahedral mesh contains Ne = 1 541 398 elements resulting in a simulation with 4.856 × 108 degrees of freedom. Table 3. Elastic properties of the sedimentary layers in the tunnel seismics model Layer  vp (m s−1)  vs (m s−1)  ρ (kg m−3)  Blue  4330  2500  2500  Grey  2000  1200  2350  Red  4650  2700  2600  Layer  vp (m s−1)  vs (m s−1)  ρ (kg m−3)  Blue  4330  2500  2500  Grey  2000  1200  2350  Red  4650  2700  2600  View Large The dimensions of the entire model box are 64 m × 100 m × 64 m. The tunnel is approximately 50 m long, ending in the fast bottom layer. Its diameter is 12 m. The model is meshed with tetrahedral elements for the NDG simulation with an edge size of h = 2 m for both upper and lower layer and with a size of h = 1 m for the middle slow layer. This choice allows stable simulations up to fc = 500 Hz. To save computational effort, some mesh coarsening was applied in the outer regions of the mesh. In total, Ne = 1 541 398 tetrahedral elements are required to discretize the model. Polynomial order is 4 and number of Lagrange anchor points Np is 35. With Nu = 9 field variables we end up with NpNuNe = 4.8 × 108 degrees of freedom for this simulation. Absorbing boundary conditions based on flux cancellation are applied on all boundaries except for the top boundary where a free surface is assumed. We also attempted to mesh this model with hexahedral elements, but failed to produce a mesh of acceptable quality. For this reason, a computation with SPECFEM could not be performed. Seismic waves in the model are excited by a single force residing at the sidewall of the tunnel and acting normally to it. Its precise location is $$\mathbf {x}_s = (34.53 \, \text{m}, 37.73 \, \text{m}, 32.00\, \text{m})$$. The source time function is a Ricker wavelet with centre frequency of 500 Hz. The time step is Δt = 3.414 × 10−6 s and Nt = 20 000 time steps were calculated. A line of 46 receivers is placed on the right tunnel sidewall (viewed from the tunnel entry). Receiver spacing is 1 m (Fig. 14). The simulation was carried out on 192 CPU cores. Different stages of the wavefield are displayed as snapshots in the xy-plane in Fig. 15. The source emits P- and S-waves and a high-amplitude surface wave attached to the tunnel sidewall. The wave speed discontinuity due to the first thrust fault is well recognized in Fig. 15(c). When the tunnel Rayleigh wave detaches from the tunnel front face, it converts into an S-wave which propagates away from the tunnel (Figs 15c–f). When hitting the second thrust fault (Fig. 15e), the curvature of the transmitted S-wave front changes and a reflection is produced (Fig. 15f). On arrival at the tunnel front face, this reflection converts into a tunnel Rayleigh wave (Fig. 15g) and runs along the entire tunnel side wall (Figs 15h–l). This scenario of a tunnel Rayleigh wave converting into an S-wave and back into a tunnel Rayleigh wave after being reflected has been described earlier by Bohlen et al. (2007). Figure 15. View largeDownload slide Snapshots of the wavefield for the curved tunnel model. The absolute value of the displacement field is shown in the xy−plane (see also Fig. 14b). The source is marked as red dot while the receivers are highlighted by green points. The time between the individual frames is Δt = 3.41 × 10−3 s. A single force normal to the tunnel sidewall (a) excites a high amplitude surface wave propagating along the tunnel (b) which is converted into an S-wave (c–e) propagating ahead of the tunnel. A reflection at the second thrust fault occurs (f) and a part of the signal propagates back to the tunnel where it is converted to a tunnel surface wave (g). This signal can be observed in the seismograms and propagates all along the free surface of the tunnel (h)–(l). Figure 15. View largeDownload slide Snapshots of the wavefield for the curved tunnel model. The absolute value of the displacement field is shown in the xy−plane (see also Fig. 14b). The source is marked as red dot while the receivers are highlighted by green points. The time between the individual frames is Δt = 3.41 × 10−3 s. A single force normal to the tunnel sidewall (a) excites a high amplitude surface wave propagating along the tunnel (b) which is converted into an S-wave (c–e) propagating ahead of the tunnel. A reflection at the second thrust fault occurs (f) and a part of the signal propagates back to the tunnel where it is converted to a tunnel surface wave (g). This signal can be observed in the seismograms and propagates all along the free surface of the tunnel (h)–(l). Synthetic seismograms of the y-component at the receivers are shown in Fig. 16. The traces are ordered from bottom to top according to their numbering which starts at 1 on the tunnel front face. The source resides close to receiver 11 and the first thrust fault lies between receivers 19 and 20. Close to the source (receivers 1–19), we recognize a superposition of the S-wave and the Rayleigh wave (R1). P-waves are very weak owing to the fact that the force acts normally to the tunnel sidewall. When arriving at the first thrust fault, a converted P-wave is excited (P) and the tunnel Rayleigh wave changes speed (R2). In addition, the tunnel Rayleigh wave is reflected from this discontinuity and propagates back to the tunnel front face (R4 and R3). The Rayleigh wave generated by the S-wave coming back from the second thrust fault is denoted by R5. It arrives first at the tunnel front face and then propagates all along the receiver line. Phase L is a reflection from the interface above the tunnel whereas ‘surf’ is a reflection from the free surface. Phase ‘AB’ is a reflection from the boundary that contains the tunnel entry owing to a not perfectly working absorbing boundary condition. Figure 16. View largeDownload slide Seismic section of the y-component (along axis at tunnel entry) of the displacement field. Receivers are distributed at the right-hand side of the tunnel sidewall (see Fig. 15, green dots, and Fig. 14b, white dots). Receiver 1 (bottom of section) resides at the tunnel front face. Receiver 46 (top of section) sits at the tunnel entry. The source is very close to receiver 11. ‘R1’ denotes the direct surface wave excited by a single force acting normal to the tunnel sidewall. ‘R2’ is its continuation after having crossed the first thrust fault which intersects the tunnel. ‘R3’ is a surface wave reflected from the thrust fault. ‘R4’ is a surface wave reflected from the tunnel front face. R1 also excites a P wave (‘P’) when hitting the thrust fault. ‘R5’ is the so-called RSSR-wave, a surface wave generated through conversion from an S-wave which is reflected from the second thrust fault ahead of the tunnel front face. ‘L’ denotes a reflection from the layer interface above the tunnel while ‘surf’ is the signal from the top free surface. At the tunnel entry, a nonphysical reflection owing to imperfect absorbing boundary conditions evolves and is marked as ‘AB’. Figure 16. View largeDownload slide Seismic section of the y-component (along axis at tunnel entry) of the displacement field. Receivers are distributed at the right-hand side of the tunnel sidewall (see Fig. 15, green dots, and Fig. 14b, white dots). Receiver 1 (bottom of section) resides at the tunnel front face. Receiver 46 (top of section) sits at the tunnel entry. The source is very close to receiver 11. ‘R1’ denotes the direct surface wave excited by a single force acting normal to the tunnel sidewall. ‘R2’ is its continuation after having crossed the first thrust fault which intersects the tunnel. ‘R3’ is a surface wave reflected from the thrust fault. ‘R4’ is a surface wave reflected from the tunnel front face. R1 also excites a P wave (‘P’) when hitting the thrust fault. ‘R5’ is the so-called RSSR-wave, a surface wave generated through conversion from an S-wave which is reflected from the second thrust fault ahead of the tunnel front face. ‘L’ denotes a reflection from the layer interface above the tunnel while ‘surf’ is the signal from the top free surface. At the tunnel entry, a nonphysical reflection owing to imperfect absorbing boundary conditions evolves and is marked as ‘AB’. 6 DISCUSSION The DG method is extraordinarily suitable for wave propagation simulations in complicated geological environments due to its ability to work with tetrahedral elements. With tetrahedra, a high-quality mesh is much easier to achieve than with the hexahedral elements used by the SEM. Creating hexahedral meshes for complicated models may take significantly more time and experience than a comparable tetrahedral mesh. The curved tunnel example demonstrates that it is extremely difficult to produce a high-quality hexahedral mesh, while for the same model a tetrahedral mesh can be constructed semi-automatically. Another big advantage of tetrahedral meshes is the fact that element sizes can be easily adapted to the physical properties of the medium without affecting the accuracy of the simulation. In this way, adaptive mesh construction can partially compensate the higher numerical costs of the DG method. This advantage pays off particularly well in complex geological situations, which often occur in unconsolidated sediments and are thus relevant for near surface seismic modelling. Nevertheless, it is worthwhile to carefully consider which method is used for a given simulation problem. In many cases, the SEM may be the method of choice. The realization of absorbing or radiation boundary conditions by just omitting incoming fluxes at boundary elements is not fully satisfactory especially for waves with grazing incidence. For this reason, NPML boundary conditions, a variation of PML conditions, are implemented. Spurious reflections from the absorbing boundary nearly vanish for the simulations using NPML. However, PML conditions are only weakly stable for the seismic wave equation (Fichtner et al. 2011) leading to cases where the solution grows exponentially with time. Especially large material contrast close to the absorbing boundary may favour instable behaviour. This problem is overcome by a switching mechanism. When using NPML boundary conditions, the total energy of the medium is monitored and if the energy, after a certain time, begins to grow, the NPML is switched to the classical, flux-controlled boundary conditions thereby saving the simulation. Exponential growth of energy in the solution is detected with an STA/LTA trigger, originally developed for earthquake detection. Viscoelastic attenuation is modelled using a series of parallel Maxwell bodies. The comparison of a viscoelastic simulations in 2-D with an analytical solution exhibits very good agreement. The numerical example also shows that three parallel Maxwell bodies are sufficient to obtain nearly constant quality factors and a log-linear dependency of wave velocities across a wide frequency range. Introducing more mechanisms is normally not required and increases the numerical costs significantly. The current implementation of the NDG method successfully passes several accuracy tests. For example, synthetic seismograms for Lamb’s problem very well match seismograms computed from an analytical solution. Further comparisons with SEM also show very good agreement. The success of these tests also validates the parallel implementation of the NDG code using MPI. A convergence test, where simulations are done with increasing spatial order, underline the successive reduction of misfit to a reference solution with increasing order. It turns out that a spatial order of N = 4 is a viable compromise between numerical effort and accuracy for the NDG method. Very high order simulations require a correspondingly small time step making the computations inefficient. The scaling of the parallel implementation with increasing numbers of processors (not described in this paper) could only be tested on a single machine with 48 CPUs owing to the heterogeneity of the available computer cluster. At least for this machine, the computation time decreases quasi-linearly with CPU number. We expect that this behaviour is also preserved for a larger number of CPUs. 7 CONCLUSIONS The NDG method is adapted to 1-D, 2-D and 3-D elastic and anelastic seismic wave propagation. Formulations of numerical fluxes for the elastic and anelastic case are given which are derived from an exact solution of the heterogeneous Riemann problem, and honour possibly differing material properties in adjacent elements. Absorbing boundary conditions relying on the Nearly Perfectly Matched Layer approach are adapted to the NDG as an alternative to flux based boundary conditions. All features are implemented in the software package NEXD allowing efficiently parallelized numerical simulations. The accuracy and implementation is successfully validated through different test cases such as Lamb’s problem by either comparing with analytical solutions or spectral element simulations. The capabilities of the method are demonstrated by a challenging example case derived from tunnel reconnaissance with complicated geological structures and internal bodies of non-trivial geometric form. The NDG method is particularly attractive for such scenarios due to its high order approximation of the seismic wavefield as well as its use of triangular and tetrahedral meshes for the model discretization. The code may be used as forward solver for full waveform inversion as it offers the calculation of adjoint wavefields and provides an interface to the computation of waveform sensitivity kernels. The code is available for download under GitHub (https://github.com/seismology-RUB, last accessed 1 December 2017). ACKNOWLEDGEMENTS This work was supported by the Collaborative Research Center 837 ‘Interaction models in mechanized tunneling’ funded by the Deutsche Forschungsgemeinschaft (DFG). REFERENCES Aki K., Richards P.G., 1980. Quantivative Seismology, Theorie and Methods , University Science Books. Ayachit U., 2015. The ParaView Guide: A Parallel Visualization Application , Kitware. Bao H., Bielak J., Ghattas O., Kallivokas L., O’Hallaron D., Shewchuk J., Xu J., 1998. Large-scale simulation of elastic wave propagation in heterogeneous media on parallel computers, Comput. Methods Appl. Mech. Eng. , 152(1–2), 85– 102. https://doi.org/10.1016/S0045-7825(97)00183-7 Google Scholar CrossRef Search ADS   Bassi F., Rebay S., 1997. A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier–Stokes equations, J. Comput. Phys. , 131( 2), 267– 279. https://doi.org/10.1006/jcph.1996.5572 Google Scholar CrossRef Search ADS   Berg P., If F., Nielsen P., Skovgaard O., 1994. Analytical reference solutions, Model. Earth Oil Explor. , 1, 421– 427. Bohlen T., 2002. Parallel 3-D viscoelastic finite difference seismic modelling, Comput. Geosci. , 28, 887– 899. https://doi.org/10.1016/S0098-3004(02)00006-7 Google Scholar CrossRef Search ADS   Bohlen T., Lorang U., Rabbel W., Mueller C., Giese R., Lueth S., Jetschny S., 2007. Rayleigh-to-shear wave conversion at the tunnel face—from 3D-FD modeling to ahead-of-drill exploration, Geophysics , 72( 6), T67– T79. https://doi.org/10.1190/1.2785978 Google Scholar CrossRef Search ADS   Boxberg M. S., Heuel J., Friederich W., 2017. A nodal discontinuous Galerkin solver for modeling seismic wave propagation in porous media, in Poromechanics VI , eds Vandamme M., Dangla P., Pereira J-M., Ghabezloo S., American Society of Civil Engineers. Carcione J. M., Kosloff D., Kosloff R., 1988. Wave propagation simulation in a linear viscoelastic medium, Geophys. J. Int. , 95( 3), 597– 611. https://doi.org/10.1111/j.1365-246X.1988.tb06706.x Google Scholar CrossRef Search ADS   Chen X.-F., 2007. Generation and propagation of seismic SH waves in multi-layered media with irregular interfaces, Adv. Geophys. , 48, 191– 264. Google Scholar CrossRef Search ADS   Cockburn B., Karniadakis G.E., Shu C.-W., 2000. The Development of Discontinuous Galerkin Methods , Springer. Google Scholar CrossRef Search ADS   Cockburn B., Li F., Shu C.-W., 2004. Locally divergence-free discontinuous Galerkin methods for the Maxwell equations, J. Comput. Phys. , 194( 2), 588– 610. https://doi.org/10.1016/j.jcp.2003.09.007 Google Scholar CrossRef Search ADS   Courant R., Friedrichs K., Lewy H., 1928. Über die partiellen Differenzengleichungen der mathematischen Physik, Math. Ann. , 100( 1), 32– 74. https://doi.org/10.1007/BF01448839 Google Scholar CrossRef Search ADS   Cummer S.A., 2003. A simple, nearly perfectly matched layer for general electromagnetic media, IEEE Microw. Wirel. Compon. Lett. , 13( 3), 128– 130. https://doi.org/10.1109/LMWC.2003.810124 Google Scholar CrossRef Search ADS   Cupillard P., Delavaud E., Burgos G., Festa G., Vilotte J.-P., Capdeville Y., Montagner J.-P., 2012. Regsem: a versatile code based on the spectral element method to compute seismic wave propagation at the regional scale, Geophys. J. Int. , 188( 3), 1203– 1220. https://doi.org/10.1111/j.1365-246X.2011.05311.x Google Scholar CrossRef Search ADS   De Hoop A., 1960. A modification of Cagniard’s method for solving seismic pulse problems, Appl. Sci. Res., Sec. B , 8( 1), 349– 356. https://doi.org/10.1007/BF02920068 Google Scholar CrossRef Search ADS   de la Puente J., 2008. Seismic Wave Propagation for Complex Rheologies , VDM Verlag Dr. Müller. Delcourte S., Fezoui L., Glinsky-Olivier N., 2009. A high-order discontinuous Galerkin method for the seismic wave propagation, in ESAIM: Proceedings , vol. 27, pp. 70– 89, EDP Sciences. Google Scholar CrossRef Search ADS   Dumbser M., Käser M., 2006. An arbitrary high order discontinuous Galerkin method for elastic waves on unstructured meshes II: the three-dimensional isotropic case, Geophys. J. Int. , 167( 1), 319– 336. https://doi.org/10.1111/j.1365-246X.2006.03120.x Google Scholar CrossRef Search ADS   Emmerich H., Korn M., 1987. Incorporation of attenuation into time-domain computations of seismic wave fields, Geophysics , 52( 9), 1252– 1264. https://doi.org/10.1190/1.1442386 Google Scholar CrossRef Search ADS   Etienne V., Chaljub E., Virieux J., Glinsky N., 2010. An hp-adaptive discontinuous Galerkin finite-element method for 3-D elastic wave modelling, Geophys. J. Int. , 183( 2), 941– 962. https://doi.org/10.1111/j.1365-246X.2010.04764.x Google Scholar CrossRef Search ADS   Fichtner A., Kennett B.L., Igel H., Bunge H.-P., 2009. Full seismic waveform tomography for upper-mantle structure in the Australasian region using adjoint methods, Geophys. J. Int. , 179( 3), 1703– 1725. https://doi.org/10.1111/j.1365-246X.2009.04368.x Google Scholar CrossRef Search ADS   Fichtner A., Bleibinhaus F., Capdeville Y., 2011. Full Seismic Waveform Modelling and Inversion , Springer. Google Scholar CrossRef Search ADS   Furumura T., Kennett B., Furumura M., 1998. Seismic wavefield calculation for laterally heterogeneous whole earth models using the pseudospectral method, Geophys. J. Int. , 135( 3), 845– 860. https://doi.org/10.1046/j.1365-246X.1998.00682.x Google Scholar CrossRef Search ADS   Geodynamics.org, 2009. Computational infrastructure for geodynamics, http://www.geodynamics.org, last accessed 1 December 2017. Gottlieb S., Shu C.-W., 1998. Total variation diminishing Runge-Kutta schemes, Math. Comput. Am. Math. Soc. , 67( 221), 73– 85. https://doi.org/10.1090/S0025-5718-98-00913-2 Google Scholar CrossRef Search ADS   Heinecke A. et al.  , 2014. Petascale High Order Dynamic Rupture Earthquake Simulations on Heterogeneous Supercomputers, in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis SC14 , pp. 3– 14, IEEE, New Orleans, LA, USA, Gordon Bell Finalist. Hesthaven J.S., Warburton T., 2008. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications , vol. 54, Springer. Igel H., 2016. Computational Seismology—A Practical Introduction , Oxford Univ. Press. Google Scholar CrossRef Search ADS   Jameson A., Schmidt W., Turkel E., 1981. Numerical solution of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes, AIAA 14th Fluid and Plasma Dynamic Conference , June 23–25, 1981, Palo Alto, California, paper 1981-1259. Karypis G., Kumar V., 1999. A fast and highly quality multilevel scheme for partitioning irregular graphs, SIAM J. Sci. Comput. , 20, 359– 392. https://doi.org/10.1137/S1064827595287997 Google Scholar CrossRef Search ADS   Käser M., Dumbser M., 2006. An arbitrary high order discontinuous galerkin method for elastic waves on unstructured meshes. I: The two-dimensional isotropic case with external source terms, Geophys. J. Int. , 166( 2), 855– 877. https://doi.org/10.1111/j.1365-246X.2006.03051.x Google Scholar CrossRef Search ADS   Käser M., Dumbser M., de La Puente J., Igel H., 2007. An arbitrary high-order Discontinuous Galerkin method for elastic waves on unstructured meshes—III. Viscoelastic attenuation, Geophys. J. Int. , 168( 1), 224– 242. https://doi.org/10.1111/j.1365-246X.2006.03193.x Google Scholar CrossRef Search ADS   Komatitsch D., Tromp J., 1999. Introduction to the spectral-element method for 3-D seismic wave propagation, Geophys. J. Int. , 139( 3), 806– 822. https://doi.org/10.1046/j.1365-246x.1999.00967.x Google Scholar CrossRef Search ADS   Komatitsch D., Vilotte J.P., 1998. The spectral-element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures, Bull. seism. Soc. Am. , 88( 2), 368– 392. Komatitsch D., Vilotte J.P., Vai R., Castillo-Covarrubias J.M., Sánchez-Sesma F.J., 1999. The spectral-element method for elastic wave equations: application to 2D and 3D seismic problems, Int. J. Numer. Methods Eng. , 45( 9), 1139– 1164. https://doi.org/10.1002/(SICI)1097-0207(19990730)45:9<1139::AID-NME617>3.0.CO;2-T Google Scholar CrossRef Search ADS   Komatitsch D., Tsuboi S., Ji C., Tromp J., 2003. A 14.6 billion degrees of freedom, 5 teraflops, 2.5 terabyte earthquake simulation on the Earth Simulator, Proceedings of the ACM/IEEE Supercomputing SC’2003 conference , pp. 4– 11, Gordon Bell Prize winner article. LeVeque R.J., 2002. Finite Volume Methods for Hyperbolic Problems , vol. 31, Cambridge Univ. Press. Google Scholar CrossRef Search ADS   Li Y., 2011. Development of numerical simulation method for nonlinear elastodynamic: application to acoustic imaging of defect with the help of cavity chaotic transducer, PhD thesis , Ecole Centrale de Lille. Marfurt K., 1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave-equations, Geophysics , 49( 5), 533– 549. https://doi.org/10.1190/1.1441689 Google Scholar CrossRef Search ADS   Maufroy E. et al.  , 2015. Earthquake ground motion in the Mygdonian Basin, Greece: the E2VP verification and validation of 3D numerical simulation up to 4 Hz, Bull. seism. Soc. Am. , 105, 1398– 1418. https://doi.org/10.1785/0120140228 Google Scholar CrossRef Search ADS   Mercerat E.D., Glinsky N., 2015. A nodal high-order discontinuous Galerkin method for elastic wave propagation in arbitrary heterogeneous media, Bull. seism. Soc. Am. , 201, 1101– 1118. Moczo P., Robertsson J.O., Eisner L., 2007. The finite-difference time-domain method for modeling of seismic wave propagation, Adv. Geophys. , 48, 421– 516. Google Scholar CrossRef Search ADS   Moczo P., Kristek J., Galis M., 2014. The Finite-Difference Modelling of Earthquake Motions: Waves and Ruptures , Cambridge Univ. Press. Google Scholar CrossRef Search ADS   Patera A.T., 1984. A spectral selement method for fluid dynamics: laminar flow in a channel expansion, J. Comput. Phys. , 54, 468– 488. https://doi.org/10.1016/0021-9991(84)90128-1 Google Scholar CrossRef Search ADS   Pelties C., de la Puente J., Ampuero J.-P., Brietzke G.B., Käser M., 2012. Three-dimensional dynamic rupture simulation with a high-order discontinuous Galerkin method on unstructured tetrahedral meshes, J. geophys. Res. , 117, B02309, doi:10.1029/2011JB008857. https://doi.org/10.1029/2011JB008857 Peter D. et al.  , 2011. Forward and adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes, Geophys. J. Int. , 186( 2), 721– 739. https://doi.org/10.1111/j.1365-246X.2011.05044.x Google Scholar CrossRef Search ADS   Reed W.H., Hill T.R., 1973. Triangular mesh methods for the neutron transport equation, Los Alamos Report LA-UR-73-479. Robertsson J. O.A., Blanch J.O., Symes W.W., 1994. Viscoelastic finite-difference modeling, Geophysics , 59( 9), 1444– 1456. https://doi.org/10.1190/1.1443701 Google Scholar CrossRef Search ADS   Sánchez-Sesma F., Ramos-Martinez J., Campillo M., 1993. An indirect boundary element method applied to simulate the seismic response of alluvial valleys for incident P, S and Rayleigh waves, Earthq. Eng. Struct. Dyn. , 22( 4), 279– 295. https://doi.org/10.1002/eqe.4290220402 Google Scholar CrossRef Search ADS   Schroeder W., Martin K., Lorensen B., 2006. The Visualization Toolkit , Kitware. Schumacher F., Friederich W., Lamara S., 2016. A flexible, extendable, modular and computationally efficient approach to scattering-integral-based seismic full waveform inversion, Geophys. J. Int. , 204( 2), 1100– 1119. https://doi.org/10.1093/gji/ggv505 Google Scholar CrossRef Search ADS   Seriani G., Priolo E., 1994. Spectral element method for acoustic wave simulation in heterogeneous media, Finite Elem. Anal. Des. , 16(3–4), 337– 348. https://doi.org/10.1016/0168-874X(94)90076-0 Google Scholar CrossRef Search ADS   Tautges T.J., 2001. The generation of hexahedral meshes for assembly geometry: survey and progress, Int. J. Numer. Methods Eng. , 50( 12), 2617– 2642. https://doi.org/10.1002/nme.139 Google Scholar CrossRef Search ADS   Tessmer E., Kosloff D., 1994. 3-D elastic modeling with surface-topography by a Chebyshev spectral method, Geophysics , 59( 3), 464– 473. https://doi.org/10.1190/1.1443608 Google Scholar CrossRef Search ADS   Virieux J., 1986. P–SV wave propagation in heterogeneous media: velocity-stress finite-difference method, Geophysics , 51( 4), 889– 901. https://doi.org/10.1190/1.1442147 Google Scholar CrossRef Search ADS   Warburton T., 2006. An explicit construction of interpolation nodes on the simplex, J. Eng. Math. , 56( 3), 247– 262. https://doi.org/10.1007/s10665-006-9086-6 Google Scholar CrossRef Search ADS   Wilcox L.C., Stadler G., Burstedde C., Ghattas O., 2010. A high-order discontinuous Galerkin method for wave propagation through coupled elastic–acoustic media, J. Comput. Phys. , 229( 24), 9373– 9396. https://doi.org/10.1016/j.jcp.2010.09.008 Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Mar 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off