A stable discontinuous Galerkin method for linear elastodynamics in 3D geometrically complex media using physics based numerical fluxes

TL;DR

This paper introduces a stable, high-order discontinuous Galerkin method for 3D elastic wave simulations in complex geometries using physics-based fluxes, enabling accurate and adaptive seismic modeling.

Contribution

It develops a physics-based, energy-stable DG scheme on curvilinear meshes for 3D elastodynamics, incorporating adaptive mesh refinement and complex boundary handling.

Findings

Demonstrates stability and high-order accuracy in heterogeneous media.

Validates the method with 3D wave propagation experiments.

Shows potential for regional seismic simulations with complex topography.

Abstract

High order accurate and explicit time-stable solvers are well suited for hyperbolic wave propagation problems. As a result of the complexities of real geometries, internal interfaces and nonlinear boundary and interface conditions, discontinuities and sharp wave fronts may become fundamental features of the solution. Thus, geometrically flexible and adaptive numerical algorithms are critical for high fidelity and efficient simulations of wave phenomena in many applications. Adaptive curvilinear meshes hold promise to minimise the effort to represent complicated geometries or heterogeneous material data avoiding the bottleneck of feature-preserving meshing. To enable the design of stable DG methods on three space dimensional (3D) curvilinear elements we construct a structure preserving anti-symmetric coordinate transformation motivated by the underlying physics. Using a physics-based numerical penalty-flux, we develop a 3D provably energy-stable discontinuous Galerkin finite element approximation of the elastic wave equation in geometrically complex and heterogenous media. By construction, our numerical flux is upwind and yields a discrete energy estimate analogous to the continuous energy estimate. The ability to treat conforming and non-conforming curvilinear elements allows for flexible adaptive mesh refinement strategies. The numerical scheme has been implemented in ExaHyPE, a simulation engine for parallel dynamically adaptive simulations of wave problems on adaptive Cartesian meshes. We present 3D numerical experiments of wave propagation in heterogeneous isotropic and anisotropic elastic solids demonstrating stability and high order accuracy. We demonstrate the potential of our approach for computational seismology in a regional wave propagation scenario in a geologically constrained 3D model including the geometrically complex free-surface topography of Mount Zugspitze, Germany.

Tables

Table 1. Table 1: Receiver positions of the 3D HHS and LOH1 problems relative to the epicenter.

Receiver	1	2	3	4	5	6	7	8	9
y[km]	0	0	0	0.490	3.919	7.348	0.577	4.612	8.647
z[km]	0.693	5.542	10.392	0.490	3.919	7.348	0.384	3.075	5.764

Table 2. Table 2: Accuracy levels

Accuracy level	Envelop Misfit	Phase Misfit
A	$\le 5\%$	$\le 5\%$
B	$\le 10\%$	$\le 10\%$
C	$\le 20\%$	$\le 20\%$

Table 3. Table 3: Accuracy for polynomial degree P = 3 , 5 , 7 𝑃 3 5 7 P=3,5,7 .

polynomial degree	Envelop Misfit	Phase Misfit	Accuracy level
$P=3$	$\le 20.6\%$	$\le 8\%$	C
$P=5$	$\le 1.8\%$	$\le 1\%$	A
$P=7$	$\le 1.5\%$	$\le 0.6\%$	A

Table 4. Table 4: The density and elastic constants for the crystal Apatite. With c 55 = c 44 subscript 𝑐 55 subscript 𝑐 44 c_{55}=c_{44} c 66 = ( c 11 − c 12 ) / 2 subscript 𝑐 66 subscript 𝑐 11 subscript 𝑐 12 2 c_{66}=\left(c_{11}-c_{12}\right)/2

$\rho$	$c_{11}$	$c_{12}$	$c_{13}$	$c_{23}$	$c_{22}$	$c_{33}$	$c_{44}$
$3190$ kg/m3	$167$ GPa	$13.1$ GPa	$66$ GPa	$66$ GPa	$167$ GPa	$140$ GPa	$66.3$ GPa

Equations

$$\begin{array}{r}{\mathbf{P}}^{-1}\frac{\partial \mathbf{Q}}{\partial t}=\sum _{\xi =x,y,z}{\mathbf{A}}_{\xi }\frac{\partial \mathbf{Q}}{\partial \xi },\end{array}$$

$$\mathbf{Q}(x,y,z,0)={\mathbf{Q}}^0(x,y,z)\in L^2(\Omega ),$$

$$\frac{dE(x,y,z,t)}{dxdydz}=\frac{1}{2}\left[{\mathbf{Q}}^T{\mathbf{P}}^{-1}\mathbf{Q}\right]>0,E(t)={\int }_{\Omega }dE(x,y,z,t)>0.$$

$${\int }_{\Omega }{\mathbf{\varphi }}^T{\mathbf{P}}^{-1}\frac{\partial \mathbf{Q}}{\partial t}dxdydz$$

$${\int }_{\Omega }{\mathbf{Q}}^T{\mathbf{P}}^{-1}\frac{\partial \mathbf{Q}}{\partial t}dxdydz=\frac{1}{2}{\oint }_{\Gamma }{\mathbf{Q}}^T\left(\sum _{\xi =x,y,z}{n}_{\xi }{\mathbf{A}}_{\xi }\right)\mathbf{Q}dS.$$

$$\frac{d}{dt}E(t)=0.$$

$$\mathbf{Q}\left(x,y,z,t\right)=\left[\begin{array}{c}\mathbf{v}(x,y,z,t)\\ \mathbf{\sigma }(x,y,z,t)\end{array}\right],$$

$${\mathbf{A}}_{\xi }=\left(\begin{array}{cc}{0}_3& {\mathbf{a}}_{\xi }\\ {\mathbf{a}}_{\xi }^T& 0_6\end{array}\right),{\mathbf{a}}_x=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\end{array}\right),{\mathbf{a}}_y=\left(\begin{array}{cccccc}0& 0& 0& 1& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),{\mathbf{a}}_z=\left(\begin{array}{cccccc}0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\\ 0& 0& 1& 0& 0& 0\end{array}\right),$$

$$\mathbf{P}=\left(\begin{array}{cc}{\rho }^{-1}1& 0\\ 0^T& \mathbf{C}\end{array}\right),1=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right),0=\left(\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right),$$

$$\frac{dE}{dxdydz}:=\frac{1}{2}[{\mathbf{Q}}^T{\mathbf{P}}^{-1}\mathbf{Q}]=\frac{\rho }{2}\left(v_x^2+v_y^2+v_z^2\right)+\frac{1}{2}{\mathbf{\sigma }}^T\mathbf{S}\mathbf{\sigma }>0,$$

$$\mathbf{C}=\left(\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& 0& 0& 0\\ c_{12}& c_{22}& c_{23}& 0& 0& 0\\ c_{13}& c_{23}& c_{33}& 0& 0& 0\\ 0& 0& 0& c_{44}& 0& 0\\ 0& 0& 0& 0& c_{55}& 0\\ 0& 0& 0& 0& 0& c_{66}\end{array}\right).$$

$$\mathbf{C}=\left(\begin{array}{cccccc}2\mu +\lambda & \lambda & \lambda & 0& 0& 0\\ \lambda & 2\mu +\lambda & \lambda & 0& 0& 0\\ \lambda & \lambda & 2\mu +\lambda & 0& 0& 0\\ 0& 0& 0& \mu & 0& 0\\ 0& 0& 0& 0& \mu & 0\\ 0& 0& 0& 0& 0& \mu \end{array}\right).$$

$${\mathbf{e}}_x={\left(1,0,0\right)}^T,{\mathbf{e}}_y={\left(0,1,0\right)}^T,{\mathbf{e}}_z={\left(0,0,1\right)}^T,$$

$$\mathbf{v}=\left(\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right),{\mathbf{T}}^{(\xi )}=\left(\begin{array}{c}{T}_x^{(\xi )}\\ T_y^{(\xi )}\\ T_z^{(\xi )}\end{array}\right)=\overline{\stackrel{ˉ}{\sigma }}{\mathbf{e}}_{\xi }={\mathbf{a}}_{\xi }\overline{\stackrel{ˉ}{\sigma }},\xi =x,y,z,\overline{\stackrel{ˉ}{\sigma }}=\left(\begin{array}{ccc}{\sigma }_{xx}& {\sigma }_{xy}& {\sigma }_{xz}\\ {\sigma }_{xy}& {\sigma }_{yy}& {\sigma }_{yz}\\ {\sigma }_{xz}& {\sigma }_{yz}& {\sigma }_{zz}\end{array}\right).$$

$${\tilde{\mathbf{A}}}_{\xi }=\mathbf{P}{\mathbf{A}}_{\xi }.$$

$$Z_{p\xi }{v}_{\xi }\mp T_{\xi }^{(\xi )},Z_{sh\xi }{v}_{\eta }\mp T_{\eta }^{(\xi )},Z_{sv\xi }{v}_{\theta }\mp T_{\theta }^{(\xi )},$$

$$c_{px}=\sqrt{\frac{c_{11}}{\rho }},c_{shx}=\sqrt{\frac{c_{44}}{\rho }},c_{svx}=\sqrt{\frac{c_{55}}{\rho }},$$

$$c_{py}=\sqrt{\frac{c_{22}}{\rho }},c_{shy}=\sqrt{\frac{c_{66}}{\rho }},c_{svy}=\sqrt{\frac{c_{44}}{\rho }},$$

$$c_{pz}=\sqrt{\frac{c_{33}}{\rho }},c_{shz}=\sqrt{\frac{c_{55}}{\rho }},c_{svz}=\sqrt{\frac{c_{66}}{\rho }}.$$

$$c_{p\xi }=c_p=\sqrt{\frac{\lambda +2\mu }{\rho }},c_{sh\xi }=c_{sv\xi }=c_s=\sqrt{\frac{\mu }{\rho }}.$$

$${\mathbf{P}}^{-1}\frac{\partial \mathbf{Q}}{\partial t}=\sum _{\xi =x,y,z}{\mathbf{A}}_{\xi }\frac{\partial \mathbf{Q}}{\partial \xi }=\text{divergence}\mathbf{F}\left(\mathbf{Q}\right)+\sum _{\xi =x,y,z}{\mathbf{B}}_{\xi }\left(\text{gradient}\mathbf{Q}\right),\mathbf{P}={\mathbf{P}}^T>0,$$

$$\sum _{\xi =x,y,z}{\mathbf{A}}_{\xi }\frac{\partial \mathbf{Q}}{\partial \xi }=\text{divergence}\mathbf{F}\left(\mathbf{Q}\right)+\sum _{\xi =x,y,z}{\mathbf{B}}_{\xi }\left(\text{gradient}\mathbf{Q}\right),$$

$${\mathbf{Q}}^T{\mathbf{B}}_{\xi }\left(\text{gradient}\mathbf{Q}\right)-\frac{\partial {\mathbf{Q}}^T}{\partial \xi }{\mathbf{F}}_{\xi }\left(\mathbf{Q}\right)=0.$$

$${\mathbf{F}}_{\xi }\left(\mathbf{Q}\right)=\left(\begin{array}{c}{e}_{\xi x}{\sigma }_{xx}+e_{\xi y}{\sigma }_{xy}+e_{\xi z}{\sigma }_{xz}\\ e_{\xi x}{\sigma }_{xy}+e_{\xi y}{\sigma }_{yy}+e_{\xi z}{\sigma }_{yz}\\ e_{\xi x}{\sigma }_{xz}+e_{\xi y}{\sigma }_{yz}+e_{\xi z}{\sigma }_{zz}\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right),{\mathbf{B}}_{\xi }\left(\text{gradient}\mathbf{Q}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ e_{\xi x}\frac{\partial v_x}{\partial \xi }\\ e_{\xi y}\frac{\partial v_y}{\partial \xi }\\ e_{\xi z}\frac{\partial v_z}{\partial \xi }\\ e_{\xi y}\frac{\partial v_x}{\partial \xi }+e_{\xi x}\frac{\partial v_y}{\partial \xi }\\ e_{\xi z}\frac{\partial v_x}{\partial \xi }+e_{\xi x}\frac{\partial v_z}{\partial \xi }\\ e_{\xi z}\frac{\partial v_y}{\partial \xi }+e_{\xi y}\frac{\partial v_z}{\partial \xi }\end{array}\right),{\mathbf{e}}_{\xi }={\left(e_{\xi x},e_{\xi y},e_{\xi z}\right)}^T,$$

$${\mathbf{Q}}^T{\mathbf{B}}_{\xi }\left(\text{gradient}\mathbf{Q}\right)-\frac{\partial {\mathbf{Q}}^T}{\partial \xi }{\mathbf{F}}_{\xi }\left(\mathbf{Q}\right)=0.$$

$${\mathbf{Q}}^T{\mathbf{F}}_{\xi }\left(\mathbf{Q}\right)={\mathbf{v}}^T\left(\overline{\stackrel{ˉ}{\sigma }}{\mathbf{e}}_{\xi }\right)={\mathbf{v}}^T{\mathbf{T}}^{(\xi )}.$$

$$\frac{d}{dt}E\left(t\right)={\oint }_{\Gamma }{\mathbf{v}}^T\mathbf{T}dS,E\left(t\right)={\int }_{\Omega }dE>0.$$

$${\int }_{\Omega }{\mathbf{Q}}^T{\mathbf{P}}^{-1}\frac{\partial }{\partial t}\mathbf{Q}dxdydz={\int }_{\Omega }{\mathbf{Q}}^T\text{divergence}\mathbf{F}\left(\mathbf{Q}\right)dxdydz+\sum _{\xi =x,y,z}{\int }_{\Omega }{\mathbf{Q}}^T{\mathbf{B}}_{\xi }\left(\text{gradient}\mathbf{Q}\right)dxdydz.$$

$$\frac{d}{dt}E\left(t\right)=\sum _{\xi =x,y,z}{\int }_{\Omega }\left[{\mathbf{Q}}^T{\mathbf{B}}_{\xi }\left(\text{gradient}\mathbf{Q}\right)-\frac{\partial {\mathbf{Q}}^T}{\partial \xi }{\mathbf{F}}_{\xi }\left(\mathbf{Q}\right)\right]dxdydz+{\oint }_{\Gamma }{\mathbf{v}}^T\mathbf{T}dS.$$

$$\frac{d}{dt}E\left(t\right)={\oint }_{\Gamma }{\mathbf{v}}^T\mathbf{T}dS.$$

Full text

1 Introduction

Numerical algorithms based on the discontinuous Galerkin (DG) method [1, 2, 3, 4, 5] have shown to be flexible, high order accurate, provably stable, and well suited for complex large scale wave propagation problems [5, 6]. In computational seismology, they have been successfully applied to extreme-scale simulations of grand-challenge scenarios, often exploiting the largest-available supercomputers [7, 8, 9, 10].

In this study, we will pay special attention to scattering of linear elastic waves in heterogeneous, isotropic and anisotropic, solid Earth models with complex free surface topography. Accurate and efficient numerical simulation of seismic surface and interface waves, and scattering of high frequency waves by complex nonplanar topography are critical for assessing and quantifying seismic risks and hazards [11, 12, 13]. Surface and interface waves [14] are often the largest amplitude waves modes. On a regional scale, physics-based 3D ground motion simulations which include realistic three-dimensional Earth structure and topography are now able to resolve frequencies that are relevant for building response (from static displacements at zero frequency up to 5-10 Hz or higher, [15, 16, e.g.]), however, many computational studies are forced to use highly smoothed representations of surface topography due to limitations either in terms of algorithmic restrictions (e.g., spatial discretisations limited to structured and or hexahedral approaches, limited flexibility of free-surface boundary conditions, limited applicability of graded meshes) and/or computational efficiency. For efficient representation of complex geometries, in this study, we use adaptive boundary conforming curvilinear elements. The elastic wave equation is transformed from Cartesian coordinates to curvilinear coordinates using a structure preserving coordinate transformation. The transformation is local within the element, and inside the element the curvilinear elements have logical Cartesian coordinates. Essentially, we solve the equations on adaptive Cartesian meshes, and the complex geometries are moved into variable metric terms and are used to define transformed variable material parameters. Our numerical scheme allows adaptive non-conforming refinement of the hexahedral elements, which enables its implementation on tree- or block-structured Cartesian meshes. We implemented the scheme in the ExaHyPE engine for solving hyperbolic PDE systems [17], which realises DG with ADER time stepping [18]. The Peano adaptive mesh refinement framework [19, 20] provides dynamically adaptive tree-structured Cartesian meshes and parallelisation in shared and distributed memory systems.

DG and physics-based numerical fluxes

The DG method is an increasingly attractive method for approximating partial differential equations (PDEs), and has been successfully used in computational seismology applications including geometrically and rheologically complicated wave propagation and dynamic rupture simulations [21, 22, 23, 24, 25, e.g.,]. An important component of the DG method is the numerical flux [26, 27]. For hyperbolic PDEs this flux is based on approximate or exact solutions of the Riemann problem [28, 29]. The choice of a numerical flux and (approximate) solutions of the Riemann problem are critical for accuracy and stability of the DG method [30, 31, 32]. The Rusanov flux [29] (also called local Lax-Friedrichs flux) is widely used, because of its simplicity and robustness. However, the Rusanov flux might be inappropriate for simulating seismic surface waves [33]. The Godunov flux with the exact solution to the Riemann problem has been demonstrated as a reasonable choice for seismological applications [34, 35]. However, the Godunov flux requires a complete eigenvector and eigenvalue decomposition of the coefficient matrices of the spatial operator. In general anisotropic media, the eigen-decomposition can be nontrivial. The eigen-decomposition will become even more cumbersome if elements are curved. In [36, 37] numerical penalty fluxes are introduced, penalizing the normal components of the primitive variables across inter-element faces, and avoiding the eigen-decomposition of coefficient matrices. However, many dynamic boundary conditions such as empirical constitutive friction laws [38, 39, 40] are formulated using derived quantities, such as tractions in local coordinates. We argue that it will be more natural to anchor numerical fluxes in elastic solids on derived quantities such as local velocity and traction vectors, rotated into local orthogonal coordinates. This will enable the development of a unified provably stable and robust adaptive DG framework in complex geometries for the numerical treatment of 1) nonlinear frictional sliding in elastic solids, 2) for coupling classical DG inter-element interfaces in elastic solids where slip is not permitted, and 3) numerical enforcement of external well-posed boundary conditions modeling various solid mechanics and geophysical phenomena.

In [33], we introduce (in 1D and 2D) an alternative and accurate approach to couple locally adjacent DG elements, using physical conditions such as friction. The nascent physics-based numerical flux obeys the eigen-structure of the PDE and the underlying physics at the internal and external DG element boundaries in a provably stable manner. Our formulation does not require a complete eigenvector-eigenvalue decomposition of the spatial coefficient matrices and can be easily adapted to model linear and nonlinear boundary and interface wave phenomena. In spirit, our approach is analogous to the method used in a finite difference framework [41] to model frictional sliding during dynamic earthquake rupture. However, static and/or dynamic adaptive mesh refinement in a finite difference setting is arduous.

Extension to 3D curvilinear elements

In this study, we will extend the physics-based numerical flux [33] to 3D geometrically complex elastic media. We extract tractions and particle velocities on an element face in Cartesian coordinates and rotate them into the curvilinear coordinates. The rotated tractions and particle velocities, in conjunction with the material impedance, can be used to extract characteristics, such as plane shear waves and compressional waves, propagating along the boundary surface. Then, we construct boundary and interface data by solving a Riemann-like problem and constrain the solutions against the physical conditions acting at the element faces. The physics based fluctuations are constructed by penalisation of boundary and interface data against the incoming characteristics. These fluctuations are then appended to the discrete equations with physically motivated penalty weights chosen such that the semi-discrete approximation satisfies an energy estimate analogous to the continuous energy estimate. The energy estimate proves the asymptotic stability of the semi-discrete approximation.

The semi-discrete DG approximation is integrated in time using the Arbitrary DERivative (ADER) time integration [42]. The ADER time discretisation is summarised in A. The cell-wise local character of DG can be readily combined with an ADER scheme leading to high-order accuracy in time within a single-step. The combination of the DG approximation in space and the ADER time integration is often referred to as the ADER-DG scheme [34]. For the ADER-DG scheme, the numerical flux fluctuation is evaluated only once for any order of accuracy. The implication is that most of the computations are performed within the element to evaluate spatial derivatives.

We will here present 3D ADER-DG numerical simulations verifying accuracy and stability of the method, using community developed benchmark problems [43, 44, 45, 46] and complex geologically constrained geometries. All of our simulations employ provable stable PML boundary conditions [47, 48] to effectively prevent artificial boundary conditions from contaminating the simulations. A stable numerical implementation of the PML for 3D linear elastodynamics is nontrivial and allows the generation of high quality seismograms. We will perform error analysis and compute error parameters relevant to computational seismology [43, 44].

The remaining part of the paper will proceed as follows. In the next section, we introduce a general model for linear elastodynamics. Curvilinear coordinates and structure preserving curvilinear transformations are presented in section 3. In section 4, we present physical boundary and interface conditions and derive energy estimates. In section 5, the physics based numerical flux and the algorithms for solving the Riemann problem in heterogeneous media and arbitrary curvilinear coordinates are presented. Numerical experiments based on the implementation in ExaHyPE are presented in section 8 verifying accuracy and stability, and demonstrating the potentials of the method in simulating complex wave phenomena. In the last section we draw conclusions and suggest future work.

2 First order linear hyperbolic PDE

Consider the 3D first order linear hyperbolic system in a source–free heterogeneous medium

[TABLE]

where $\mathbf{P}=\mathbf{P}^{T}$ with $\mathbf{Q}^{T}\mathbf{P}\mathbf{Q}>0$ and $\mathbf{A}_{\xi}=\mathbf{A}_{\xi}^{T}$, $\xi=x,y,z$. Here, $t\geq 0$ denotes time and $(x,y,z)\in\Omega$ are the Cartesian coordinates of the spatial domain $\Omega\subset\mathbb{R}^{3}$. At the initial time, $t=0$, we set the initial condition

[TABLE]

belonging to the space of square integrable functions. In general, the symmetric positive definite matrix $\mathbf{P}$ depends on the spatial coordinates $x,y,z$, and encodes the material parameters of the underlying medium. The non-dimensional constant matrices $\mathbf{A}_{\xi}$ encapsulate the underlying linear conservation law and the corresponding linear constitutive relation.

We introduce the energy density and the physical energy, defined by

[TABLE]

Let $\Gamma$ denote the boundary of the domain $\Omega$, and $\mathbf{n}=\left(n_{x},n_{y},n_{z}\right)^{T}$ the outward unit normal on the boundary. We multiply (1) with $\boldsymbol{\phi}^{T}(x,y,z)$ from the left, where $\boldsymbol{\phi}(x,y,z)\in L^{2}\left(\Omega\right)$ is an arbitrary test function, and integrate over the whole spatial domain, $\Omega$. We have

[TABLE]

In the right hand side of (4), integrating-by-parts, and using the fact that the coefficient matrices are constant and symmetric, $\mathbf{A}_{\xi}=\mathbf{A}_{\xi}^{T}$, gives

[TABLE]

Replacing $\boldsymbol{\phi}$ with $\mathbf{Q}$ in (5), in the right hand side, the volume terms vanish, having

[TABLE]

The decay condition, $|\mathbf{Q}|\to 0$ at the boundary $(x,y,z)\in\Gamma$, yields the energy equation

[TABLE]

The energy is conserved, $E(t)=E(0)$ for all $t\geq 0$.

Depending on the coefficient matrices $\mathbf{P},\mathbf{A}_{\xi}$ the system (1) can describe acoustic waves, electromagnetic waves, linear MHD waves, linear elastic waves, or the interaction of acousto–elastic waves propagating in a heterogeneous medium. Our general interest is the development of a provably stable DG approximations for first order linear hyperbolic PDEs of the form (1), in heterogeneous and geometrically complex elastic solids, subject to well-posed boundary conditions.

In this study, we will pay special attention to scattering of linear elastic waves in heterogeneous, isotropic and anisotropic, elastic solids with complex free surface topography. We note, that with limited modifications the method we develop can be extended to other problems modeled by (1) in future work.

2.1 Waves in elastic solids

To describe wave propagation in elastic solids, we introduce the unknown wave fields

[TABLE]

with the particle velocity vector, $\mathbf{v}(x,y,z,t)=\left[v_{x},v_{y},v_{z}\right]^{T}$, and the stress vector,

$\boldsymbol{\sigma}(x,y,z,t)=\left[\sigma_{xx},\sigma_{yy},\sigma_{zz},\sigma_{xy},\sigma_{xz},\sigma_{yz}\right]^{T}$. The symmetric constant coefficient matrices $\mathbf{A}_{\xi}$ describing the conservation of momentum and the constitutive relation, defined by Hooke’s law, are given by

[TABLE]

where $\mathbf{0}_{3}$ and $\mathbf{0}_{6}$ are the $3$-by-$3$ and $6$-by-$6$ zero matrices.

The symmetric positive definite material parameter matrix $\mathbf{P}$ is defined by

[TABLE]

where $\rho(x,y,z)>0$ is the density of the medium, and $\mathbf{C}=\mathbf{C}^{T}>0$ is the symmetric positive definite matrix of elastic constants. With the unknown wave fields prescribed by (8), and the coefficient matrices defined in (9) and (10), the first three equations in (1) are the conservation of momentum and the last six equations are the time derivative of the constitutive relation, defined by Hooke’s law, relating stress fields to strains where the constant of proportionality is the stiffness matrix of elastic coefficients $\mathbf{C}$.

The mechanical energy density is the sum of the kinetic energy density and the strain energy density

[TABLE]

where $\mathbf{S}=\mathbf{C}^{-1}$ is the compliance matrix.

In a general anisotropic medium the stiffness matrix $\mathbf{C}$ is described by 21 independent elastic coefficients. In an orthotropic anisotropic medium the stiffness matrix has 9 independent elements,

[TABLE]

In the isotropic case, the medium is described by two independent elastic coefficients, the Lamé parameters $\mu>0$, $\lambda>-\mu$, thus we have

[TABLE]

However, as shown in [49] (pages 241–243), for strong ellipticity we must have $\mu>0$, $\lambda>-1/2\mu$.

We introduce the 3D canonical basis vectors

[TABLE]

and define the Cartesian components of the velocity and traction vectors

[TABLE]

Next, we define the coefficient matrices

[TABLE]

In the absence of boundaries and discontinuous interfaces, the elastic wave equation supports two families of solutions, primary (p-waves) and secondary (s-waves) waves. The p-wave and s-wave modes are related to the nontrivial eigenfunctions of $\widetilde{\mathbf{A}}_{\xi}$, and the nontrivial eigenvalues of $\widetilde{\mathbf{A}}_{\xi}$ are $\pm c_{p\xi},\quad\pm c_{sh\xi},\quad\pm c_{sv\xi},$ and correspond to the p-wave and s-wave speeds. Note that $c_{p\xi}>0$ are the p-wave speeds, $c_{sh\xi}>0$ are the wave speeds of the horizontally polarized s-wave and $c_{sv\xi}>0$ are the wave speeds of the vertically polarized s-wave. The negative and positive going p-wave and s-wave modes are given by

[TABLE]

where $Z_{p\xi}=\rho c_{p\xi}$, $Z_{sh\xi}=\rho c_{sh\xi}$, $Z_{sv\xi}=\rho c_{sv\xi}$ are the impedances. Here, $\xi,\eta,\theta=x,y,z$, $\eta\neq\xi$, and $\theta\neq\xi,\eta$. The wave modes defined in (17) are the plane p-waves and plane s-waves propagating along the ${\xi}$-axis, and are related to the 1D Riemann invariant. The eigenvalues of the matrices $\widetilde{\mathbf{A}}_{\xi}$ can be easily determined. For example in orthotropic anisotropic media, with $\mathbf{C}$ defined in (12), the eigenvalues are given by

[TABLE]

Note that for the same wave mode the wave speed can vary in all directions. In isotropic media, the eigenvalues are uniform in all directions

[TABLE]

Thus, in isotropic media a wave mode propagates with a uniform wave speed in all directions. In particular, the vertically polarized and the horizontally polarized s-waves have identical wave-speed $c_{sh\xi}=c_{sv\xi}=c_{s}$ in all directions.

We recall our main interest being the development of an energy stable DG approximation of the equation of motion (1) defined by (8), (9) and (10), in heterogeneous and geometrically complex elastic solid with complicated nonplanar free-surface topography. To this end, we will next introduce the anti-symmetric split form. This is necessary for the development of provably stable approximations on curvilinear meshes for complex geometries and generally arbitrary heterogeneous media.

2.2 Anti-symmetric splitting

To enable effective numerical treatments, we introduce the split form of the equation of motion (1) defined by

[TABLE]

for some flux function $\mathbf{F}\left(\mathbf{Q}\right)=[\mathbf{F}_{x}\left(\mathbf{Q}\right),\mathbf{F}_{y}\left(\mathbf{Q}\right),\mathbf{F}_{z}\left(\mathbf{Q}\right)]^{T}$ and coefficient matrices $\mathbf{B}_{\xi}$.

Definition 1**.**

Consider the equation of motion (1). Suppose that the spatial operator can be split into:

[TABLE]

*where the first term, with $\divergence\mathbf{F}\left(\mathbf{Q}\right)$, is called the conservative flux term and the second term, with $\mathbf{B}_{\xi}\left(\gradient\mathbf{Q}\right)$, is the non-conservative-products flux term. The split operators, given in (21), are anti-symmetric if *

[TABLE]

The anti-symmetric split-form, (20) with (22), can be useful when designing provably stable approximations on curvilinear meshes for complex geometries, and variable material properties.

There are several ways of casting a hyperbolic PDE, such as (1), in the anti-symmetric form, see the examples in [50, 51, 52, 53, 54]. Here, our choice of the anti-symmetric form is motivated by the underlying physics. For linear elasticity we will use

[TABLE]

where the conservative flux comes from conservation of momentum, that is Newton’s second law of motion. The non-conservative product term comes from Hooke’s law, the constitutive relation relating stress and strain. In a Cartesian coordinate system, $\mathbf{e}_{\xi}$ are the canonical bases defined in (14). In general curvilinear coordinates the bases $\mathbf{e}_{\xi}$ are arbitrary nonzero vectors, with $|\mathbf{e}_{\xi}|>0$.

Lemma 1**.**

Consider the split form of the equation (20), with the split operators defined in (23). For arbitrary bases vectors $\mathbf{e}_{\xi}$, with $|\mathbf{e}_{\xi}|>0$, the conservative flux term and the non-conservative-products flux term satisfy the anti-symmetric property (22), that is

[TABLE]

Proof.

The proof of the lemma follows from direct calculations, by using (23) and evaluating the products in (22). ∎

Another important consequence of the choice of the split operators (23) is the fact

[TABLE]

Theorem 1**.**

*Consider the split form of equation of motion (20) for the linear elastic wave equation, defined by (23) in the spatial domain $(x,y,z)\in\Omega$ where $\Gamma$ denotes the boundary of the domain and $\mathbf{n}$ is the outward unit normal on the boundary. Let $\mathbf{T}=\bar{\bar{\sigma}}\mathbf{n}$ denote the traction on the boundary, and the energy density $dE$ defined by (11) for an elastic medium. The solutions of the anti-symmetric split form (20) satisfy *

[TABLE]

Proof.

From the left, multiply the split form of equation of motion (20) by the transpose of the solution $\mathbf{Q}^{T}$ and integrate over the whole domain $\Omega$, having

[TABLE]

On the left hand side of (42) we recognize time derivative of the energy $E\left(t\right)$. On the right hand side of (26), integrate–by–parts the first (conservative flux) term, we have

[TABLE]

By the anti-symmetric property (22), the volume terms vanish having

[TABLE]

∎

In the coming section below, we will introduce curvilinear coordinates and transformations to model geometrically complex elastic solids.

3 Curvilinear coordinates and structure preserving coordinate transformations

To simplify the presentation, we consider two DG elements filled with heterogeneous elastic media separated by an interface at $\left(\widetilde{x}\left(y,z\right),y,z\right)$, where $\widetilde{x}\left(y,z\right)$ is an arbitrarily smooth level surface describing the surface of the interface. The interface can be a pre-existing fault in the medium, or the interface between two adjacent DG elements.

For geometrically complex models, numerical treatments can be simplified if there is a smooth coordinate transformation from the physical space to a reference element, the unit cube $\left(q,r,s\right)\in\widetilde{\Omega}=[0,1]\times[0,1]\times[0,1]$. We split the elements along the interface and map each element to the reference element, defined by

[TABLE]

Note that the mapping and coordinate transformations are element local.

Let us denote the metric coefficients and the Jacobian of the transformation (29) by

[TABLE]

[TABLE]

Here, the subscripts denote partial derivatives, that is $x_{q}=\frac{\partial x}{\partial q}$, $x_{r}=\frac{\partial x}{\partial r}$, etc.

To construct structure preserving coordinate transformation, we will use two different transformations of the spatial derivatives in the transformed coordinates [41], the conservative form:

[TABLE]

and the non-conservative form:

[TABLE]

Note that in the continuous setting, the transformed derivatives (30) and (31) are mathematically equivalent. In nontrivial geometries, discrete approximations of the transformed derivatives (30) and (31) will yield two different discrete spatial operators.

In order to preserve the anti-symmetric structure, in (20), it becomes natural to transform the derivatives in the conservative flux term using the conservative transformation (30), and the derivatives in the non-conservative-products term using the non-conservative transformation (31). The elastic wave equation (20) in the transformed curvilinear coordinates $(q,r,s)$ is

[TABLE]

Here, $\gradient=\left(\partial/\partial q,\partial/\partial r,\partial/\partial s\right)^{T}$ is the gradient operator, $\widetilde{\mathbf{P}}=J^{-1}\mathbf{P}$, where $\mathbf{P}$ is the material matrix defined in (9), and

[TABLE]

where $\xi=q,r,s$, with $\xi_{\eta}=\partial\xi/\partial\eta$ and $\eta=x,y,z$.

Note that in (32) all spatial derivatives for the stress fields are transformed using the conservative form (30) and all spatial derivatives for the velocity fields are transformed using the non-conservative form (31). When discrete approximations are introduced, this is crucial in order minimize the number of floating point operations, and also prove numerical stability [41].

We will show that the transformed equation of motion, (32) with (33), preserves the anti-symmetric property (22).

Lemma 2**.**

Consider the transformed equation of motion (32), in curvilinear coordinates, with the conservative flux terms and non-conservative-products flux terms given by (33). The corresponding spatial operators satisfy the anti-symmetric property (22), that is

[TABLE]

Proof.

With $\mathbf{e}_{\xi}=J\left({\xi}_{x},{\xi}_{y},{\xi}_{z}\right)^{T}$, the proof of Lemma 2 follows directly from Lemma 1. ∎

Let $\xi=q,r,s$, and $\Gamma$ denote a boundary face at $\xi=0$ or $\xi=1$. The the positively pointing unit normals on the boundary are given by

[TABLE]

Note again that

[TABLE]

where $\mathbf{v}=[v_{x},v_{y},v_{z}]^{T}$ is the velocity vector, and $\mathbf{T}=[T_{x},T_{y},T_{z}]^{T}=\bar{\bar{\sigma}}\mathbf{n}$ is the traction vector.

Let us introduce the reference boundary surface $\widetilde{\Gamma}=[0,1]\times[0,1]$, and define the boundary term

[TABLE]

Another important consequence of the the transformed anti-symmetric split form (32) with (33) is the following lemma

Lemma 3**.**

*Consider the conservative flux term defined in (33). We have *

[TABLE]

where the boundary term $\mathrm{BTs}\left(v,T\right)$ is defined in (3)

Proof.

Consider

[TABLE]

and integrate-by-parts, we have

[TABLE]

Using the fact (35), and the boundary term $\mathrm{BTs}\left(v,T\right)$ defined in (3) completes the proof. ∎

Introduce the energy density in the transformed space

[TABLE]

Analogous to Theorem 1, we also have

Theorem 2**.**

*Consider the transformed equation of motion (32), in curvilinear coordinates, $(q,r,s)\in\widetilde{\Omega}$, with the flux terms and non-conservative products terms given by (33). Let $\xi=q,r,s$, and $\Gamma$ denote a boundary face at $\xi=0$ or $\xi=1$, where $\mathbf{T}=\bar{\bar{\sigma}}\mathbf{n}$ is the traction vector on the boundary. The solutions of the transformed equation (32) satisfy *

[TABLE]

*where $\mathrm{BTs}\left(v,T\right)$ is the boundary term defined in (3). *

Proof.

From the left, multiply the transformed equation of motion (32) by $\mathbf{Q}^{T}$ and integrate over the whole domain $\widetilde{\Omega}$, having

[TABLE]

On the left hand side of (42) we recognize time derivative of the energy. On the right hand side of (42), we use Lemma 3, and replace the conservative flux term with (37), we have

[TABLE]

where $\mathrm{BTs}$ is the boundary term defined in (3). By Lemma 2, the volume terms in (43) vanish, having

[TABLE]

∎

The boundary term, $\mathrm{BTs}$, is the rate of the work done by the traction, $\mathbf{T}=\bar{\bar{\sigma}}\mathbf{n}$, against the boundary. Note that the energy rate is controlled by the boundary term, $\mathrm{BTs}$. In a bounded domain, well-posed boundary conditions are designed such that the boundary term is never positive, $\mathrm{BTs}\leq 0$.

4 Boundary and interface conditions

We will now introduce physical boundary and interface conditions prescribed at element faces. In Figure 1, for each of the two blocks, there are five external boundaries, and one internal boundary connecting the two blocks. First, we will consider boundary conditions posed at external element boundaries, and proceed later to interface conditions acting at the internal element face, connecting the two locally adjacent elements.

4.1 Boundary conditions

Here, we present linear well-posed boundary conditions closing the external boundaries. Consider the 5 external boundaries of each elastic block, see Figure 1. To define well-posed boundary conditions, we introduce the basis vectors denoted by $\mathbf{n}=\left(n_{x},n_{y},n_{z}\right),\mathbf{m}=\left(m_{x},m_{y},m_{z}\right),\textbf{{l}}=\left(l_{x},l_{y},l_{z}\right)$, where, $\mathbf{n}$ is the unit normal, defined by (34), on the boundary pointing in the positive $\xi$-direction, with $\xi=q,r,s$.

Note that the unit vectors are locally defined on each point on the boundary, see also Figure 2. Given the unit normal $\mathbf{n}$, defined by (34), we can construct the other two basis vectors $\mathbf{l},\mathbf{m}$, as follows:

[TABLE]

Here, $\mathbf{m}_{0}\left(\neq\pm\mathbf{n}\right)$ is an arbitrary unit vector. Let the local impedances at the boundary be denoted by $Z_{\eta},\eta=l,m,n,$ where $Z_{n}=\rho c_{n}$ is the p–wave impedance and $Z_{m}=\rho c_{m}$, $Z_{l}=\rho c_{l}$ are the s–wave impedances. Here, $c_{n},c_{m},c_{l}$ are the corresponding effective wavespeeds defined by

[TABLE]

In anisotropic media with a geometrically complex curvilinear coordinate system, the effective wavespeeds depend on the orientation of the normal vector $\mathbf{n}$. Note that ${n_{x}^{2}+n_{y}^{2}+n_{z}^{2}}=1$. Thus, in an isotropic medium the effective wavespeeds are given by $c_{n}=c_{p},\quad c_{m}=c_{l}=c_{s}.$

On each point on the boundary, we denote the particle velocity vector, traction vector, and the local rotation matrix on the boundary by

[TABLE]

where $\det{(\mathbf{R})}\neq 0$ and $\mathbf{R}^{-1}=\mathbf{R}^{T}$.

Next, rotate the particle velocity and traction vectors into the local orthonormal basis, $\mathbf{l}$ , $\mathbf{m}$ and $\mathbf{n}$, having

[TABLE]

The corresponding in and out of the domain characteristics at the boundary are

[TABLE]

Here, the characteristics defined in (49) are plane p–waves and plane s–waves propagating along the normal vector $\mathbf{n}$ on the boundary.

At the boundary $\xi=1$ $(\xi=0)$, if $Z_{\eta}>0$ then ${q}_{\eta}$ (${p}_{\eta}$) are the characteristics going into the domain and ${p_{\eta}}$ (${q_{\eta}}$) the characteristics going out of the domain. The number of boundary conditions must correspond to the number of characteristics going into the domain, see [41, 55]. We consider linear boundary conditions,

[TABLE]

where the reflection coefficients $\gamma_{\eta}$ are real numbers with $0\leq|\gamma_{\eta}|\leq 1$. The boundary conditions (50) specify the ingoing characteristics on the boundary in terms of the outgoing characteristics. In an elastic medium, we have $Z_{\eta}>0$ for all $\eta=l,m,n$, and there are three boundary conditions at each boundaries $\xi=1$, ($\xi=0$). The boundary condition (50), can describe several physical situations. We have a free-surface boundary condition if $\gamma_{\eta}=1$, an absorbing boundary condition if $\gamma_{\eta}=0$ and a clamped boundary condition if $\gamma_{\eta}=-1$. For later use in deriving energy estimates, we note that

[TABLE]

Lemma 4**.**

Consider the well-posed boundary conditions (50) with $|\gamma_{\eta}|\leq 1$. The boundary term $\mathrm{BTs}$ defined in (3) is negative semi-definite, $\mathrm{BTs}\leq 0$, for all $Z_{\eta}>0$.

Proof.

With $\mathbf{v}^{T}\mathbf{T}=\left(\mathbf{R}\mathbf{v}\right)^{T}\left(\mathbf{R}\mathbf{T}\right)=\sum_{\eta=l,m,n}{v_{\eta}T_{\eta}}$, the boundary term $\mathrm{BTs}\left(v,T\right)$ defined in (3) can be written as

[TABLE]

Finally, the identity (4.1) completes the proof of the lemma. ∎

Using the energy method we can now prove:

Theorem 3**.**

*Consider the transformed equation of motion (32) subject to the boundary condition (50), with $|\gamma|\leq 1$. The solutions of the transformed equation (32) subject to the boundary condition (50) satisfy *

[TABLE]

Proof.

The proof of Theorem 3 follows the same steps as in the proof of Theorem 2, arriving at (44). We complete the proof using Lemma 4, and by ensuring that the boundary terms are never positive, $\mathrm{BTs}\leq 0$. ∎

4.2 Interface conditions

We introduce physical interface conditions acting at internal DG elements boundaries, in elastic solids. These physical interface conditions will connect two adjacent elements elastic media. One objective of this study is to use the physical conditions to patch DG elements together [33]. Consider the interface, as in Figure 1 and denote the corresponding fields and material parameters in the positive/negative sides of the interface with the superscripts $+/-$. To define the interface conditions we rotate the particle velocity vector and the traction vector on the boundary into the local orthogonal coordinates $\mathbf{n},\mathbf{m},\mathbf{l}$, as in (48), having

[TABLE]

We define the jumping condition in the velocity fields as

[TABLE]

As before, to ensure well–posedness, the number of interface conditions must be equal to the number of ingoing characteristics at the interface.

We consider two elastic solids in a locked contact, with $Z_{\eta}^{\pm}>0$ for all $\eta=l,m,n$. The interface is locked, that is, there is no opening/gap, no inter-penetration and no slip. As shown in figure 3, there are a total of $6$ (3 characteristics going into the negative element and 3 characteristics going into the positive element) ingoing characteristics at the interface. Therefore, we will need exactly 6 conditions specifying the relationships of the fields across the interface. The interface conditions are force balance, and vanishing opening and slip velocities

[TABLE]

Now, we will show that the interface condition (55) with the equation of motion (32), conserve the total mechanical energy. Introduce the interface term

[TABLE]

Thus, by (55) we have ${T}_{\eta}[\![{{v}_{\eta}]\!]}=0$, thus, the interface term vanishes identically, $IT_{s}(v^{\pm},T^{\pm},Z^{\pm})\equiv 0$. In particular we can prove the theorem:

Theorem 4**.**

*Consider the modelling domain decomposed into two elements, $\Omega=\Omega_{-}\cup\Omega_{+}$, as in Figure 1, and denote the corresponding fields, material parameters and energies in the positive/negative sides of the interface with the superscripts $+/-$. Let each element be mapped to the reference element $\widetilde{\Omega}=[0,1]^{3}$, with the transformed equation of motion (32) defined on each element and subject to the interface conditions (55). The sum of the energies satisfies *

[TABLE]

Proof.

Using the energy method, from (42)–(44), we again have

[TABLE]

where $\mathrm{BTs}\left(v,T\right)$ is the boundary term defined in (3). Collecting contributions from both elements and ignoring all other boundaries, excepting the boundary at the interface, we have

[TABLE]

Using the fact (56) that the interface term vanishes identically, $IT_{s}(v^{\pm},T^{\pm},Z^{\pm})\equiv 0$, completes the proof. ∎

Theorems (3) and (4) prove that the corresponding IBVPs are well-posed and asymptotically stable. The challenge, for the DG scheme, is how to incorporate the boundary conditions and the interface conditions in a seamless and provably stable manner, for the discrete problem, in general heterogeneous media with complex geometries. To succeed, we will introduce hat-variables, so that we can simultaneously construct data for the velocity fields and traction fields, at internal and external element faces.

5 Hat-variables and physics based fluxes

We will now reformulate the boundary condition (50) and interface condition (55) by introducing transformed (hat-) variables so that we can simultaneously construct (numerical) boundary/interface data for particle velocities and tractions. The hat-variables encode the solution of the IBVP on the boundary/interface. The hat-variables will be constructed such that they preserve the amplitude of the outgoing characteristics and exactly satisfy the physical boundary conditions [41]. To be more specific, the hat-variables are solutions of the Riemann problem constrained against physical boundary conditions (50), and the interface condition (55). We refer the reader to [33] for more detailed discussion. Once the hat-variables are available, we construct physics based numerical flux fluctuations by penalizing data against the incoming characteristics (49) at the element faces.

5.1 Boundary data

For $Z_{\eta}>0$, we define the characteristics

[TABLE]

Here, $q_{\eta}$ are the left going characteristics and $p_{\eta}$ are the right going characteristics. We will construct boundary data which satisfy the physical boundary conditions (50) exactly and preserve the amplitude of the outgoing characteristics $q_{\eta}$ at $\xi=0$, and $p_{\eta}$ at $\xi=1$.

To begin, define the hat-variables preserving the amplitude of outgoing characteristics

[TABLE]

Since hat-variables also satisfy the physical boundary condition (50), we must have

[TABLE]

The algebraic problem for the hat-variables, defined by equations (61) and (62), has a unique solution, namely

[TABLE]

The expressions in (63) define a rule to update particle velocities and tractions on the physical boundaries $\xi=0,1$. That is

[TABLE]

The computation of the hat-variables for the external boundaries is summarized in Algorithm 1 below.

By construction, the hat-variables $\widehat{v}_{\eta},\widehat{T}_{\eta}$, satisfy the following algebraic identities:

[TABLE]

The algebraic identities (65a)–(65c) will be crucial in proving numerical stability. Please see [33] for more details.

Lemma 5**.**

Consider the boundary term $\mathrm{BTs}$ defined in (3), where $v_{\eta}=\widehat{v}_{\eta}$, $T_{\eta}=\widehat{T}_{\eta}$, with $\widehat{T}_{\eta},\widehat{v}_{\eta}$ defined in (63). The boundary term $\mathrm{BTs}$ is never positive for all $|\gamma_{\eta}|\leq 1$, that is $\mathrm{BTs}\left(\widehat{v}_{\eta},\widehat{T}_{\eta}\right)\leq 0$ where

[TABLE]

Proof.

With $Z_{\eta}>0$, for any $|\gamma_{\eta}|\leq 1$, from (65c) we have

[TABLE]

[TABLE]

We must have $\mathrm{BTs}\left(\widehat{v}_{\eta},\widehat{T}_{\eta}\right)\leq 0$ for all $|\gamma_{\eta}|\leq 1$ . ∎

Lemma 5 is completely analogous to Lemma 4. We will now formulate a result equivalent to Theorem 3.

Theorem 5**.**

*Consider the transformed equation of motion (32) subject to the boundary condition (64). The solutions of the transformed equation (32) with (64) satisfy *

[TABLE]

5.2 Interface data

To begin, define the outgoing characteristics at the interface

[TABLE]

where $Z_{\eta}^{\pm}>0$ are the impedances. We define the hat-variables preserving the amplitude of the outgoing characteristics at the interface

[TABLE]

The hat-variables also satisfy the interface conditions (55) exactly. For each setup, given ${q}_{\eta}^{+},{p}_{\eta}^{-}$, the procedure will solve (71), and (55) for the hat-variables, $\widehat{v}_{\eta}^{\pm},\widehat{T}_{\eta}^{\pm}$. We consider an interface separating two elastic solids. The hat-variables must satisfy (71) and the interface conditions, (55), that is, force balance, no opening and no slip conditions. Combining the two equations in (71) and ensuring force balance, $\widehat{T}_{\eta}^{-}=\widehat{T}_{\eta}^{+}=\widehat{T}_{\eta}$, we have

[TABLE]

Furthermore, enforcing no opening and no slip conditions, $[\![{\widehat{v}_{\eta}]\!]}=0$, in (72) gives

[TABLE]

We can now define explicitly the hat-variables corresponding to the particle velocities and tractions

[TABLE]

Note that we have equivalently redefined the physical interface condition (55) as follows

[TABLE]

The procedure to compute the hat-variables for the internal element boundaries is given in Algorithm 2.

By construction, the hat-variables $\widehat{v}_{\eta}^{\pm}$, $\widehat{T}_{\eta}^{\pm}$ satisfy the following algebraic identities:

[TABLE]

The identities (76a)–(76b) can be easily verified, see [33], and will be useful in the prove of numerical stability.

We can now formulate a result equivalent to Theorem 4.

Theorem 6**.**

*Consider the modeling domain decomposed into two elements, $\Omega=\Omega_{-}\cup\Omega_{+}$, as in Figure 1. Denote the corresponding fields, material parameters and energies in the positive/negative sides of the interface with the superscripts $+/-$. Let the transformed equation of motion (32) be defined on each element, with the two adjacent elements connected at the interface through the interface condition (75). The sum of the energies satisfies *

[TABLE]

5.3 Physics based flux fluctuations

The next step is to construct fluctuations by penalizing data, that is hat-variables, against the ingoing characteristics only. If $Z_{\eta}>0$, then we have

[TABLE]

The fluctuations are computed in the transformed coordinates $l,m,n$ . We will now rotate them to the physical coordinates $x,y,z$, having

[TABLE]

Note that

[TABLE]

The physics based flux fluctuations obeying the eigen–structure of the elastic wave equation are given by

[TABLE]

Note that

[TABLE]

The elemental weak form reads

[TABLE]

Here, the variable $\xi=q,r,s$, indicates the directions where the flux is computed. The flux fluctuation vectors, $\mathbf{FL},\mathbf{FR}$ enforce weakly the boundary conditions (64) and interface interface conditions (75) at the element faces. Note, however, that since we have not introduced any numerical approximations, the fluctuations vanish identically, ${G}_{\eta}^{\pm}=\widetilde{G}_{\eta}^{\pm}=0$, for the exact solutions that satisfy the PDE and the boundary and interface conditions, we have $\mathbf{FL}\equiv 0$, $\mathbf{FR}\equiv 0$.

Introduce the fluctuation term

[TABLE]

For the two elements model we introduce the external boundary terms

[TABLE]

and the interface term

[TABLE]

Remark 1**.**

Note that by construction the interface terms vanish identically $IT_{s}(\widehat{v}^{\pm},\widehat{T}^{\pm})\equiv 0$, by inspection the fluctuation terms are negative semi-definite $F_{luc}\left({G},Z\right)\leq 0$, and vanishes identically $F_{luc}\left({G},Z\right)\equiv 0$ for the exact solutions satisfying the PDE, the boundary and interface conditions. By Lemma 5 the boundary terms are negative semi-definite $\mathrm{BTs}\left(\widehat{v}_{\eta}^{-},\widehat{T}_{\eta}^{-}\right)\leq 0$.

We can now prove the theorem.

Theorem 7**.**

The elemental weak form (5.3) satisfies the energy equation

[TABLE]

Proof.

In (5.3) replace the test function $\boldsymbol{\phi}$ by the solution $\mathbf{Q}$, and integrate the conservative flux term only by parts. The volume term vanishes, having

[TABLE]

Using (79) in the right hand side of (86) gives

[TABLE]

Collecting contributions from the two elements on both sides of the interface, and using the identities (65) and (76) gives the energy equation (85). ∎

Remark 2**.**

It is particularly noteworthy that we have completely avoided explicit eigen-decomposition of the compound coefficient matrices

[TABLE]

which will be necessary for some DG methods [6, 24] based on classical Godunov numerical flux [28]. Here $\mathbf{P}$ is the parameter matrix defined in (10), ${\mathbf{A}}_{\xi}$ are the non-dimensional coefficient matrices defined in (9) and $\mathbf{n}=[n_{x},n_{y},n_{z}]^{T}$ is the outward unit normal on the boundary defined in (34). The eigen-decomposition of the matrices $\widetilde{\mathbf{A}}_{n}$ on each element face and for every element can be cumbersome in general curvilinear elements, and in particular in heterogeneous anisotropic elastic medium.

In the next section, we will introduce the discontinuous Galerkin approximation and prove numerical stability by deriving discrete energy estimates analogous to (85).

6 The discontinuous Galerkin approximation

Inside the transformed element $(q,r,s)\in\widetilde{\Omega}=[0,1]^{3}$, approximate the elemental solution and the conservative flux term by polynomial interpolants of degree $P$, and write

[TABLE]

where $\mathbf{F}_{ijk}\left(\bar{\mathbf{Q}}(t)\right)=\mathbf{F}\left(\bar{\mathbf{Q}}_{ijk}(t)\right)$, $\bar{\mathbf{Q}}_{ijk}(t)$ are the evolving elemental degrees of freedom to be determined, and $\boldsymbol{\phi}_{ijk}(q,r,s)$ are the $ijk$-th interpolating polynomial. We consider tensor products of nodal basis with $\boldsymbol{\phi}_{ijk}(q,r,s)=\mathscr{L}_{i}(q)\mathscr{L}_{j}(r)\mathscr{L}_{k}(s)$, where $\mathscr{L}_{i}(q)$, $\mathscr{L}_{j}(r)$, $\mathscr{L}_{k}(s)$, are one dimensional nodal interpolating Lagrange polynomials of degree $P$, with

[TABLE]

The interpolating nodes $q_{m}$, $m=1,2,\dots,P+1$, are the nodes of a Gauss quadrature with

[TABLE]

where $h_{i}>0$, $h_{j}>0$, $h_{k}>0$, are the quadrature weights. We will only use quadrature rules such that for all polynomial integrand $f(q)$ of degree $\leq 2P-1$, the corresponding one space dimensional rule is exact, $\sum_{i=1}^{P+1}f(q_{i})h_{i}=\int_{0}^{1}f(q)dq.$ Admissible candidates are the Gauss-Legendre-Lobatto quadrature rule with GLL nodes, the Gauss-Legendre quadrature rule with GL nodes and the Gauss-Legendre-Radau quadrature rule with GLR nodes. While both endpoints, $q=0,1$, are part of GLL quadrature nodes, the GLR quadrature contains only the first endpoint $q=0$ as a node. Lastly, for the GL quadrature, both endpoints, $q=0,1$, are not quadrature nodes. Note that when an endpoint is not a quadrature node, $q_{1}\neq 0$ or $q_{P+1}\neq 1$, extrapolation is needed to compute numerical fluxes at the element boundaries, $q=0,1$. We also remark that the GLL quadrature rule is exact for polynomial integrand of degree $2P-1$, GLR quadrature rule is exact for polynomial integrand of degree $2P$, and GL quadrature rule is exact for polynomial integrand of degree $2P+1$.

6.1 The spectral difference approximation

Introduce the square matrices $H,A\in\mathbb{R}^{(P+1)\times(P+1)}$, defined by

[TABLE]

Note that the matrix

[TABLE]

is a one space dimensional spectral difference approximation of the first derivative.

Using the fact that the quadrature rule is exact for all polynomial integrand of degree $\leq 2P-1$ implies that

[TABLE]

Equations (91) and (92) are the discrete equivalence of the integration-by-parts property. If both endpoints $q=0,1$ are quadrature nodes and we consider nodal bases then we have $B=\text{diag}[-1,0,0,\dots,0,1].$ The matrix $B$ projects the nodal degrees of freedom to element faces.

The one space dimensional derivative operator (91) can be extended to higher space dimensions using the Kronecker products $\otimes$, having

[TABLE]

[TABLE]

Here, $I$ is the $(P+1)\times(P+1)$ identity matrix, and $I_{9}$ is the $9\times 9$ identity matrix. Note that since $\mathbf{H}_{q},\mathbf{H}_{r},\mathbf{H}_{s}$ are diagonal matrices then $\mathbf{H}$ is also a diagonal matrix. The matrix $\mathbf{H}$ is independent of the order of the matrix products, that is, $\mathbf{H}=\mathbf{H}_{q}\mathbf{H}_{r}\mathbf{H}_{s}=\mathbf{H}_{r}\mathbf{H}_{q}\mathbf{H}_{s}=\mathbf{H}_{s}\mathbf{H}_{r}\mathbf{H}_{q}$.

We also introduce the projection matrices

[TABLE]

where

[TABLE]

6.2 The semi–discrete approximation

We will now make a classical Galerkin approximation, by choosing test functions in the same space as the basis functions. Thus, replacing ${\mathbf{Q}}(q,r,s,t)$ by $\bar{\mathbf{Q}}(q,r,s,t)$ and $\mathbf{F}\left({\mathbf{Q}}\right)$ by $\mathbf{F}\left(\bar{\mathbf{Q}}\right)$ in (5.3), and approximating all integrals with the corresponding quadrature rules yields the semi-discrete equation,

[TABLE]

for the evolving degrees of freedom, $\bar{\mathbf{Q}}(t)=[\bar{\mathbf{Q}}_{ijk}(t)]$. The numerical flux fluctuation term $\mathbf{Flux}\left(\bar{\mathbf{Q}}\right)$ implements the boundary conditions (50), and the interface conditions (55), at the element faces, and it is defined by

[TABLE]

with

[TABLE]

The elemental degrees of freedom have been arranged row-wise as a single vector of length $9(P+1)^{3}$. Note the close similarity between the semi-discrete approximation (95) and the continuous analogue (32). The discrete operator $\gradient_{D}=\left(\mathbf{D}_{q},\mathbf{D}_{r},\mathbf{D}_{s}\right)^{T}$ is also analogous to the continuous gradient operator $\gradient=\left(\partial/\partial q,\partial/\partial r,\partial/\partial s\right)^{T}$. In $\gradient_{D}$ we have replaced the continuous derivative operators in $\gradient$ with their discrete counterparts,

[TABLE]

where the spatial derivative operators, $\mathbf{D}_{q},\mathbf{D}_{r},\mathbf{D}_{s}$, are given in (93).

6.3 Numerical stability

To prove the stability of the semi-discrete approximation (95), we will derive a semi–discrete energy estimate analogous to (85). To begin, approximate the continuous energy in each element by the quadrature rule (89), having

[TABLE]

We also approximate the surface integrals in the boundary and interface terms (56), (5) and (5.3). We will now prove the discrete anti-symmetric property (22) and a discrete analogue of Lemma 3. We have

Lemma 6**.**

Consider the semi-discrete approximations of the transformed equation of motion (95), in curvilinear coordinates, with the flux terms and non-conservative product terms given by (33). For a polynomial approximation of degree $P$, if the quadrature rule is exact for all polynomial integrand $f(\xi)$ of degree $\leq 2P-1$, $\sum_{m=1}^{P+1}f(\xi_{m})h_{m}=\int_{0}^{1}f(\xi)d\xi$, then the corresponding spatial discrete operators satisfy the discrete anti-symmetric property

[TABLE]

and

[TABLE]

Proof.

Equation (99) follows directly from our choice of the anti-symmetric form (20). Now consider

[TABLE]

and use the discrete integrate-by-parts principle, (91) and (92), we have

[TABLE]

Using the fact that

[TABLE]

completes the proof. ∎

A direct consequence of the discrete anti-symmetric form (99) is the discrete identity

[TABLE]

For the two elements model we introduce the external boundary terms

[TABLE]

the interface term

[TABLE]

and the fluctuation term

[TABLE]

Remark 3**.**

As in the continuous setting, it is noteworthy that by construction the interface terms vanish identically $\mathscr{IT}_{s}(\widehat{v}^{\pm},\widehat{T}^{\pm})\equiv 0$, by inspection the fluctuation terms are negative definite $\mathscr{F}_{luc}\left({G},Z\right)<0$, and will vanish in the limit of mesh resolution. By Lemma 5 the boundary terms are negative semi-definite $\mathscr{BT}_{s}\left(\widehat{v}_{\eta}^{-},\widehat{T}_{\eta}^{-}\right)\leq 0$.

We can now prove the numerical stability of the semi-discrete approximation (95).

Theorem 8**.**

*Consider the semi-discrete DG approximation (95). If the quadrature rule $\sum_{m=1}^{P+1}f(\xi_{m})h_{m}=\int_{0}^{1}f(\xi)d\xi$ is exact for all polynomial integrand $f(\xi)$ of degree $\leq 2P-1$, then the numerical solution satisfies the energy equation

[TABLE]

Proof.

We will use the energy method. From the left, multiply (95) by $\bar{\mathbf{Q}}^{T}{\mathbf{H}}$

[TABLE]

On the left hand side of (105) we recognize time derivative of the semi-discrete energy $\mathscr{E}(t)$ defined in (98). With $\sum_{m=1}^{P+1}f(q_{m})h_{m}=\int_{0}^{1}f(q)dq$, then the summation-by-parts principle, (92), (91) holds. On the right hand side of (105), we use Lemma 6, that is replace the conservative flux term with (100), we have

[TABLE]

Using the fact (103), the volume terms vanish, remaining only the surface terms. Simplifying the flux term $\bar{\mathbf{Q}}^{T}{\mathbf{H}}\mathbf{Flux}\left(\bar{\mathbf{Q}}\right)$ with (97) having

[TABLE]

Using (79) in the right hand side of (107) gives

[TABLE]

Collecting contributions from both elements, and using the identities (65) and (76) give the energy equation (104). ∎

Note the similarities between Theorem 8 and Theorem 7. The only differences are that we have replaced the solution $\mathbf{Q}$ by its polynomial interpolant $\bar{\mathbf{Q}}$, and the integrals by quadrature rules.

By (104), the semi-discrete DG approximation (95) is asymptotically stable. This means that the solutions can never grow in time. We can now integrate (95) in time using any suitable explicit ODE scheme, such as high order explicit Runge-Kutta methods. In this paper we will integrate the (95) in time using the ADER scheme [42, 6]. The ADER time discretisation is summarized in A.

7 Implementation in the ExaHyPE engine: ExaSeis

We implement the here presented curvilinear method in the ExaHyPE-Engine [17] which is publicly available at www.exahype.org. The ExaHyPE-Engine solves linear and non-linear hyperbolic partial differential equations in first order formulation. In the linear case, ExaHyPE uses the ADER-DG method as summarised in A. To solve for any given hyperbolic PDE in ExaHyPE, users are required to implement all PDE specific terms in a generated C++ interface. For the method we present in this paper, these are the conservative and non-conservative flux terms of Eqs. (22) and (23) and the eigenvalues in Eq. (18). While the engine provides the Rusanov flux as default numerical flux, we here replace it with the numerical flux introduced in (96), to properly treat boundary conditions and material interfaces. As initial conditions for the experiments in Sec. 8 we implement problem specific material parameters and moment-tensor point sources. We generate curvilinear meshes based on a simple k-d-tree approach, which is publicly available as extracted library at https://gitlab.lrz.de/ExaHyPE-Seismic/curvi.

The elastic wave propagation model with curvilinear meshes is part of the ExaSeis application collection (https://gitlab.lrz.de/ExaHyPE-Seismic). ExaSeis includes several approaches related to computational seismology, allowing to study topography scattering effects or multi-physics dynamic rupture models. It also allows to use perfectly matching layers (PML) [47]. To allow reproducibility of all results of Sec. 8 we provide a repository with prepared scripts at https://gitlab.lrz.de/ExaHyPE-Seismic/ExaSeis-Benchmarks.

8 Numerical experiments

We will now present numerical simulations in complex 3D isotropic and anisotropic elastic solids. The experiments are designed to verify accuracy and numerical stability of the method, for both, body waves and elastic surface waves, as well as for high frequency scattered waves. We will also demonstrate the potential of our method for applications including geometrically complex free surface topography.

We will consider three benchmark problems with different levels of difficulty [43, 44, 56, 57, 45, 46, 58]. The benchmark problems are i) the 3D layer over a homogeneous half-space (LOH1) problem, ii) 3D elastic surface waves in an anisotropic medium, and iii) a 3D regional strong topography contrast simulation including scattered high frequency waves and a geologically informed complex surface geometry (Zugspitze model). These benchmark problems are designed to quantify and assess the accuracy of simulation codes for seismic surface and interface waves, and effective resolution of high frequency wave modes generated by scattering from the complex non-planar topography.

The numerical experiments will be conducted in 3D elastic media. We demonstrate spectral convergence by performing simulations on moderately fine 3D meshes where we vary the degree of the polynomial approximations, $P=3,5,7$. The simulated solutions are compared to analytical and reference solutions of community defined benchmark problems [43, 44, 45, 46]. We will perform error analysis and compute error parameters relevant to computational seismology. The reader is also referred to [33], where detailed grid convergence studies for the physics based flux for elastodynamics are performed in 1D and 2D.

In many of the simulations, we have used PML boundary conditions to efficiently absorb outgoing waves. The PML modelling error is set to $\approx 1\%$, which is a moderately small error tolerance. Details of the discretization and implementation of the PML is reported in [47]. We have used GLL nodes for most of the experiments presented here. However, it is important to note that equivalent results have been obtained using GL nodes.

We introduce the effective sub-element resolution, grid size $h=\Delta{x}/(P+1)$, comparable to a finite difference grid size. We will use the global time-step

[TABLE]

where $P$ is the degree of the polynomial approximation, $d=3$ is the spatial dimension, $\mathrm{CFL}=0.9$ and

[TABLE]

Note that $\sum_{\xi=x,y,z}n_{\xi}^{2}=1$, and in an isotropic medium the effective normal p-wave speed is $c_{n}=c_{p}$.

We will consider seismic sources defined by the singular moment tensor point source

[TABLE]

Here, $\mathbf{M}$ is the symmetric second order moment tensor, $\delta_{\eta}(\eta)$ are the one dimensional Dirac delta function, $(x_{0},y_{0},z_{0})$ is the source location, and $g(t)$ is the source time function.

8.1 Layer over a half-space (LOH1)

We consider the 3D LOH1 benchmark problem, which has an analytical solution [43, 44]. The LOH1 benchmark problem consists of a planar free surface and an internal planar interface separating a thin homogeneous soft layer and hard half-space. The material properties in the medium are given by

[TABLE]

The benchmark considers homogeneous initial conditions for all fields, and generates waves by adding the double-couple moment tensor point source at depth, x=2\leavevmode\nobreak\km below the free surface. The source time function is given by

[TABLE]

and the moment tensor $\mathbf{M}$ which is zero except for the shear components $M_{yz}=M_{zy}=M_{0}$, and $M_{0}=10^{18}$ Nm is the moment magnitude. We place 9 receivers on the free-surface, at $x=0$, where the solutions are sampled. Table 1 shows the positions of the receivers relative to the epicenter. Receiver 6 and Receiver 9, which are about 10.39\leavevmode\nobreak\km away from the epicenter, are the two farthest receivers from the source.

Note that the domain is unbounded at depth and in the tangential directions. To perform numerical simulations we consider the bounded computational cube $(x,y,z)=[0,16.333\leavevmode\nobreak\ \text{km}]\times[-2.287\leavevmode\nobreak\ \text{km},14.046\leavevmode\nobreak\ \text{km}]\times[-2.287\leavevmode\nobreak\ \text{km},14.046\leavevmode\nobreak\ \text{km}]$, and the source is located at $(x_{0},y_{0},z_{0})=(2\leavevmode\nobreak\ \text{km},0,0)$. The truncated boundaries of the domain were surrounding by the PML to absorb outgoing waves. First we discretize the computational domain uniformly with 25 DG elements in each spatial direction. We consider degrees $P=3,5,7$ polynomial approximations, and run the simulation until the final time t=9\leavevmode\nobreak\s.

The solutions are displayed in Figure 4 for Receiver 6 and in Figure 5 for Receiver 9, and for $P=3,5,7$. With increasing polynomial degree, $P=3,5,7$, the numerical solutions converge spectrally to the analytical solution. For $P=5,7$, the numerical solutions match the exact solution excellently well, with about $\sim 1$% relative error.

8.1.1 Error analysis

We now analyse the numerical error for the LOH1 problem, and compute error parameters relevant to seismological applications. To quantitatively assess the accuracy of the numerically simulated seismograms, we compare the numerical seismograms for the LOH1 benchmark simulations with the exact reference solution using the time-frequency (TF) misfit criteria proposed by [43, 44]. We will briefly describe the technique here, and refer the reader to [43, 44] for a more elaborate discussion. The misfit criteria are based on the Time-Frequency Representation (TFR) of the seismogram, denoted as $W(t,f)$. A local TF envelope difference is defined as:

[TABLE]

and local TF phase difference is defined as

[TABLE]

where $|W(t,f)|$ and $|Wref(t,f)|$ are the TFR of the numerical and reference data, respectively. To obtain a single-valued measure of the EM or PM between these two seismograms, we use the following normalized formulars

[TABLE]

We use the software package Obspy [59] to post process the seismograms, by passing the seismograms through a band limited filter, and evaluate the accuracy level of the numerical seismograms at all stations.

Table 2 classifies the quality of the numerical seismograms, where $A$, with EM, PM $\leq 5\%$, indicates the highest quality numerical seismograms [43, 44].

The results of the analysis, for Receiver 6 with $P=3,5,7$, are displayed in Figure 6.

As shown in Figure 6 and Table 3, at this resolution the seismograms for polynomial approximations of degree $P=5,7$, belong to the A-class, are of the highest quality. We can also see that the seismograms for polynomial approximations of degree $P=3$ belong to the C-class. For $P=3$, the accuracy can be improved by increasing the mesh resolution (h-refinement).

8.1.2 Non-conforming adaptive mesh

Further, we investigate static adaptive mesh refinement, and stability for non conforming elements. Now the domain is discretized adaptively with uniform elements of $\Delta{x}=17/9$ km in the bedrock when $x>4$ km, $\Delta{x}=17/(9\times 3)$ km in the source region when $1<x\leq 4$ km, and $\Delta{x}=17/(9\times 9)$ km in the topmost layer $0\leq x\leq 1$ km (see the left panel in Figure 7). Note that in the topmost layer, closest to the surface, where surface waves are present the mesh is 3 times finer. We used GL nodes and degree $P=3$ polynomial approximation.

The effective grid resolution comparable to a finite difference grid size is $h=53$ m in the topmost layer $0\leq x\leq 1$ km, $h=3\times 53$ m in the source region $1<x\leq 4$ km, $h=9\times 53$ m in the bedrock $x>4$ km. We evolve the wave fields with the time-step $\Delta{t}=0.0026$ s until the final time of 9 s.

As above, we post-process the seismograms, and compute the EM and the PM. For Receiver 3 and the $z$ component of the particle velocity $v_{z}$, these quantities are displayed in the right panel of Figure 7. Note that EM$\sim 5\%$, and PM$\sim 2\%$, which are reasonably moderate error tolerances.

8.2 Surface waves in anisotropic media

We consider simulations of elastic waves propagating in an orthotropic anisotropic medium [45, 46]. Wave propagation in anisotropic elastic media is more complex than in isotropic elastic media. From both analytical and numerical perspectives, anisotropy introduces additional challenges, since a wave mode can propagate with different wave speeds in different directions. In fact, analytical solutions exist for only few problems. Numerical simulations become inevitable for studying wave propagation in anisotropic elastic media. We will study the accuracy of the method for surface waves in an anisotropic crystal, where analytical solutions exists. The crystal is apatite and the set-up is the same as in [46]. The density of the medium and elastic constants are given in Table 4, with $c_{55}=c_{44}$, $c_{66}=\left(c_{11}-c_{12}\right)/2$.

From (18), we can determine the wave speeds

[TABLE]

The computational domain is the cube $\left(x,y,z\right)\in[0,20\leavevmode\nobreak\ \text{cm}]\times[0,20\leavevmode\nobreak\ \text{cm}]\times[0,20\leavevmode\nobreak\ \text{cm}]$. Waves are excited by adding the point source

[TABLE]

to the first component of the momentum equation, where the source time function is

[TABLE]

$\delta_{\eta}(\eta-\eta_{0})$ are the one dimensional Dirac delta function. The source is located at $x_{0}=10\leavevmode\nobreak\ \text{cm}$, $y_{0}=10\leavevmode\nobreak\ \text{cm}$, $z_{0}=0$, the surface of the crystal. At the surface $z=0$, we impose the free surface boundary condition and the absorbing boundary condition at the other 5 boundaries of the domain.

We discretize the computational domain uniformly with 27 elements in each spatial direction, with a polynomial approximation of degree $P=5$. We place a receiver at $x_{0},=10\leavevmode\nobreak\ \text{cm}$, $y_{0},=10\leavevmode\nobreak\ \text{cm}$, $z_{0}=15\leavevmode\nobreak\ \text{cm}$, and evolve the wave fields, with the global time step $\Delta{t}=0.015\leavevmode\nobreak\ \mu$s, until the final time, $t=50\leavevmode\nobreak\ \mu$s. In Figure 8, we compare the analytical solution with the numerically derived solution.

The numerical solution matches the analytical solution very well, with less than $1\%$ relative error.

8.3 Wave propagation in complex geometries

We will now demonstrate the potential of the method in modeling elastic wave propagation in a 3D domain with geometrically complex free surface topography. We consider a 3D setup incorporating the strong topography contrasts of Mount Zugspitze, Germany.

Mount Zugspitze is reaching about $3$ km above sea level and is the highest peak of the Wetterstein mountain range as well as the highest mountain in Germany. The topography of this region is complex, with large deviations from planarity across the Earth’s surface, and constitute a significant challenge for numerical solvers. Effective numerical resolution of the scattered wavefield is very important for application purposes, and for example may affect the resolution of regional seismic tomography performed by the AlpArray project [60]. In particular, due to the high frequency wave modes generated by scattering from the complex non-planar topography, accurate and stable numerical simulation of seismic wave propagation in this region is both numerically and computationally arduous. The Zugspitze model has no analytical solutions. We will verify accuracy by making comparisons with the reference data in [58] simulated with WaveQLab3D [41], a petascale finite difference elastic wave solver.

The topography data spans a rectangular surface area $(x,z)\in[0,80\leavevmode\nobreak\ \text{km}]\times[0,80\leavevmode\nobreak\ \text{km}]$, with 100\leavevmode\nobreak\m resolution of the elevation data. Down dip we truncate the domain at y=80\leavevmode\nobreak\ \km, so that the computational domain corresponds to the modulated cuboid $(x,y,z)\in[0,80\leavevmode\nobreak\ \text{km}]\times[\widetilde{y},80\leavevmode\nobreak\ \text{km}]\times[0,80\leavevmode\nobreak\ \text{km}]$, where $\widetilde{y}(x,z)$ is the elevation data. We have processed the topography data by passing it through a band limited filter, see Figure 9, obtaining 3 realisations: Realisation A, Realisation B and Realisation C. Realisation A is a flat topography and includes low frequencies, Realisation B includes intermediate frequencies and Realisation C includes high frequencies.

At each truncated boundary in $x$-axis and $z$-axis, and down dip at $y=80$ m, we have included a 5 km absorbing layer where PML boundary conditions [47] are implemented to prevent artificial reflections from the computational boundaries from contaminating the solution. A stable numerical implementation of the PML for 3D linear elastodynamics is nontrivial and allows the generation of high quality seismograms.

We consider the homogenous material properties

[TABLE]

The domain is discretized uniformly with $241$ elements in each spatial direction and degree $P=3$ polynomial approximation. The effective sub-cell grid resolution comparable to a finite difference grid size is h=\Delta{x}/(P+1)=93.36\leavevmode\nobreak\m. We evolve the solution with the time step (109) until the final time $t=30$ s.

Snapshots of numerical simulations are shown in Figure 10 and the seismograms are shown in Figures 11–13, for the 3 receiver stations and the 3 topography realisations. For a low frequency topography model, Realisation A, note that the seismograms are dominated by the direct arrivals and waves propagate coherently. The seismograms for Realisations B and C are dominated by high frequency scattered waves, which are particularly present in the coda waves. We compare seismograms from ExaSeis and WaveQLab3D with $h=100$ m finite difference grid spacing. We observe a near perfect agreement of the seismograms at all frequencies, see Figures 11–13. However at high frequencies there are tiny differences, which will diminish as we increase the resolution.

9 Summary and outlook

In this paper, we present a new energy-stable discontinuous Galerkin approximation of the elastic wave equation in general and geometrically complex 3D media, using the physics-based numerical penalty-flux [33]. As opposed to the Godunov flux, the physics-based flux does not require a complete eigenvector-eigenvalue decomposition of the spatial coefficient matrices. By construction, our numerical flux is upwind and yields a discrete energy estimate analogous to the continuous energy estimate. The discrete energy estimate holds for conforming and non-conforming curvilinear elements. The ability to handle non-conforming curvilinear meshes allows for flexible adaptive mesh refinement strategies. The numerical scheme has been implemented in the ExaSeis application in ExaHyPE [17], a simulation engine for hyperbolic PDEs on adaptive structured meshes. Numerical experiments are presented in 3D isotropic and anisotropic media demonstrating stability and accuracy. Finally, we present numerical verification in a regional geophysical wave propagation problem in an Earth model with geometrically complex free-surface topography.

We utilize a recent extension of the method to discretize PML boundary conditions for elastodynamics [47]. We expect, that the proposed numerical method will extend to model linear and nonlinear boundary and interface wave phenomena. In a forthcoming paper we will extend the method to nonlinear friction problems, and present numerical simulations of nonlinear dynamic earthquake ruptures on dynamically adaptive meshes, embedded in 3D geometrically complex solid earth models.

Acknowledgments

The work presented in this paper was enabled by funding from the European Union’s Horizon 2020 research and innovation program under grant agreements no. 671698 (ExaHyPE), no. 823844 (ChEESE) and no. 852992 (TEAR). The authors also acknowledge support by the German Research Foundation (DFG) (grants no. GA 2465/2-1, GA 2465/3-1), by KAUST-CRG (grant no. ORS-2017-CRG6 3389.02) and by KONWIHR (project NewWave). A.L. is supported by the Swiss Federal Institute of Technology grant (project ETH-10 17-2). Computing resources were provided by the Institute of Geophysics of LMU Munich [61], the Leibniz Supercomputing Centre (SuperMUC-NG project pr63qo) and the KAUST Shaheen Supercomputing Laboratory (project k1488).

The first author KD would like to thank Dimitri Komatitsch (1970–2019) for his help, fruitful discussions and for providing the analytical solution for elastic surface waves in an anisotropic medium.

Appendix A The Arbitrary DERivative (ADER) time integration

In this section we will summarize the ADER time-stepping scheme. For more elaborate discussions, we refer the reader to [42, 6, 23]. To begin, we rewrite the semi-discrete approximation (95) as a system first order ODEs

[TABLE]

with

[TABLE]

where $D\bar{\mathbf{Q}}$ is the discrete spatial operator is split into the derivative term, emanating from the PDE, and $F\bar{\mathbf{Q}}$ is the numerical flux fluctuation term, incorporating the boundary and interface conditions.

The numerical flux fluctuation is a very small term, $F\bar{\mathbf{Q}}\approx 0$, and will vanish $F\bar{\mathbf{Q}}\to 0$ in the limit of mesh refinement $\Delta{t}\to 0$.

We now introduce the discrete time variables $t_{k}\leq t\leq t_{k+1}$, $\Delta{t}_{k}=t_{k+1}-t_{k}$, and the pseudo time variable $\tau=t-t_{k}$ such that $0\leq\tau\leq\Delta{t}_{k}$, and ${d}/{d\tau}={d}/{dt}$. Going from the current time level $\tau=0$ to the next time level $\tau=\Delta{t}_{k}$, we integrate the ODE (113), exactly having

[TABLE]

where the second equality follows from linearity. If we can evaluate the integrals $\int_{0}^{\Delta{t}_{k}}{\bar{\mathbf{Q}}}d\tau$ in (115) exactly, then the time integration in (115) is exact. However, exact time integration is possible only in the most trivial case where the right hand side of (113) vanish identically for all components. Now, we will make an important approximation. We assume that the time step $\Delta{t}$ is sufficiently small, such that

[TABLE]

are reasonable approximations. Next we construct the predictor, $\widetilde{\bar{\mathbf{Q}}}(\tau)$, by Taylor expansions of the solution around $\tau=0$ and replace the time derivatives with spatial operator in (116), we have

[TABLE]

where $P$ is the polynomial degree used in the spatial approximation. We can now approximate the integrals in (115) using the predictor. The result of this integration is called the time average, $\bar{\bar{\mathbf{Q}}}(0)$,

[TABLE]

By replacing the integrals in (115) with the time average ${\bar{\bar{\mathbf{Q}}}(0)}$, we derive a high order accurate, explicit, one-step, time integration scheme

[TABLE]

Note that the numerical flux fluctuations, $F\bar{\bar{\mathbf{Q}}}(0)$, is evaluated only once for any order of approximation. This is opposed to Runge–Kutta methods or standard Taylor series methods where the numerical flux fluctuation is included in the spatial operator, $D^{m}\to(D+{F})^{m}$, to approximate higher time derivatives in the Taylor series terms. For the ADER scheme, the fact that the numerical flux fluctuation is evaluated only once for any order implies that most of the computations are performed within the element to compute the predictor in (117) and the time average in (118). This has a huge impact in high performance computing applications, since we can design efficient communication avoiding parallel algorithms. Since the predictor $\widetilde{\bar{\mathbf{Q}}}(\tau)$ is defined in the entire time interval $0\leq\tau\leq\Delta{t}_{k}$, the ADER scheme, (115)–(119), is also easily amenable to local time-stepping methods. When the ADER time stepping scheme is combined with the DG spatial approximation the fully discrete scheme is often referred as the ADERDG method [6, 23]. For a DG polynomial approximation of degree $P$, a stable ADERDG method is $(P+1)$th order accurate in both space and time.