NUMERICAL MODELING OF SEISMIC WAVES BY DISCONTINUOUS SPECTRAL ELEMENT METHODS

. We present a comprehensive review of Discontinuous Galerkin Spectral Element (DGSE) methods on hybrid hexahedral/tetrahedral grids for the numerical modeling of the ground motion induced by large earthquakes. DGSE methods combine the ﬂexibility of discontinuous Galerkin meth-ods to patch together, through a domain decomposition paradigm, Spectral Element blocks where high-order polynomials are used for the space discretization. This approach allows local adaptivity on discretization parameters, thus improving the quality of the solution without aﬀecting the computational costs. The theoretical properties of the semidiscrete formulation are also revised, including well-posedness, stability and error estimates. A discussion on the dissipation, dispersion and stability properties of the fully-discrete (in space and time) formulation is also presented. Here space discretization is obtained based on employing the leap-frog time marching scheme. The capabilities of the present approach are demonstrated through a set of computations of realistic earthquake scenarios obtained using the code SPEED ( http://speed.mox.polimi.it ), an open-source code speciﬁcally designed for the numerical modeling of large-scale seismic events jointly developed at Politecnico di Milano by The Laboratory for Modeling and Scientiﬁc Computing MOX and by the Department of Civil and Environmental Engineering.


Introduction
In the last decades, the scientific research on elastic waves propagation problems modeled through the (visco)elastodynamics equation has experienced a constantly increasing interest in the mathematical, geophysical and engineering communities. In particular, the use of numerical simulations is a powerful tool to study the ground motion induced by earthquakes in regions threatened by seismic hazards. The distinguishing features of a numerical method designed for seismic wave simulations are: accuracy, geometric flexibility and scalability. To be accurate, numerical methods must keep dissipative and dispersive errors low. Geometric flexibility is required since the computational domain usually features complicated geometrical shapes as well as sharp discontinuities of mechanical properties.. Additionally, highly realistic earthquake models are typically characterized by domains whose dimension, ranging from hundreds to thousands square kilometers, is very large compared with the wavelengths of interest. This typically leads to a discrete problem featuring several millions of unknowns. As a consequence, parallel algorithms must be scalable in order to efficiently exploit high performance computers.
Historically, one of the first numerical methods employed for the solution of the elastodynamics equation was the finite difference (FD) discretization, see e.g. [8] and [53,84] for its application to real test cases. Dispersion errors of low order FDs can be mitigated by introducing a staggered grid approach for the velocity-stress formulation of the wave equation as in [77,111], or by using fourth-order FD schemes as in [22,35,75]. An overview of the stability and dispersion properties of three dimensional fourth-order FD scheme can be found in [83] for example. A major challenge for FD methods is related to the treatment of complicated domains and of free-surface or absorbing boundary conditions. Suitable techniques were developed in order to account for surface topography in [99]. We refer the reader to [82] for a comprehensive review of FD schemes applied to seismic waves propagation.
Another class of numerical methods that has been extensively used for the approximation of the elastodynamics equation is that of spectral methods. Introduced nearly fifty years ago for fluid dynamics problem, spectral schemes based on truncated Fourier series have been firstly applied to hyperbolic equations in [74] and to seismic waves propagation in [71]. To handle more general boundary conditions, the original version has been modified by using Chebychev or Legendre polynomials, leading to the so-called pseudo-spectral formulation (see for example [64,70,114]). An alternative approach to Spectral Methods consists in using iso-parametric approximations (e.g. quadratic mapping [27,66]) but this approach can be applied only to smooth surfaces. Spectral element (SE) methods have been designed to improve the flexibility of spectral methods, yet preserving their high accuracy. Firstly introduced in [86] for fluid dynamics problem, they are based on local polynomial approximations very similar to that provided by the p-version of the finite element method [20]. Successively, they have been extensively applied to the elastodynamics equation [50,51,67,90,92,102]. The success of spectral element discretization can be ascribed to their accuracy, flexibility and efficiency. Indeed, the use of high-order polynomial expansions not only provides exponential convergence to smooth solutions but also reduces drastically dispersion and dissipation errors. Moreover, the possibility to locally refine the grid allows to deal with complex geometries and sharp material discontinuities. Finally, the use of Gauss-Legendre-Lobatto interpolation points leads to a diagonal mass matrix and makes explicit time-integration schemes computationally efficient. Nowadays, thanks to the aforementioned properties, SE methods probably provide the most successful tool in computational seismology, in particular for large scale applications, see for example [29,32,68,69,109]. SE schemes are based on a fully continuous approximation which employs the same polynomial degree in all the cells. For this reason, in highly heterogeneous media or in regions where the geometry requires local refinement of the grid, the spread of the number of unknowns and the need of too restrictive conditions on the time step, can lead to high computational costs. Furthermore, SE methods are historically based on discretizations made by tensor product elements (i.e., quadrilateral elements in 2D and hexahedral elements in 3D). However, generating hexahedral grids for complex geometries may require a huge computational effort. Indeed, automatic mesh procedures with hexahedral elements are not able to accurately reproduce complex interfaces between different layers, especially when small angles are involved. The issues concerning the grid generation can be overtaken by using tetrahedral and/or hybrid meshes. Indeed, different (both commerical and open-source) softwares are available for the automatic generation of 3D unstructured tetrahedral grids, able to deal with very complex domains. To exploit this flexibility in the meshing process, high-order or spectral element methods have been extended to tetrahedral elements based on either nodal or modal approximations. Examples of nodal (or Lagrangian) basis functions can be found in [60,110,112], while a modal basis approximation has been originally proposed in [46]. The latter exploits a warped tensor product structure and has the very interesting property of being orthogonal with respect to the L 2 scalar product, but unfortunately it is not suitable for a global continuous approximations. This basis has been modified in [105] by introducing the so called boundary adapted basis in order to enforce the continuity between the elements. Differently from their hexahedral counterpart, both the nodal and the modal approach lead to mass matrices that are no longer diagonal. Therefore, when used with explicit time integrators, the inversion of a global mass matrix at each time-step has an impact on the computational efficiency. For this reason SE on tetrahedral grids encountered up to now a limited interest in the seismology community. An example of a triangular spectral element approximation based on nodal basis functions can be found in [80].
Another numerical method that in recent years has been extensively used for elastic waves propagation is the Discontinuous Galerkin (DG) method [43,61,95]. DG methods employ a discontinuous approximation of the solution based on standard finite elements formulation combined with suitable approximations of interface conditions between the elements, cf [43,61,95]. Thanks to their intrinsically local nature, DG methods are particularly well suited to handle highly heterogeneous media, e.g. strong contrasts in the soil properties, or in soil-structure interaction problems, where local refinements are needed to resolve the different spatial scales. Indeed, it is possible to tailor the local mesh size and the local polynomial approximation degree to the region of interest, according to the mechanical properties or the geometrical features of the computational domain. In this framework, the two main approaches developed in literature rely on the displacement formulation or on the stress-velocity formulation, respectively. In the displacement formulation the DG approximation is usually based on the interior penalty technique, both in its symmetric and non-symmetric version, see for example [41,98]. See also [9] for symmetric high-order DG methods for the elastodynamics equation in the displacement formulation. These schemes have been further improved to take into account viscoelastic effects in [96]. The stress-velocity formulation is usually combined with a flux approximation at the interfaces, similar to the one of the finite volume methods. Examples can be found in [47,65], where the flux term is treated by using an upwind method. It is worth mentioning that in these works the DG scheme is combined with an Arbitrary high order DERivatives (ADER) time integration scheme in order to obtain high-order accuracy in the space-time domain. The ADER-DG method has been successfully employed in several applications, from seismic wave propagation problems [32,48] to dynamic rupture simulation [87]. Moreover, in the ADER-DG framework, a discontinuous discretization on two dimensional non-conforming hybrid grids has been introduced in [59]. In order to reduce dissipative effects introduced in the aforementioned method, in [42,49] a DG method based on centered fluxes is proposed. A slightly modified version of this approach is considered in [81] to account for arbitrary heterogeneous media. DG approaches (both in its high-order or spectral element version) provide accurate solutions and are well suited for parallel implementation. All these methods are designed on general triangular and tetrahedral grids and employ different polynomial approximations strategies. For example, the ADER-DG scheme is based on a modal approximation that exploits the Dubiner basis, while the method presented in [81] employs a nodal basis built from the warp and blend interpolation nodes [61,112]. To mitigate the proliferation of unknowns that characterize DG methods, a Discontinuous Galerkin Spectral Element (DGSE) approach based on a domain decomposition paradigm has been introduced in [17]. More precisely, the discontinuities are imposed only at the interfaces between suitable non-conforming macro-regions, so that the flexibility of the DG methods is preserved without loosing the accuracy and efficiency of SE methods. Very recently, high-order DG methods for elastic wave propagation have been extended to arbitrarily-shaped polygonal/polyhedral grids [16], to further enhance the geometrical flexibility of the DG approach. One of the crucial aspects of numerical methods for waves propagation is their capacity to limit dispersive and dissipative effects affecting the discrete solution. DGSE methods feature very low dispersion and dissipation errors thanks to the their "SE engine". They require a small number of grid points per wavelength to retain the same level of accuracy [39]. The dispersion properties of high-order methods for the scalar wave equations have been analysed in [3] and in [4], where different kinds of DG approximations were considered. In the elastic case, the dispersive behaviour of spectral element methods has been analysed in [101] using a Rayleigh quotient approximation of the eigenvalue problem resulting from the dispersion analysis. DGSE approximations on quadrilateral grids have been investigated using a plane wave analysis in [41] and in [17]. A similar approach has been presented for spectral elements and DGSE approximations on triangular grids [76,78], where different sets of interpolating nodes have been compared, and in [15], where the authors have used the modal boundary adapted functions proposed in [106]. All these works deal with two-dimensional model problems and show that triangular spectral elements feature dispersion and dissipation properties similar to those of the standard tensor product spectral elements. The extension to three dimension is presented in [52], and a similar conclusion can be drawn.
The aim of this work is to present a comprehensive review on numerical modeling of seismic waves based on DGSE methods on hybrid hexahedral/tetrahedral grids. These methods combine the flexibility of discontinuous Galerkin methods to connect together, through a domain decomposition paradigm, Spectral Element blocks where high-order polynomials are used. DGSE methods are implemented in SPEED (SPectral Elements in Elastodynamics with Discontinuous Galerkinhttp://speed.mox.polimi.it), an open-source code aims at simulating large-scale seismic events in three-dimensional complex media: from far-field to near-field including soil-structure interaction effects. SPEED is jointly developed at Politecnico di Milano by The Laboratory for Modeling and Scientific Computing (MOX) of the Department of Mathematics and the Department of Civil and Environmental Engineering. More precisely, we will review the main advantages of this approach in terms of stability, accuracy, dispersion and dissipation properties, as well as computational costs.
The rest of this paper is organized as follows. In Section 1 we introduce the physical problem and governing equations. A particular emphasis is given to the review of the key ingredients needed to set up a highly realistic model for the numerical simulation of the ground motion induced by large earthquakes, namely, the modeling of the seismic source, absorbing boundaries and suitable damping terms. In Section 2 we present the Discontinuous Spectral Element method for the space discretization coupled with a leap-frog time marching scheme. In Section 3 we revise the theoretical results for the semi-discrete formulation and include the stability and the corresponding error estimates for the semi-discrete problem, whereas Section 4 is concerned with the dissipation, dispersion and stability analysis for the fully discrete scheme. In Section 5 we report some numerical computations of realistic earthquake scenarios in order to fully exploit the capabilities of the present approach. These results have been obtained using the open source code SPEED. Finally, some technical results needed for the theoretical analysis are contained in Appendix A.

Physical problem and governing equations
An earthquake is the result of a sudden release of energy in the outer layers of the Earth's crust. The strain energy slowly accumulated over the time is transformed into kinetic energy and radiated through the Earth's layers. By studying the main properties of motion, using suitable numerical techniques, it is possible characterize those features of ground shaking useful for a better determination of seismic actions for design. Seismic waves can be subdivided into body and surface waves. Body waves travel through the interior of the Earth and they are further split into pressure (P) and shear (S) waves. P-waves generate a ground motion that is aligned with the direction of the wave-field, they can travel through all media (solids and liquids) and are characterized by the highest wave velocities, typically ranging from 0.5 to 6 km/s. S-waves are transversal waves that induce a ground motion perpendicular to the wave propagation field and, differently from pressure waves, they can travel only through solid media. S-waves are slower than P-waves, with a wider range of propagation velocities, typically ranging from about 100m/s to 3.5Km/s. When body waves reach the Earth's surface, they generate surface waves, such as Love or Rayleigh waves, that induce a rolling and a circular ground motion, respectively. Surface waves are characterized by the slowest wave velocities and have a longer period of oscillation and a larger amplitude with respect to body waves.
In the following we introduce the mathematical model that describes the propagation of seismic wdmagniaves within the Earth's media. We consider a bounded domain Ω ⊂ R 3 (representing the portion of the ground where we study wave propagation), and assume that its boundary is decomposed into three disjoint portions Γ D , Γ N and Γ N R , where we impose the values of displacement, of tractions, and of the fictitious tractions introduced to avoid unphysical reflections, respectively.
Considering a temporal interval (0, T ], with T > 0, the dynamic equilibrium equation for a viscoelastic medium subject to an external force leads to the following formulation where ρ is the medium density, u = u(x, t) is the displacement field, σ(u) is the stress tensor and f = f (x, t) is a given external load (representing e.g. a seismic source). On the boundary we impose a rigidity condition on Γ D , a traction t = t(x, t) on Γ N and a fictitious traction t * = t * (x, t) on Γ N R . The exact expression of the term t * will be discussed later on. Finally, u 0 and v 0 are smooth initial values for the displacement and the velocity field, respectively. We assume that ρ is a uniformly bounded, strictly positive function. In (1) the parameter ζ ≥ 0 (usually referred to as damping factor ) is a decay factor with the dimension of the inverse of time and it is supposed to be piecewise constant in order to model media with sharp elastic impedance and damping behaviors. The term 2ρζu + ρζ 2 u can also be seen as an alternative or a complement to absorbing boundary conditions [72]. Usually, in engineering applications, if ζ > 0, the viscoelastic behavior of a material is expressed through the adimensional quality factor Q, defined as Q = Q 0 f /f 0 , where f 0 is a reference frequency and Q 0 = πf 0 /ζ.
For the stress tensor σ we use the following constitutive equation (Hooke's law): where (u) = (∇u + ∇ T u)/2 is the strain tensor, I is the identity tensor, λ = λ(x) ∈ L ∞ (Ω) and µ = µ(x) ∈ L ∞ (Ω) are the first and second Lamé elastic coefficients, respectively, and tr(·) is the trace operator. By introducing the fourth order Hooke's tensor D, relation (2) can also be written as with D a uniformly bounded and symmetric and positive definite tensor. We recall that the Lamé coefficients and the mass density are related to the compressional and shear wave velocities through the following relations respectively.

Modelling the seismic source
In this section we discuss the nature of the source term f in (1). In practical applications f can be a point-wise force acting on a point x 0 in the i th direction, i.e.
where e i is the unit vector of the i th Cartesian axis, δ(·) is the delta distribution, and f (·) is a function of time.
The expression of f (·) can be selected among different waveforms. A relevant example is the Ricker wavelet [94], defined as where f 0 is the wave amplitude, f p is the peak frequency of the signal and t 0 is a fixed reference time.
Another possible choice of the seismic input is a vertically incident plane wave. In particular, a uniform distribution of body forces along the plane z = z 0 of the form f (x, t) = f (t)e i δ(z − z 0 ) generates a displacement in the i th direction of the formū where H(·) is the Heaviside function and c (= c P or c S ) is the wave velocity, see [57]. Taking the derivative with respect to time of (7) and evaluating the result at z = z 0 we can express f (t) as Finally, one of the most important seismic input in earthquake simulation is the double-couple source force. Its mathematical representation is based on the seismic moment tensor m(x, t), defined in [5] as where n F and s F denote the fault normal and the rake vector along the fault, respectively. M 0 (x, t) describes the time history of the moment release at the point x and V is the elementary volume of the force. The equivalent body force distribution is finally obtained through the relation f (x, t) = ∇ · m(x, t), see [50]. We recall that the seismic moment M 0 , measured in N m, can be used to obtain the value of the moment magnitude M w of the earthquake as The moment magnitude scale, expressed through the adimensional value M w , is used by seismologists to classify the size of earthquakes and it is related to the strain energy released. A variation of one in the moment magnitude scale corresponds to a 2 times increase in terms of the energy released, so that an earthquake of magnitude M w ≈ 8 will release 2 times the energy of an earthquake of magnitude M w ≈ 7.

Modelling absorbing boundaries
One crucial aspect to be taken into account for the accurate simulation of (visco)elastic waves propagation in unbounded domains is how to properly choose absorbing boundary conditions at the artificial boundaries of the numerical model. A possible approach consists in modelling the absorbing boundary layers by introducing a fictitious traction term on Γ N R , consisting of a linear combination of displacement space and time derivatives. Here, we consider the local P3 paraxial condition presented in [107], which is sufficiently accurate if c P /c S ≤ 2, as in the application under consideration. Let n = [n x , n y , n z ] T be the unit normal vector to Γ N R and let τ 1 = [τ 1,x , τ 1,y , τ 1,z ] T and τ 2 = [τ 2,x , τ 2,y , τ 2,z ] T be a couple of mutually orthogonal unit vectors that lie on the plane tangent to the boundary, the P3 paraxial absorbing conditions read as .
The stress σ * = σ * (u) defined on the absorbing boundary in the local coordinate system (τ 1 , τ 2 , n) has the following expression By inserting the above espression in (11), the traction term Finally, the expression of t * in the global coordinate system can be recovered by writing the normal and the tangential derivatives as and then projecting the resulting vector on the global coordinate system by Details on the practical choice of the vectors (τ 1 , τ 2 ) can be found in [28].

Numerical discretization
In this section we describe the numerical approximation of the (weak formulation) of (1) through a Discontinuous Spectral Element method coupled with a leap-frog time marching scheme.
For an open, bounded, polygonal domain D ⊂ R 3 and for a non-negative integer s we denote by H s (D) the L 2 Sobolev space of order s and by · s,D and | · | s,D the corresponding norm and semi-norm. In addition, for s = 0, we write L 2 (D) in place of H 0 (D). We set H s (D) = [H s (D)] 3 and H s (D) = [H s (D)] 3×3 , and denote by (·, ·) D and ·, · ∂D the corresponding L 2 and L 2 inner products. Finally, we introduce the following space for time dependent functions defined on the interval (0, T ) cf. [1]. For vector-and tensor-valued functions the above space is defined analogosuly. The space C 0 (0, T ; H s (D)) and C 1 (0, T ; H s (D)) are defined analogously.
supplemented with the initial conditions If Γ N R = ∅ and ζ = 0, the above problem is well-posed and its unique solution

Partitions and trace operators
We consider a (not necessarily conforming) decomposition T Ω of Ω into L nonoverlapping polyhedral subdomains Ω , i.e., Ω = ∪ Ω , Ω ∩ Ω = ∅ for = . On each Ω , we built a conforming, quasi-uniform computational mesh T h of granularity h > 0 made by open disjoint elements K j , and suppose that each K j ∈ Ω is the affine image through the map F j : K −→ K j of either the unit reference hexahedron or the unit reference tetrahedron K. Given two adjacent regions Ω ± , we define an interior face F as the non-empty interior of ∂K + ∩ ∂K − , for some K ± ∈ T h ± , K ± ⊂ Ω ± , and collect all the interior faces in the set F I h . Moreover, we define F D h , F N h and F N R h as the sets of all boundary faces where displacement, traction or fictitious tractions are imposed, respectively. Implicit in this definition is the assumption that each boundary face can belong to exactly one of the sets Finally, we collect all the boundary faces in the set F b h . To carry out the analysis, we suppose that the following mesh assumption holds, see [55,88], (see also [44,45] for the case of highly discontinuous coefficients).
Assumption 2.1. Mesh assumption. For any element K ∈ T h and for any face F ⊂ ∂K, it holds h K h F . We refer to [25], for example, for the weakening of the above assumption and to [11,24] for for the use of polyhedral-shaped elements. Let K ± ∈ T h ± , K ± ⊂ Ω ± be two elements sharing a face F ∈ F I h , and let n ± be the unit normal vectors to F pointing outward to K ± , respectively. For (regular enough) vector and tensor-valued functions v and τ , we denote by v ± and τ ± the traces of v and τ on F , taken within the interior of K ± , respectively, and set where

Discontinuous Galerkin Spectral Element approximation
To each subdomain Ω we assign a nonnegative integer N , and introduce the finite dimensional space for where we have used the short-hand notation where, as before, if there is only one subdomain, F I h = ∅ and the formulation (16) correspond to the classical Spectral Element method. In this case, we will denote the corresponding discrete space as V CG , to highlight that the corresponding semi-discrete solution is globally continuous. On the other hand if L > 1 the discrete solution is piecewise discontinuous across macroelements and (weak) continuity is enforced based on employing, at a subdomain level, the symmetric interior penalty DG (SIPG) method [17]. Notice that here the essential boundary conditions are enforced strongly. In the context of a "fully" DG approach, i.e., the DG paradigm is employed elementwise, relevant applications of the SIPG method to elastic wave propagation problems can be found in [41,97], while in [98] the same problem is solved based on employing the non-symmetric version of the interior penalty method (NIPG) [113]. Other DG methods based on the flux approximation at element interfaces have been considered for the elastodynamics problem in [42,47,81] for the stress-velocity formulation and in [10] for both the displacement and the stress-displacement formulation. We refer to [10] for a unified analysis of the h-version of the method and to [13,15] for the hp-version of the method and its analysis.

Fully discrete formulation
In this section we address the issue of integrating in time the semi-discrete formulation (16). By fixing a basis for the discrete space V DG , the semi-discrete algebraic formulation of problem (16), reads as supplemented with suitable initial conditions U(0) = U 0 andU(0) = V 0 . Here, denoting by N dof be total number of degrees of freedom, the vector U = U(t) ∈ R N dof contains, for any time t, the expansion coefficients of the semi-discrete solution u h (t) ∈ V DG in the chosen set of basis functions. Analogously, M, M 2 , A, M 3 are the matrices representations of the bilinear forms , respectively, cf. (16). When absorbing boundary conditions are included in the model, the matrices S and R take into account the boundary term t * (u), v F N R h , otherwise they are identically equal to the null matrix.
Finally, F is the vector representation of the linear functional F h (·) cf. (16).
For the time integration of the system of second order ordinary differential equations (19), we employ the leap-frog method [91], that it widely employed time integration scheme for the numerical simulation of elastic waves propagation, see for example [21,31,67,81]. Other time integration techniques based on high-order time marching scheme can be employ to discretize in time the (visco)elasdodynamics equation, such as Runge-Kutta methods [64], the ADER-DG method [65] or the space-time DG discretization [12]. With this aim, we subdivide the time interval (0, T ] into N T subintervals of amplitude ∆t = T /N T and we denote by U i ≈ U(t i ) the approximation of U at time t i = i∆t, i = 1, 2, . . . , N T . Whenever ζ = 0, i.e. no dumping is present and Γ N R = ∅, then M 2 , M 3 and R are null and the leap-frog method reads as We notice that (20) involves a linear system with the matrix M to be solved at each time step. The choice the basis functions spanning the space V DG strongly influences the structure of the mass matrix M and, therefore, the computational costs related to the solution of the linear system. Furthermore, the leap-frog method being an explicit second order accurate scheme, to ensure its numerical stability a Courant -Friedrich -Levy (CFL) condition has to be satisfied (see [91]). In this general case, i.e. ζ = 0 or Γ N R = ∅, we set K = M 2 + S and Q = A + M 3 + R. Then, the leap-frog method consists in computing 3. Analysis of the semi-discrete formulation: stability and error estimates In this section we review the results on stability, convergence, and corresponding error estimates of the semidiscrete problem (16). For the sake of presentation we assume null tractions, i.e., t = 0, and lack of absorbing boundaries, i.e., Γ N R = ∅. Moreover, we assume that each macro region Ω , = 1, 2, . . . , consists of a single element, either a hexahedron or a tetrahedron, having size h and characterized by a polynomial of degree N . Notice that in this case the subdomain partition T Ω coincide with the mesh T h , and the skeleton F I h includes all the non-empty intersections between neighboring elements Ω ± . In the following, to avoid the proliferation of constants, we use the notation a b to represent the inequality a ≤ Cb, where C is a given positive constant, independent of the discretization parameters, but that can depend on the material properties ρ, λ, µ.
Before proceeding we properly define facewise the stabilization function η ∈ L ∞ (F I h ) in (18) as where α is a positive constant to be properly chosen. For a piecewise constant tensor D, {D} H is defined as

Well-posedness
Denoting by H s (T Ω ) the space H s (Ω ) of piecewise H s vector-valued functions, for any real s ≥ 0, we recall some useful inequalities.
We refer, for example, to [2,23,26,100,104] for further details and for the proof.
We next introduce the following (discretization parameters dependent) norms with the convention that w where the second estimate holds provided that the parameter α appearing in the definition of the stabilization function cf. (22) is chosen sufficiently large. Moreover,

Stability
In this section we present a stability analysis of the semi-discrete formulation (16) in the following (mesh dependent) energy norm cf. [10,13] for the case ζ = 0. For t = 0, the above definition modifies as where u 0,h , v 0,h ∈ V DG are suitable approximations of the initial data u 0 and v 0 , respectively, cf. (1). Differently than what done in [96][97][98], the present analysis holds without including an extra term that penalizes the jump of the time derivative of the displacement.  (22), u h (t) ∈ V DG be the solution of problem (16). If f ∈ L 2 (0, T ; L 2 (Ω)), then Proof. Our proof makes use of some technical results that are proved in Appendix A. By taking v =u h in (16), we obtain where we have used the definition of the energy norm (25). Integrating the above inequality over the time interval (0, t) we have Since where the first bound holds provided that the stability parameter α appearing in the definition of the penalization function (22) is chosen large enough. Then The Cauchy Schwarz inequality leads to The thesis follows based on employing the Gronwall's lemma A.1.

Semi-discrete error estimates
In this section we derive a priori error estimates for the semi-discrete problem (16) in the energy norm (25). We start by recalling the following interpolation estimates (see [19]) The above result stands at the basis for the following interpolation bounds.
Moreover, if in addition v,v ∈ H s (Ω ), s ≥ 2, = 1, . . . , L, then Proof. We first show (33). With this aim, use Lemma 3.4 to bound each contribution appearing in the definition of the DG norm |||·||| DG This yields To show (34), we recall the definition (25) of energy norm · E and employ again the interpolation estimates of Lemma 3.4 to get Summing up the above contributions we obtain and the proof is complete.
These bounds are optimal in h and suboptimal in N by a factor N 1/2 . We refer to [62,88] for analogous bounds for stationary (scalar) second order elliptic problems and to [13] for the case of elastic wave propagation problems. Optimal error estimates with respect to the polynomial approximation degree N can be shown either using the projector of [56] (provided the solution belongs to a suitable augmented space), or whenever a continuous interpolant can be built; cf. [108].
Proof. (Proof of Theorem 3.6.) The discrete formulation (16) is strongly consistent, that is the exact solution u of (1) satisfies equation (16) for any time t ∈ (0, T ] Subtracting (16) from the above identity and setting e = u − u h , we obtain the error equation We next decompose the error as e = e I − e h , with e I = u − u I and e h = u h − u I , u I ∈ V DG being the interpolant defined as in Lemma 3.5, and rewrite the above identity as By taking as test function v =ė h , we get Observing that ρ 1/2 ζ 1/2ė h 0,Ω ≥ 0, integrating in time between 0 and t, and using that e h (0) = u h (0)−u I (0) = 0 Employing Lemma A.3 we obtain We next estimate the last two terms on the right hand side. For the first one, the integration by parts formula (78) together with the observation that e h (0) = u h (0) − u I (0) = 0, and Lemma 3.2 lead to where in the last step we have used the definition of the energy norm (25). Analogously, for the second term By substituting the above two bounds in (45) we get From the Young's inequality, it follows and therefore, we can therefore chose = 4C being C the hidden constant in (51) and obtain from (51)  we obtain The Therefore we obtain and the proof is complete.

Dispersion, dissipation and numerical stability analysis
The aim of this section is to present an analysis of the dispersion, dissipation and stability properties for the fully discrete numerical scheme.
If the numerical wave shows a phase leg with respect to the physical one we have a dispersion effect, whereas if we observe a decrease in amplitude we are in presence of dissipation (cf. Figure 1, left). A traveling wave in a three dimensional elastic medium can be decomposed into: i) a pressure wave (P-wave) with velocity c P inducing a displacement in the same direction of the propagating wave; ii) a shear wave (S-wave) with velocity c S inducing a displacement in a direction transversal to the propagating wave. Additionaly, S-waves can be further decomposed into a vertical component (SV-wave), inducing a motion on a plane perpendicular to the wave direction, and a horizontal component (SH-wave), inducing a transversal motion on a horizontal plane containing the wave direction, see [5,103] and Figure 1 (right). The Von Neumann analysis [33,63] investigates the properties of numerical solutions approximating plane waves of the form propagating in an unbounded domain and looks for quantitative estimates of the dispersion and the dissipation errors. Here, a = [a 1 , a 2 , a 3 ] T ∈ R 3 represents the amplitude of the wave, ω the angular frequency, and k = 2π/L(cos θ cos φ, sin θ cos φ, sin φ) the wave vector, being L the wavelength, and θ and φ the angles between the propagating direction and the coordinate axes. The physical wave can be recovered by taking the real part of (59).
To comply with unboundedness, we consider problem (1) posed over a reference domain E C and impose periodic boundary conditions on its boundary, see Figure 2. In the case of a hexahedral tessellation, the smallest periodic grid is made by a single cube, i.e. E C = (−1, 1) 3 whereas for tetrahedral grids it is composed by six tetrahedra, Figure 2. Given that we are considering a DG spectral element method, the interelement jump and average contributions are assembled at the interfaces between E C and its neighbors (periodic reference pattern), see Figure 2. Assuming, without loss of generality, no seismic source (i.e. f ≡ 0), no dumping factor (i.e. ζ = 0) and no reflecting boundaries (i.e. Γ N R = ∅), the semi-discrete problem (19) becomes where U 0 = ae i(k·x) and V 0 = −iωae i(k·x) . Follow the approach of [15,17,41,52,76,78] we impose periodic boundary conditions by introducing a suitable projection matrix Π and we obtain from (60) where M = Π T M 1 Π and A = Π T AΠ. We next consider the fully discrete formulation based on employing the leap-frog time integration scheme (20) to (61); the analogous results for the semidiscrete case can be found in [52]. Following [6], we substitute (59) into (61) and we obtain The above system can be rewritten as Now, using the following relation we obtain a generalized eigenvalue problem of the form where the eigenvalues Λ are related to the angular frequency ω at which the wave travels in the grid through the relation We will use this after solving the eigenvalue problem in order to derive the grid-dispersion relations as it will be shown later on. We remark that for three dimensional seismic modeling only three eigenvalues in (64) have a physical meaning as they are related to shear and pressure waves, cf. [39,41] for the bidimensional case. All the other eigenvalues correspond to nonphysical modes, see e.g. [33] for the one dimensional case.

Dispersion errors
The relative dispersion errors are given by where c P and c S are the P-and S-wave velocities (see (4)) whereas c P,h and c S,h are the corresponding numerical wave velocities whose expression is given by where δ = h/(N L) is the sampling ratio, i.e., δ −1 is the number of interpolation points per wavelength, h is the mesh size and ω P,h and ω S,h are the numerical angular frequencies computed through (65) for the P-wave and the S-wave, respectively. In practice we first solve numerically (64) to obtain the eigenvalues in (65) that represent the best approximations of the angular frequencies of the travelling waves. For a 3D plane wave we can distinguish between the frequency ω P , which is related to the longitudinal displacement, and the two frequencies ω SV = ω SH , which are related to the transversal displacement of vertically polarized (SV) and horizontally polarized (SH) shear waves, respectively. In general, as pointed out before, the number of eigenvalues obtained through (64) exceeds the number of physical modes. Therefore, we identify the numerical eigenvalues Λ P and Λ S , corresponding to the physical frequencies, by computing the numerical velocities obtained for each eigenvalue and comparing them to the real values of c P and c S , respectively. We remark that the computed   eigenvalues approximating ω SV and ω SH are not exactly the same but their difference is negligible, cf. [101,115]. In the following we will select Λ S as that eigenvalue, between the two physically relevant ones, that leads to the worst approximation of c S .
We next present some numerical results aiming at measuring the dispersion error, versus the polynomial approximation degree N , the sampling ratio δ, and the time step ∆t. For the sake of brevity we do not discuss here the dependence of the dispersion error on the angles of incidence, and refer to [52] for more details. We first address the behavior of the dispersion error by varying the polynomial degree N . In particular, we compute the dispersion errors versus N for different time steps ∆t, fixing the values of velocities, sampling ratio and incident angles as c P = 1, c S = 0.5, δ = 0.2 and φ = θ = π/4, respectively. In Figure 3 we observe that, as ∆t goes to zero, e P and e S decay exponentially fast until the term ∆t 2 becomes dominant, cf. [15,52]. As a comparison, Figure 3 also report the corresponding results for the semi-discrete approximation, i.e. using a sufficiently small time step in such a way that exact time integration can be assumed. We also notice that the same level of accuracy is obtained on both hexahedral and tetrahedral grids for a polynomial degree N ≥ 5. We have repeated the same set of numerical experiments varying the sampling ration δ, with N = 1. The results are reported in Figure 4. We can clearly observe that for δ ≤ 0.1, i.e., ten points per wavelength, the same level of accuracy is attained on both hexahedral and tetrahedral grids. Regarding the dissipation of the fully-discrete scheme, the considerations made in the semi-discrete case are still valid.

Dissipation errors
We next move to the dissipation analysis that is carried out by studying the amplitude of the numerical displacement. By considering as the exact solution of (60) the unitary amplitude plane wave, we can express its amplitude as |e i(k·x−ωt) | = e tIm(ω) . Since the physical wave satisfies Im(ω) = 0, its amplitude is equal to 1 for all times t. On the contrary, the numerical wave will have in general Im(ω h ) = 0. Then, we say that the scheme is non dissipative if Im(ω h ) = 0, whereas it is dissipative if Im(ω h ) < 0. In the generalized eigenvalues problems (64) the mass and the stiffness matrices are symmetric and positive definite. Therefore, the computed eigenvalue are all real, leading to a scheme that does not suffer from dissipation errors.

Numerical stability
In this section we briefly analyze the stability properties of the fully discrete scheme (21), focusing as before on the case where no dumping is present, i.e. ζ = 0, and no reflecting boundaries are considered, i.e. Γ N R = ∅. According to the Courant, Friedrichs and Lewy (CFL) condition, the leap-frog scheme is stable provided that the time step ∆t satisfies where h is the granularity of the computational grid and c P is the velocity of pressure waves. The constant C CF L depends both on the material properties ρ, λ, µ and on the polynomial degree N . We consider a generalized eigenvalue problem of the form (64) but rescaled on the reference element K where Λ = (h/∆t) 2 sin 2 (ω h ∆t/2) and introduce the stabilization parameter q = c P ∆t/h. The following inequality holds or, equivalently, As shown in [40], Λ depends on the wave vector k and therefore on the values of the incident angles φ and θ. Thus, (71) can be formulated as where Λ max is the maximum eigenvalue of (69), taken with respect the values of φ and θ. The constant c * depends on the Lamé parameters λ and µ as well as on the penalty paramenter α appearing in the definition of the stabilization function (22), more precisely it is proportional to α −1/2 , see [14].
We next present some numerical experiments to give a quantitative estimate of the CFL bound expressed in term of the stability parameter q. The material properties considered are the same as in the previous section and the stability parameter α = 10. In Table 1 we report the computed value of q as a function on N based on employing both the reference hexahedron and tetrahedron. We can observe that, as expected, q grows proportionally to N 2 . This result is in agreement with the analogous ones obtained in [15,17,40,52]. In addition we notice that tetrahedral elements are subjected to a more restrictive stability condition than hexahedral ones. In particular, in the case of a continuous (resp. discontinuous) approximation on a tetrahedral grid, the stability parameter q is around the 60% (resp. 80%) of the corresponding one computed using a hexahedral mesh.  Table 1. Computed stability parameter q as a function of N = 2, .., 8 on both the hexahedral and tetrahedral reference elements.

Layer over half space benchmark
In this example we apply our method to a more challenging test case, known in the literature as the LOH (Layer Over a Half-space) test. This problem has been proposed in [38] and it has been used as reference benchmark for different numerical codes for the elastodynamics [37,47]. The domain Ω consists of a rectangle of dimension (−15, 15) × (−15, 15) × (0, 17) km, with a 1 km layer (L) on the top surface (see Figure 5). The material properties of the layer and of the halfspace (HS) are summarized in Table 2. The numerical results presented in the following have been obtained with the following choices of the discretization parameters. In the half-space we employ a hexahedral grid with a fix size of h ≈ 600 m, while in the layer we consider a hybrid discretization made by tetrahedral elements in the positive quadrant and by hexahedral elements of fixed size h ≈ 300 m in the remaining portion, see  We also remark that the grid of each portion of the domain has been built independently from the others, leading to a grid characterized by non conforming interfaces on the skeleton.
To compare the results obtained with the reference solution we report the time histories of the velocity field measured at a receiver located at the point (6,8,0) km. For each component of the velocity we also compute the relative seismogram misfit E through the formula where ns is the number of time samples in the seismogram,u h (t i ) is the numerical velocity at the time t i anḋ u(t i ) is the value of a quasi-analytical solution, which can be found in [18,38]. In Figures 7, 8 and 9 we report the radial, transversal and vertical components of the velocity field recorded at the point (6,8,0) km. From the velocity histories we observe that the setting for the case A provides the worst approximation. The numerical solution, indeed, presents spurious oscillations, in particular in the transversal component. This oscillatory effect is probably due to dispersion effects caused by the low order approximation. In both the cases B and C the oscillations are mitigated and the value of the seismic misfit is below the 10% and 4%, respectively. The small oscillations that can be appreciated after 6.5 s are probably due to reflecting waves generated by the artificial boundaries. From the numerical results, we can conclude that our approach seems to have at least the same accuracy of other finite difference and finite element numerical approaches [37,48,79].
Finally, in Table 3 we report the number of elements and of degrees of freedom used in the tetrahedral portion of the layer, together with the computational times needed for the set-up, the assembly of the mass matrix and to compute the solution at each time step. The results show that the setting used for test case B (h ≈ 200 m, N = 3) gives the best compromise in term of accuracy and computational times.   Table 3. LOH test case. Numerical settings, computational times and average seismic misfit for the three mesh configuration.

Comparison with a Spectral Element discretization
In this section we want compare the overall computation cost of a DGSE discretization with respect to a SE one. In this respect, we consider for the LOH problem described in Section 5.1 two different computational models: the first one consists of a hexahedral conforming mesh, while the second one of a hexahedral non conforming mesh, cf. Figure 10. In Table 4 we report the discretization parameters employed for both computational models.
For completeness, in Figure 10 Figure 11 we can infer that the average relative error with respect to the reference solution is 1% for the conforming SE model and around 1.3% for the DGSE non-conforming model. It is important to underline that the non conforming models have a sensibly lower number of elements, nonetheless the accuracy of the solution is almost preserved (see Figure 11). Further improvements of the results can be achieved, for instance, choosing higher order polynomial degrees. Finally, in Table 5 we compare the time-to-solution associated to the SE and the DGSE discretizations. Both simulations ran on 512 cores on the Fermi cluster located at Cineca (https://www.cineca.it/en/content/fermi-bgq). Comparing the results it is possible to see that the DGSE discretization required almost the same amount of time per step respect to the SE one. For a thorough  comparison between the SE and DGSE methods as well as scalability analysis of the DGSE approximation we refer the reader to [36,79].

Earthquake scenarios for the area around Beijing
The accurate evaluation of seismic risk scenario in large urban areas is based on suitable approaches for earthquake ground motion prediction. However, the empirical tools that are most often used for this purpose neglect the specific tectonic and geological conditions in which the urban area lies, thus producing in many cases significant underestimations of the earthquake ground motion intensity. An overview of applications of physics-based numerical modelling approaches to large urban areas was presented in [85]. Specifically, SPEED has been applied so far to several case studies, including Christchurch and Wellington (New Zealand), Santiago (Chile), Istanbul (Turkey), and several urban areas in Italy. In this section, we aim at summarizing results obtained from a large set of earthquake scenario simulations in Beijing (China).
Situated on a sedimentary basin, with its more than 20 million inhabitants, Beijing is one of the many megacities around the world highly exposed to the seismic threat. As a matter of fact, in the past thousand   Table 6. Geometric parameters of the Shunyi-Qianmen-Liangxiang fault. Fault Origin is defined as the point of the fault at zero strike and zero dip.
years, destructive earthquakes have occurred in the area around Beijing, with magnitude varying from M w 6 to M w 6.5, cf. [58]. For these reasons, the characterization of strong ground motion plays a crucial role for the risk assessment studies in this large urbanized area.

Setup of the 3D computational model
The 3D numerical model for the Beijing area was set up exploiting: i) the topography model; ii) depth of sedimentary base; iii) the kinematic seismic fault model; iv) and the 3D velocity profiles. For the elevation model, a freely-available digital elevation dataset of CGIAR-CSI for the Beijing region has been extracted and downloaded from the website http://www.cgiar-csi.org (with a precision of roughly 90 × 90 m, for east-west and north-south directions around Beijing city), while the depth of sedimentary base model has been derived from the digitalization of the map proposed in [54], as shown in Figure 12.
The seismic fault considered is named Shunyi-Qianmen-Liangxiang fault, lying just through the downtown of Beijing (Figure 12, cf. also Figure 13 left.). This is a quasi-normal segmented fault (with dip angle of about 80 • ), consisting of three main segments with different strike angles. The total length of the fault is about 90 km, capable of producing up to M w 7.3 events. The geometric parameters of the Shunyi-Qianmen-Liangxiang fault, as implemented in the numerical model, are reported in Table 6.  4  3500  2100  2200  200  3  4 -12  6000  3400  2760  800  3  12 -30  6200  3500  2810  900  Table 7. Horizontally stratified crustal model, cf. [54].
In order to define the 3D soil model, the V s,30 map provided by MunichRe company, cf. [7], and the depth of sediment base obtained previously, see Figure 13, were used. In particular, in the first layer (0 to 2 km depth), we defined different velocity profiles (in m/s) as follows where z top and z sed represent the projection of a generic point with coordinate z into the surface and the sediment base, respectively. Similarly, for the soil density (kg/m 3 ) we consider the following profiles 30 >= 600 m/s ρ = 1530 + 5 |z − z top |, for V s,30 < 600 m/s and z >= z sed , ρ = 1800 + 5 |z − z top |, for V s,30 < 600 m/s and z < z sed .
(75) Figure 13 (right) shows the c s model obtained for the first layer. The properties of the underlying bedrock layers (depth > 2 km) have been selected in agreement with [54]. The quality factor Q is estimated directly by the c s values and is assumed to be proportional to frequency, see Section 1, for the target value Q = c s /10 to be obtained at f = 1 Hz. The 3D velocity model is summarized in Table 7.
The computational domain, spanning an area of 70 × 70 km 2 down to 30 km depth (see Figure 14), was built by combining all the information above. Considering a rule of thumb of 5 grid points per minimum   Table 8. Simulated scenarios grouped for magnitude.
wavelength for non-dispersive wave propagation in heterogeneous media by the SE approach (see [79]), and considering a maximum frequency f max = 1.5 Hz, the model consists of 859,677 hexahedral elements, resulting in approximately 160 million degrees of freedom, using a fourth order polynomial approximation degree. The conforming mesh has a size varying from a minimum of 150 m, on the top surface, up to 600 m at 4 km depth and reaching 1800 m in the underlying layers.

Summary of simulated scenarios
A total of 30 scenarios were simulated by varying the magnitude, from 6.5 up to 7.3, the kinematic slip distribution, the hypocenter location and the location of the rupture area. The simulations were performed on the Marconi clusterì at CINECA, Italy (http://www.cineca.it/en/content/marconi). Each simulation takes around 12 hours on 512 cores. A summary of the seismic scenarios, grouped according to the magnitude, is given in Table 8.
To automatically construct N physically constrained slip distributions for a given fault and a given earthquake magnitude, a pre-processing tool has been devised taking into account joint probability distributions of the main kinematic parameters. This ensures that the resulting scenario variability will not be affected by systematic bias in the input parameters. In particular, we considered the kinematic source rupture generator proposed by Crempien and Archuleta (2015) (see [34]). We refer the reader to the relevant publication and to the SCEC Broadband Platform for the documentation of the code. Note that for each scenario, the rupture velocity follows the built-in scheme proposed in [34], and the source time function is a simplified smoothed Heaviside function. A time step equal to 0.001 s has been chosen and a total observation time T = 60 s has been considered. In order to model the non-linear soil behavior of the soft soil deposits (V s,30 <= 400 m/s and z top <= z <= z top −300 m) a simple Non-Linear Elastic (NLE) soil model is implemented as a generalization to 3D load conditions of the classical modulus reduction (µ − γ) and damping (ζ − γ) curves used within 1D linear-equivalent approaches (e.g. [73]), where µ, ζ and γ are the shear modulus, damping ratio and 1D shear strain, respectively. Namely, to extend those curves to the 3D case, a scalar measure of shear strain amplitude is considered as follows: where ε I , ε II and ε III are the principal values of the strain tensor. Once the value of γ max is calculated at the generic position x and generic instant of time t, its effective value (set as 0.6γ max ) is introduced in the µ − γ and ζ − γ curves to update the corresponding values for the following time step, cf. Figure 14 (right). Therefore, unlike the classical linear-equivalent approach, in this non-linear elastic approach the initial values of damping and stiffness are recovered at the end of the excitation.
In Figure 15 (left column) we report the maps of the Peak Ground Velocity (PGV) obtained for the average scenario with magnitude M w 6.5 (top), M w 6.9 (middle) and M w 7.3 (bottom). We remark that: i) these maps have been recovered by averaging the complete set of scenarios sharing the same magnitude, ii) only important urbanized areas are reported. In the same figure (central column) a comparison with the ground motion prediction equation (GMPE) proposed by [30], here denoted by CAEA15, is also reported. We note that the distance metric used here is R rupt (closest distance to the fault rupture), but analogous results can be achieve with other metrics. The scattered plot have then been fitted by a nonlinear regression model and the resulting curve are reported in Figure 15 (right column). The numerical results obtained by SPEED are substantially in agreement with the proposed GMPE. However, synthetic scenarios produce higher PG values for the velocity field in the proximity of the fault rupture. This has already been noted in many recent works, see [85] and references therein, and may have a great impact on seismic hazard estimates, as standard GMPEs cannot account for such effects. Finally, we also report in Figure 16 some snapshots of the peak ground velocity wave field for a target scenario with magnitude M w 6.5. It is possible to notice that the wave field propagates in the direction south-west to north-east and that higher values of PGV are obtained closed to the projection of the fault rupture on the surface.

A. Appendix
In this section we provide some technical results needed for the proof of Theorem 3.3. The first key ingredient is the following Gronwall lemma, see [89] for example.
Based on standard DG techniques, the following preliminary bound holds.
where α is the stability parameter appearing in the definition of the penalization function (22).
Proof. Let F ∈ F I h be an interior face shared by two neighboring elements Ω ± . Using the definition of the average operator it clearly holds Recalling that σ(v ± ) = D ε(v ± ), it is therefore enough to show since the thesis then follows by summing over all the interior faces F ∈ F I h . Using the definition of the penalization function (22), the local bounded variation property, together with the trace-inverse inequality (23b) we have where the hidden constant depends on the material properties through the quantity D ∞ / {D} H .
The second bounds holds provided the stability parameter α appearing in the definition of the penalization function (22) is chosen large enough.
Proof. The first bound follows immediately based on employing the Cauchy-Schwarz inequality, Lemma A.2, and the definition of the energy norm (25). That is, To prove the second bound, it is sufficient to show that Indeed, using again the definition of the energy norm (25) and the above estimate we immediately have Therefore, we next show that (85) holds provided that the stability parameter α appearing in the definition of the penalization function (22) is chosen large enough. To this aim, we observe that by using the Young inequality we have, for any positive number , where η is the stabilization function defined in (22). Therefore, from the definition of the DG norm (24) it follows