# A spectral hybridizable discontinuous Galerkin method for elastic-acoustic wave propagation

S. Terrana, Jean-Pierre Vilotte, L. Guillot

<figure><img src="image_1.png" /></figure>

## To cite this version:

S. Terrana, Jean-Pierre Vilotte, L. Guillot. A spectral hybridizable discontinuous Galerkin method for elastic-acoustic wave propagation. Geophysical Journal International, 2018, 213, pp.574-602. 10.1093/gji/ggx557. insu-03589361

HAL Id: insu-03589361

https://insu.hal.science/insu-03589361

Submitted on 25 Feb 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# A spectral hybridizable discontinuous Galerkin method for elastic-acoustic wave propagation

S. Terrana,<sup>1,2</sup> J.P. Vilotte<sup>1</sup> and L. Guillot<sup>2</sup>

Accepted 2017 December 20. Received 2017 December 4; in original form 2017 August 2

## SUMMARY

We introduce a time-domain, high-order in space, hybridizable discontinuous Galerkin (DG) spectral element method (HDG-SEM) for wave equations in coupled elastic-acoustic media. The method is based on a first-order hyperbolic velocity-strain formulation of the wave equations written in conservative form. This method follows the HDG approach by introducing a *hybrid* unknown, which is the approximation of the velocity on the element boundaries, as the only globally (i.e. interelement) coupled degrees of freedom. In this paper, we first present a hybridized formulation of the exact Riemann solver at the element boundaries, taking into account elastic-elastic, acoustic-acoustic and elastic-acoustic interfaces. We then use this Riemann solver to derive an explicit construction of the HDG stabilization function $\tau$ for all the above-mentioned interfaces. We thus obtain an HDG scheme for coupled elastic-acoustic problems. This scheme is then discretized in space on quadrangular/hexahedral meshes using arbitrary high-order polynomial basis for both volumetric and hybrid fields, using an approach similar to the spectral element methods. This leads to a semi-discrete system of algebraic differential equations (ADEs), which thanks to the structure of the global conservativity condition can be reformulated easily as a classical system of first-order ordinary differential equations in time, allowing the use of classical explicit or implicit time integration schemes. When an explicit time scheme is used, the HDG method can be seen as a reformulation of a DG with upwind fluxes. The introduction of the velocity hybrid unknown leads to relatively simple computations at the element boundaries which, in turn, makes the HDG approach competitive with the DG-upwind methods. Extensive numerical results are provided to illustrate and assess the accuracy and convergence properties of this HDG-SEM. The approximate velocity is shown to converge with the optimal order of $k+1$ in the $L^2$-norm, when element polynomials of order $k$ are used, and to exhibit the classical spectral convergence of SEM. Additional inexpensive local post-processing in both the elastic and the acoustic case allow to achieve higher convergence orders. The HDG scheme provides a natural framework for coupling classical, continuous Galerkin SEM with HDG-SEM in the same simulation, and it is shown numerically in this paper. As such, the proposed HDG-SEM can combine the efficiency of the continuous SEM with the flexibility of the HDG approaches. Finally, more complex numerical results, inspired from real geophysical applications, are presented to illustrate the capabilities of the method for wave propagation in heterogeneous elastic-acoustic media with complex geometries.

**Key words:** Numerical solutions; Computational seismology; Guided waves; Interface waves; Wave propagation.

## 1 INTRODUCTION

In geophysics, mechanical waves are an invaluable tool to get clues about the Earth’s internal structure, the internal seismic sources, and the dynamical coupling between the solid Earth and its external fluid envelopes (oceans and atmosphere). Today acquisition systems combine seismic receivers, hydrophones, and microbarometers providing observations on the elastic, hydro-acoustic and infrasound wave fields generated by both natural and anthropic sources.

Exploiting these observations requires numerical simulation methods for coupled elastic-acoustic wave propagation that are: (i) accurate enough for very long propagations (i.e. low numerical dissipation and dispersion); (ii) able to model geophysical media involving complex

<footnote n="1"><sup>1</sup>Institut de Physique du Globe de Paris, 1 Rue Jussieu, F-75005 Paris, France. E-mail: terrana@mit.edu</footnote>

<footnote n="2"><sup>2</sup>CEA, DAM, DIF, F-91297 Arpajon, France</footnote>
interfaces geometries; (iii) able to deal with strong impedance contrasts in particular at interfaces between solid Earth and the ocean or the atmosphere, with accurate fluid-solid transmission conditions. Classical semi-analytical or asymptotic methods (e.g. Chapman 2004) can be used in simple configurations, but when considering very complex media, purely numerical methods are needed.

With increasing capacity and capability of high performance computing architectures, numerical simulation of elastic-acoustic wave propagation has become feasible and extensively used in geophysics. Mathematical and geophysical literature contains a plethora of numerical methods and theoretical approximations for dealing with elastic and acoustic wave equations (Maupin & Wu 2007; Virieux et al. 2012; Moczo et al. 2014), each with its own problem-dependant strengths-and-weaknesses.

Continuous Galerkin Spectral Element Methods (CG-SEMs), originally developed in the context of fluid mechanics by Patera (1984) and Maday & Patera (1989), were extended to solve time domain wave equations (Seriani & Priolo 1994; Faccioli et al. 1997; Komatitsch & Vilotte 1998; Komatitsch et al. 1999) and since then have gained increasing interest in geophysics for a wide range of applications (Chaljub et al. 2003, 2007; De Martin 2011; Peter et al. 2011; Cupillard et al. 2012; Bottero et al. 2016; Tsuboi et al. 2016).

Classical CG-SEM makes use of an hexahedral mesh as support for the spatial discretization. On each hexahedral element, the solution is approximated by a high-order polynomial which is a tensor-product of 1-D Lagrange polynomials associated with Gauss-Lobatto-Legendre (GLL) nodes, while the $C^0$ continuity of the approximation is ensured on the whole computational domain. Choosing quadrature points coinciding with the interpolation nodes, and using an explicit second-order energy-conserving time scheme leads to a very efficient overall numerical scheme. As such, classical CG-SEM handles naturally heterogeneous media, achieves high-order convergence and spectral accuracy in space, minimal numerical dispersion for waves modelling, and can be successfully implemented on modern hybrid and multicore computing architectures (Tsuboi et al. 2016). Moreover, the classical CG-SEM has been extended to solve coupled elastic-acoustic problems, using a generalized potential formulation of the acoustic wave equations (Komatitsch et al. 2000; Chaljub et al. 2003; Chaljub & Valette 2004; Chaljub et al. 2007).

However, geometrical flexibility has long been recognized as an issue in classical CG-SEM. The two main challenges are the generation of hexahedral meshes for complex 3-D geometries, and the difficulty of handling non-conforming meshes. A solution to the first challenge has been an extension of classical CG-SEM to unstructured simplicial (i.e. tetrahedral) meshes, which has been actively developed by Cohen et al. (2001), Komatitsch et al. (2001), Mercerat et al. (2006) and Mazzieri & Rapetti (2012). With judicious choices of interpolation nodes and quadrature points, these approaches show to achieve dispersion characteristics similar to their hexahedral counterpart (Liu et al. 2012). However, practitioners have been wary using these methods for large scale problems due to the significant extra computational cost of tetrahedra-based approximations. An ideal approach would then be using hex-dominant meshes, whose generation is now well mastered (Meyers et al. 1998; Yamakawa & Shimada 2003), but the transition from hexahedra to tetrahedra would then be troublesome for a $C^0$-continuous approach.

This leads to the second challenge: whether small elements appear to capture complex geometrical features, or whether strong velocity contrasts are present (for instance at the crust-atmosphere interface), the mesh size and the polynomial order need to be adapted in different parts of the domain in a non-conform way, leading to the so-called *h*-adaptivity and *p*-adaptivity strategies, respectively. These strategies are sometimes necessary both to provide an appropriate spatial sampling for each geophysical media, and to avoid the smallest elements to lower inadmissibly the time step used for an explicit time integration. It is possible to use these strategies in a CG-SEM framework, using a mortar element-like strategy (Chaljub et al. 2003) or a mortar-free approach (Diaz & Joly 2005). But in all cases, a global system related to the non-conforming surfaces needs to be inverted at every time step. This implies a significant increase in the implementation complexity and in the computational cost, and it may even affect the stability of the time integrator (Diaz & Joly 2005). Indeed it turns out to be rather unnatural for a continuous method, such as CG-SEM, to handle non-conforming meshes, which are discontinuous in essence.

That is one of the reasons why, over the past few years, considerable attention has been turned to time domain Discontinuous Galerkin (DG) methods, for the solution of the elastic-acoustic wave propagation. In contrast with CG methods, no interelement continuity is imposed to the solution, but elements exchange numerical fluxes. As a result, DG methods get some attractive properties. In particular, they are more flexible than CG methods for complex geometries, allow efficient *h/p* and time adaptivity strategy, and retain very good scalability properties.

Second-order DG formulations have been developed in the framework of the interior penalty formulations, (e.g. Grote et al. 2006; Chung & Engquist 2006; Ainsworth et al. 2006; de Basabe et al. 2008).

For first-order formulation of the wave equations, at least four approaches have been proposed and enjoy some popularity in geophysics due to their ability to handle elastic-acoustic coupling. The ADER-DG scheme, a high-order tetrahedra-based DG method has been developed (e.g. Käser & Dumbser 2006; de la Puente et al. 2007), based on a local modal approximation and upwind fluxes together with high-order Taylor expansions as explicit time integration schemes. This method has been proven to handle efficiently *h/p* and time adaptivity (Hermann et al. 2011; Dumbser et al. 2007). More recently low-order tetrahedra with centred fluxes, and using local nodal approximation has been developed (Delcourte et al. 2009; Etienne et al. 2010) and extended to higher orders by Mercerat & Glinsky (2015) and by Delcourte & Glinsky (2015) with high-order leap-frog time schemes. An original DG approach from Ye et al. (2016) proposes high-order tetrahedra with penalized upwind numerical fluxes and implicit-explicit time scheme. Finally, a high-order hexahedra DG method with upwind fluxes and a Runge-Kutta explicit time scheme was proposed by Wilcox et al. (2010) and Bui-Thanh & Ghattas (2012). The latter appears quite attractive due to its lower complexity compared to the modal approach (Hesthaven & Warburton 2008).

However in spite of all their nice properties DG methods have not yet made significant impact in geophysics, with some notable exceptions (Minisini et al. 2013; Wenk et al. 2013). The high computational cost and memory print are indeed major impediments to the widespread application of DG methods to real-world applications. Even though the generated linear system is more weakly coupled, DG
methods have more degrees of freedom for the same mesh and the same polynomial degree than CG methods—due to the nodal duplication at the element boundary interfaces—and more globally coupled degree of freedom since fluxes between neighbouring elements must be accounted for. Therefore, it would be desirable to develop new DG methods for elastic–acoustic wave equations in a framework that retains the nice properties of DG methods while being competitive with CG-SEMs, and eventually allowing efficient coupling with CG-SEM for practical applications.

Hybridizable DG (HDG) methods were introduced by Cockburn et al. (2009) for elliptic problems to address some of the above limitations and can be seen as the hybridization of early mixed methods (Raviart & Thomas 1977; Brezzi et al. 1985). The HDG method has two decisive advantages over the traditional CG and DG methods. First, in HDG methods a hybrid variable is introduced, which is the approximation of the solution on the element boundaries, that is, on the mesh skeleton only. This hybrid unknown is then the only global (i.e. interelement) coupling field. When solving for a steady-state problem, or when using an implicit time marching scheme, the globally coupled system involves only the hybrid degrees of freedom. It is then a significantly smaller global system when compared to CG or DG methods. Second, when polynomial degree $k \ge 1$ are used to represent the numerical solution, the approximate gradient computed with HDG has been shown to converge with an optimal order $k + 1$, allowing the solution to converge with an extra order $k + 2$ when a suitable local post-processing is used (Soon et al. 2009; Cockburn & Shi 2013; Cockburn & Stolarski 2014; Fu et al. 2015). These two attractive features have been powerful arguments for the development of HDG methods as solvers for various types of PDEs: diffusion problems (Chabaud & Cockburn 2012), incompressible flows (Cockburn et al. 2010), compressible flows (Peraire et al. 2010), Maxwell equations (Nguyen et al. 2011a), Helmholtz equations (Giorgiani et al. 2013), and finally acoustic and elastic wave equations in frequency domain (Bonnasse-Gahot 2015; Hungria et al. 2017) or in time domain with implicit time schemes (Cockburn & Quenneville-Bélair 2014; Nguyen et al. 2011b).

Recently, an HDG method with an explicit time scheme has been proposed for acoustic wave propagation by Stanglmeier et al. (2016), which turns out to be an explicit DG method with upwind fluxes for a particular choice of the stabilization function. The links between explicit DG and HDG have been also studied by Bui-Thanh (2015) for a wide class of problems. And the development in this paper can be seen as an extension of these previous works for coupled elastic–acoustic wave equations. However, when an explicit time integration is used there is no longer any global system to invert, and one may wonder if the HDG method is then still attractive. It turns out that it is. First, because the hybridization does not add any complexity to the computation of the fluxes, which makes HDG competitive with DG. Second, the post-processing allows HDG to obtain an extra accuracy at low cost. Third, Cockburn et al. (2009) have proven that the HDG can handle elegantly the h/p-adaptivities and can be straightforwardly coupled with CG. All these points make HDG attractive when compared with other DG methods.

This paper is organized as follows. A brief description of the first-order conservative formulation of the acoustic and elastic wave equations is outlined in Section 2. The hybrid velocities, tractions and pressures are introduced in Section 3 on simple two-domain decompositions. Then the transmission relations satisfied by these hybrid variables are presented for elastic, acoustic and elastic–acoustic domain decompositions. In Section 4, we outline new hybridized versions of exact Riemann solvers, in order to derive local versions of Dirichlet-to-Neumann (DtN) operators at the subdomains interfaces. Sections 3 and 4 thus provide the relations between the hybrid variables that are used to derive our HDG discrete formulation for coupled elastic–acoustic wave equations in Section 5. An interesting point is that our approach provides a physical way to derive the parameter $\tau$, usually interpreted as a numerical stability parameter in the HDG literature. A short analysis of the method in terms of local conservation and semi-discrete stability is then provided. The HDG scheme leads to a system of algebraic differential equations (ADEs) that can be rewritten as a first-order system of ordinary differential equations (ODEs). Therefore, the time integration can be performed with any classical explicit ODE solver (Section 6), such as midpoint or Runge–Kutta methods. Finally, the accuracy and the convergence of the explicit HDG method are analysed using a number of selected academic test problems in Section 7. The optimal convergence rate $k + 1$ is observed for the velocity, and a better precision is obtained when an extra post-processing step is applied. An accurate coupling with CG-SEM is outlined, and finally, the capabilities of the method are illustrated through more realistic elastic–acoustic problems.

## 2 ELASTIC AND ACOUSTIC WAVE EQUATIONS

First, we start with a description of the geometry and fix some notations that will be used throughout this paper. Let us consider a bounded physical domain, which is identified to the open subset $\Omega \subset \mathbb{R}^{n_d}$, with $n_d$ the spatial dimension. $\bar{\Omega}$ denotes the closed subset $\Omega \cup \Gamma$, where $\Gamma$ is the external boundary of $\Omega$. The physical domain is made of solid and fluid regions, which, respectively, are referred to as $\Omega_S$ and $\Omega_F$:

$$ \bar{\Omega} = \bar{\Omega}_S \cup \bar{\Omega}_F. \qquad (1) $$

The external boundary $\Gamma$ can be decomposed as $\Gamma = \Gamma^D \cup \Gamma^N \cup \Gamma^{abs}$. Here $\Gamma^D$ and $\Gamma^N$ are regions of $\Gamma$ where, respectively, Dirichlet and Neumann boundary conditions are applied. $\Gamma^{abs}$ denotes the boundary region where outgoing waves are absorbed.

### 2.1 Elastic wave equations

Let us consider the solid domain $\Omega_S$, and a time interval of interest $I = [0, T] \in \mathbb{R}_+$. The displacement field is denoted by $u(x, t) : \bar{\Omega}_S \times I \to \mathbb{R}^{n_d}$, and the velocity field by $v(x, t)$, while $x$ is the position at rest with respect to a reference coordinates system $(0, x_1, ..., x_{n_d})$.
The Euler-Lagrange equations of elastodynamics are obtained by a first-order adiabatic perturbation around a reference state of hydrostatic equilibrium without taking into account the effects of gravity or rotation, and can be written as a hyperbolic system of first-order partial differential equations in conservative form

$$ \frac{\partial \epsilon}{\partial t} = \nabla^2 \mathbf{v} \qquad \text{on } \Omega_S \times \mathbf{I}, \qquad (2a) $$

$$ \frac{\partial \rho \mathbf{v}}{\partial t} = \nabla \cdot \boldsymbol{\sigma}(\mathbf{u}) + \mathbf{f} \quad \text{on } \Omega_S \times \mathbf{I}, \qquad (2b) $$

where $\rho(\mathbf{x})$ is the reference mass density, and $\mathbf{f}(\mathbf{x}, t)$ is a generalized body force. In (2) $\nabla \mathbf{a}$, $\nabla^s \mathbf{a}$ and $\nabla \cdot \mathbf{a}$, respectively, denote the gradient, the symmetric part of the gradient, and the divergence of a given tensor field $\mathbf{a}$.

The symmetric Cauchy stress tensor is denoted by $\boldsymbol{\sigma} : \bar{\Omega}_S \times \mathbf{I} \to \mathbb{S} \subset \mathbb{R}^{n_d \times n_d}$, where $\mathbb{S}$ is the subspace of symmetric second-order tensors of dimension $n_d(n_d + 1)/2$, and is linearly related to the infinitesimal strain tensor $\epsilon \in \mathbb{S}$ through the Hooke's law:

$$ \boldsymbol{\sigma}(\mathbf{x}, t) = \mathbf{C}(\mathbf{x}) : \boldsymbol{\epsilon}(\mathbf{x}, t) \quad \forall \mathbf{x} \in \bar{\Omega}_S, \forall t \in \mathbf{I}, \qquad (3) $$

where the fourth-order stiffness tensor $\mathbf{C}$ controls the speed of propagating waves in the medium and the symbol ‘:’ stands for the tensor double product. Then, because the pre-stress is hydrostatic, $\mathbf{C}$ is symmetric, with minor and major symmetries, and positive definite restricted to $\mathbb{S}$. Furthermore, assuming an isotropic body, $\mathbf{C}$ can be fully represented by the two Lamé parameters $\lambda$ and $\mu$

$$ C_{ijkl} = \lambda \delta_{ij} + \mu (\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}), \qquad (4) $$

where $\delta$ is the Kronecker delta.

The infinitesimal strain tensor is defined as,

$$ \boldsymbol{\epsilon} = \nabla^s \mathbf{u} = \frac{1}{2} (\nabla \mathbf{u} + \nabla^T \mathbf{u}) \quad \text{on } \Omega_S \times \mathbf{I}, \qquad (5) $$

The system of first-order partial differential equations (2) can be rewritten in a compact and conservative form as,

$$ \frac{\partial \mathbf{q}}{\partial t} + \nabla \cdot \mathcal{F} \mathbf{q} = \mathbf{f}, \qquad (6) $$

where the **q**-vector, of dimension $m_d$, is associated with the state variables. When $n_d = 3$ we have

$$ \mathbf{q} = \begin{pmatrix} \boldsymbol{\epsilon} \\ \rho \mathbf{v} \end{pmatrix} = \left( \epsilon_{11} \ \epsilon_{22} \ \epsilon_{33} \ \epsilon_{13} \ \epsilon_{12} \ \rho v_1 \ \rho v_2 \ \rho v_3 \right)^T. \qquad (7) $$

$\mathcal{F} \mathbf{q}$ is then the associated $m_d \times n_d$ flux matrix. Thanks to the linearity of the hyperbolic problem, (6) can be reformulated as

$$ \frac{\partial \mathbf{q}}{\partial t} + \frac{\partial \bar{\mathbf{A}}_x \mathbf{q}}{\partial x} + \frac{\partial \bar{\mathbf{A}}_y \mathbf{q}}{\partial y} + \frac{\partial \bar{\mathbf{A}}_z \mathbf{q}}{\partial z} = \mathbf{f} \quad \text{on } \Omega_S \times \mathbf{I}, \qquad (8) $$

making use of the Jacobian matrices $\bar{\mathbf{A}}_x$, $\bar{\mathbf{A}}_y$ and $\bar{\mathbf{A}}_z$ associated with the flux functions along each of the coordinate direction. Their expression is classical and can be found for example in LeVeque (2002), chapter 22.

The expression of the flux $\mathcal{F} \mathbf{q} \cdot \mathbf{n}$ across an interface locally oriented by the unit normal $\mathbf{n} = (n_x, n_y, n_z)^T$ is simply obtained as a linear combination of matrices $\bar{\mathbf{A}}_x$, $\bar{\mathbf{A}}_y$ and $\bar{\mathbf{A}}_z$:

$$ \mathcal{F} \mathbf{q} \cdot \mathbf{n} = (n_x \bar{\mathbf{A}}_x + n_y \bar{\mathbf{A}}_y + n_z \bar{\mathbf{A}}_z) \mathbf{q} = \bar{\mathbf{A}}_n \mathbf{q}, \qquad (9) $$

where $\bar{\mathbf{A}}_n = n_x \bar{\mathbf{A}}_x + n_y \bar{\mathbf{A}}_y + n_z \bar{\mathbf{A}}_z$ is simply expressed in terms of the Lamé parameters as

$$ \bar{\mathbf{A}}_n = - \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & \frac{n_x}{\rho} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{n_y}{\rho} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{n_z}{\rho} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{n_z}{2\rho} & \frac{n_y}{2\rho} \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{n_z}{2\rho} & 0 & \frac{n_x}{2\rho} \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{n_y}{2\rho} & \frac{n_x}{2\rho} & 0 \\ (\lambda + 2\mu)n_x & \lambda n_x & \lambda n_x & 0 & 2\mu n_z & 2\mu n_y & 0 & 0 & 0 \\ \lambda n_y & (\lambda + 2\mu)n_y & \lambda n_y & 2\mu n_z & 0 & 2\mu n_x & 0 & 0 & 0 \\ \lambda n_z & \lambda n_z & (\lambda + 2\mu)n_z & 2\mu n_y & 2\mu n_x & 0 & 0 & 0 & 0 \end{pmatrix}. \qquad (10) $$

See Appendix B for a 2-D equivalent. Matrix $\bar{\mathbf{A}}_n$ is clearly composed of two non-zero blocks: an upper $6 \times 3$ block $\bar{\mathbf{A}}_{12}$ and a lower $3 \times 6$ block $\bar{\mathbf{A}}_{21}$ such that

$$ \bar{\mathbf{A}}_n = \begin{pmatrix} 0 & \bar{\mathbf{A}}_{12} \\ \bar{\mathbf{A}}_{21} & 0 \end{pmatrix}. \qquad (11) $$
## 2.2 Acoustic wave equations

Let us now consider the fluid domain $\Omega_F$. Acoustic wave is a small pressure disturbance that propagates through a compressible gas, causing infinitesimal changes in the density and pressure of inviscid fluid. Most sound waves are essentially a linear phenomenon: the magnitude of the disturbances from the background state is small enough that products or powers of the perturbation amplitude can be ignored. For the sake of simplicity we consider here the acoustic first-order adiabatic perturbation of an inviscid fluid, initially at rest, and in a state of hydrostatic equilibrium, neglecting the gravity effects.

The linearized acoustic wave equations can then be formulated as a hyperbolic system of first-order partial differential equations in a conservative form

$$ \frac{\partial \epsilon}{\partial t} = \nabla \cdot \mathbf{v} \quad \text{on } \Omega_F \times \mathbf{I}, \qquad (12a) $$

$$ \frac{\partial \rho \mathbf{v}}{\partial t} = \nabla(\kappa \epsilon) + \mathbf{f} \quad \text{on } \Omega_F \times \mathbf{I}, \qquad (12b) $$

where $\rho(\mathbf{x})$ is the reference mass density, $\mathbf{v}(\mathbf{x}, t)$ is the velocity field perturbation, and $\mathbf{f}(\mathbf{x}, t)$ is a generalized body force, together with the quantity $\epsilon(\mathbf{x}, t)$ related to the pressure via

$$ \epsilon = -\frac{1}{\kappa} p, \qquad (13) $$

where $p(\mathbf{x}, t)$ is the pressure perturbation and $\kappa(\mathbf{x})$ is the bulk modulus of the fluid.

The system of first-order partial differential equations (12) can be rewritten in compact form as,

$$ \frac{\partial \mathbf{q}}{\partial t} + \nabla \cdot \mathcal{F} \mathbf{q} = \mathbf{f} \quad \text{on } \Omega_F \times \mathbf{I} \qquad (14) $$

where the $\mathbf{q}$-vector associated with the state variables is now

$$ \mathbf{q} = \begin{pmatrix} \epsilon \\ \rho \mathbf{v} \end{pmatrix} = (\epsilon \quad \rho v_1 \quad \rho v_2 \quad \rho v_3)^T, \qquad (15) $$

and $\mathcal{F}\mathbf{q}$ is the associated flux matrix. The acoustic flux across an interface locally oriented by the normal $\mathbf{n} = (n_x, n_y, n_z)^T$ is given by $\mathcal{F}\mathbf{q} \cdot \mathbf{n} = \bar{A}_n \mathbf{q}$ where $\bar{A}_n$ is now given by

$$ \bar{A}_n = - \begin{pmatrix} 0 & \frac{n_x}{\rho} & \frac{n_y}{\rho} & \frac{n_z}{\rho} \\ \kappa n_x & 0 & 0 & 0 \\ \kappa n_y & 0 & 0 & 0 \\ \kappa n_z & 0 & 0 & 0 \end{pmatrix}. \qquad (16) $$

## 2.3 Initial and boundary conditions

The elastic-acoustic wave equations have to be completed with both initial and boundary conditions. The initial conditions are set consistently with the assumption that the medium is at rest in the reference configuration

$$ \mathbf{v}(\mathbf{x}, 0) = \mathbf{v}_0(\mathbf{x}) \quad \text{and} \quad \epsilon(\mathbf{x}, 0) = \epsilon_0(\mathbf{x}) \quad \text{on} \quad \Omega_S, \qquad (17a) $$

$$ \mathbf{v}(\mathbf{x}, 0) = \mathbf{v}_0(\mathbf{x}) \quad \text{and} \quad \epsilon(\mathbf{x}, 0) = \epsilon_0(\mathbf{x}) \quad \text{on} \quad \Omega_F. \qquad (17b) $$

Finally, Dirichlet and Neumann boundary conditions can be prescribed respectively as

$$ \mathbf{v} = \mathbf{v}^D \quad \text{on} \quad \bar{\Omega}_S \cap \Gamma^D \times \mathbf{I} \quad \text{and} \quad \mathbf{v} \cdot \mathbf{n} = \mathbf{v}_n^D \quad \text{on} \quad \bar{\Omega}_F \cap \Gamma^D \times \mathbf{I}, \qquad (18a) $$

$$ (\sigma \mathbf{n}) = \mathbf{t}^N \quad \text{on} \quad \bar{\Omega}_S \cap \Gamma^N \times \mathbf{I} \quad \text{and} \quad \mathbf{p} = \mathbf{t}^N \quad \text{on} \quad \bar{\Omega}_F \cap \Gamma^N \times \mathbf{I}. \qquad (18b) $$

## 3 ELASTIC-ACOUSTIC DOMAIN DECOMPOSITION AND TRANSMISSION RELATIONS

In this section we derive the transmission relations for heterogeneous elastic, acoustic and elasto-acoustic domain decompositions. The hybrid fields introduced in this section, will be then approximated as separate variables with the HDG method presented in Section 5.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** (a) Elastic–elastic domain decomposition: on the interface $\Gamma_{SS}$, the field $\hat{v}$ is single-valued while $\hat{\sigma}\hat{n}$ is double-valued. (b) Acoustic–acoustic domain decomposition: on $\Gamma_{FF}$ the field $\hat{v}_n$ is single-valued while $\hat{p}$ is double-valued. (c) Acoustic–elastic domain decomposition: on $\Gamma_{FS}$ the field $\hat{v}_n$ is single-valued while the hybrid pressure $\hat{p}$ and the hybrid traction $\hat{\sigma}\hat{n}$ are defined on the acoustic and elastic sides, respectively.</figcaption></figure>

## 3.1 Elastic–elastic domain decomposition

Let us consider here the decomposition of an elastic domain $\Omega_S$ into two heterogeneous subdomains separated by an interface $\Gamma_{SS}$

$$ \bar{\Omega} = \bar{\omega}^{+} \cup \bar{\omega}^{-} \quad \text{with} \quad \Gamma_{SS} = \bar{\omega}^{+} \cap \bar{\omega}^{-}. \qquad (19) $$

On the solid–solid interface $\Gamma_{SS}$, both velocities and tractions are required to be continuous—classical kinematic and dynamic boundary conditions in continuum mechanics

$$ \|\mathbf{v}(\mathbf{x}, t)\| = 0 \quad \forall \mathbf{x} \in \Gamma_{SS}, \forall t \in \mathbf{I} \qquad (20a) $$

$$ \|\sigma(\mathbf{x}, t)n(\mathbf{x})\| = 0 \quad \forall \mathbf{x} \in \Gamma, \forall t \in \mathbf{I} \qquad (20b) $$

where $[\cdot]$ stands for the jump operator across the oriented interface, that is, for any tensor $\mathbf{a}$

$$ \|\mathbf{a}\| = \mathbf{a}^{-} - \mathbf{a}^{+} \quad \text{with} \quad \left\{ \begin{array}{l} \mathbf{a}^{-} = \lim_{\varepsilon \to 0} \mathbf{a}(\mathbf{x} - \varepsilon\mathbf{n}) \\ \mathbf{a}^{+} = \lim_{\varepsilon \to 0} \mathbf{a}(\mathbf{x} + \varepsilon\mathbf{n}) \end{array} \right. \quad \text{for all } \mathbf{x} \in \Gamma_{SS}. \qquad (21) $$

### 3.1.1 Hybridization of the elastic dynamic conditions

Let us introduce two kinds of fields on the solid–solid interface: a single valued *hybrid* velocity field $\hat{v}$ defined on $\Gamma_{SS}$, and the traction $\hat{\sigma}\hat{n}$ defined on the boundary of $(\bar{\omega}^{-} \cap \Gamma_{SS}) \cup (\bar{\omega}^{+} \cap \Gamma_{SS})$, which is double-valued with $\hat{\sigma}\hat{n}^{-}$ and $\hat{\sigma}\hat{n}^{+}$ on the each side of the interface (see Fig. 1a). Continuity of the velocities across the solid-solid interface is then ensured by the single-valued definition of $\hat{v}$ while the traction continuity can now be weakly prescribed as

$$ \int_{\Gamma_{SS}} \|\hat{\sigma}\hat{n}\| \cdot \mu \, d\Gamma = 0 \quad \text{for all } \mu \in H^{1/2}(\Gamma_{SS}). \qquad (22) $$

This last relation is called here the *conservativity condition*, to be consistent with the HDG terminology. As a side note, the hybrid field is also called the *trace* in the literature since it may be seen as the result of a mathematical *trace operator*.

It is worth to note here that, for each subdomain, either the tractions $\hat{\sigma}\hat{n}$ or the velocity $\hat{v}$ can be implicitly defined respectively as a function of each other. Following our approach, we choose to define the tractions $\hat{\sigma}\hat{n}$ as the responses of the elastic bodies to the prescribed velocity $\hat{v}$, that is, we introduce the following DtN operators

$$ \hat{\sigma}\hat{n}^{-} = \mathrm{DtN}^{-}(\hat{v}) \quad \text{and} \quad \hat{\sigma}\hat{n}^{+} = \mathrm{DtN}^{+}(\hat{v}). \qquad (23) $$

When the DtN operators are invertible, enforcing the conservativity condition (22) is, ultimately, solving for $\hat{v}$. Once $\hat{v}$ is known, $\hat{\sigma}\hat{n}^{\pm}$ can be determined through (23), making the local resolution of elastic wave problem independent on each subdomain $\omega^{\pm}$.
### 3.1.2 Hybrid variables and the relaxation of dynamic conditions

We briefly recall, in a loose manner, some important notions related to finite element methods. The goal of this paragraph is to underline how and why hybrid variables ($\hat{v}$, $\hat{\sigma n}$) and hybridization of methods, has arisen in the finite element context; and why it leads to a relaxation of the dynamic condition (20b).

In classical CG finite element methods, one seeks approximate solutions of (2) by solving a weak form of the sole (2b). As underlined in the general introduction of the paper, $C_0$-continuity of the primal variable (velocity $v$), is then ensured by construction (and therefore, (20a) is respected); nothing is said about its dual (stress $\sigma$ or strain $\epsilon$), and (20b) is only approximated.

An alternative approach is to seek approximate solutions for both primal and dual variables, by solving weak forms of (2b) and (2a); this leads to the so-called mixed finite element methods (e.g. Brezzi & Fortin 1991). Unfortunately, the discretization of these mixed methods leads to a non-definite, large linear algebraic system, the solution of which is computationally expensive. One reason for this is the continuity requirement for the tractions, precisely (20b). One may choose to relax this constraint by working with non-conforming elements (Cohen 2002), or by introducing an additional field on the element boundaries—in such a way that static elimination of primal and dual variables can be performed at the element level, yielding a much smaller matrix equation to solve for the additional field. This is precisely the hybrid variable—equivalent of our field $\hat{v}$—introduced for the first time by Fraejis de Veubeque (1965), and regarded for a long time as an implementation trick only. Later, Brezzi & Fortin (1991) noted that this hybrid field can be seen as the Lagrange multiplier associated with the constraint of weak continuity for the tractions, namely (22). Moreover, Arnold & Brezzi (1985) showed that hybridization was more than a trick, and that the hybrid variable, being the trace of the primal variable, carries extra-information on the solution, that can be used to post-process the volumetric solution, in order to obtain super-convergence properties.

Hybridization was traditionally (Arnold & Brezzi 1985) carried out with rather involved algebraic manipulations where both primal, dual and hybrid variables were implied; Cockburn & Gopalakrishnan (2004) provided a clever procedure to lighten that task, and make it local and explicit. Furthermore, Cockburn et al. (2009) bridged the gap between different variational frameworks, and showed that both mixed and some of DG methods could be hybridized. For the time being, we just retain that building an hybridized method, leads to relaxed dynamic boundary conditions (22), and that the addition of an extra, hybrid variable can offer superconvergence properties.

## 3.2 Acoustic-acoustic domain decomposition

Let us now consider an acoustic domain $\Omega_F$ with the same domain decomposition (19), with now a fluid-fluid interface $\Gamma_{FF}$.

Across the fluid-fluid interface, only the pressure and the normal component of the velocity vector have to be continuous, leading to the conditions:

$$ \|v(x, t) \cdot n(x)\| = 0 \quad \forall x \in \Gamma_{FF}, \forall t \in I \qquad (24a) $$

$$ \|p(x, t)\| = 0 \quad \forall x \in \Gamma_{FF}, \forall t \in I. \qquad (24b) $$

A similar strategy leads to introduce a single-valued hybrid scalar field $\hat{v}_n$, defined on $\Gamma_{FF}$, as the normal component of the velocity; and the hybrid pressure $\hat{p}$ defined on $(\bar{\omega}^- \cap \Gamma_{FF}) \cup (\bar{\omega}^+ \cap \Gamma_{FF})$, which is double-valued on the respective side of $\Gamma_{FF}$ (see Fig. 1b).

Continuity of the normal component of the velocity is again insured by the single-valued definition of $\hat{v}_n$ while the pressure continuity can now be weakly enforced as

$$ \int_{\Gamma_{FF}} \|\hat{p}\| \mu \, d\Gamma = 0 \quad \text{for all } \mu \in H^{1/2}(\Gamma_{FF}), \qquad (25) $$

defining the conservativity condition for the fluid-fluid interface.

The pressures $\hat{p}^{\pm}$ can be implicitly expressed as a function of $\hat{v}_n$ through some acoustic DtN operators

$$ \hat{p}^{-} = \mathrm{DtN}^{-}(\hat{v}_{n}) \quad \text{and} \quad \hat{p}^{+} = \mathrm{DtN}^{+}(\hat{v}_{n}). \qquad (26) $$

When the DtN operators are invertible, enforcing the conservativity condition (25) is, ultimately, solving for $\hat{v}_n$. And once $\hat{v}_n$ is known, the $\hat{p}^{\pm}$ can be retrieved using (26), making the local resolution independent on each subdomain $\omega^{\pm}$.

## 3.3 Acoustic-elastic domain decomposition

Let us consider finally an elasto-acoustic domain $\Omega$ decomposed into an acoustic subdomain $\omega^-$ and an elastic subdomain $\omega^+$ separated by an interface $\Gamma_{FS}$.

On fluid-solid interfaces, an inviscid fluid can slip along the solid surface, but cannot be affected by solid shear stresses. Therefore, the normal component of the velocity and the pressure should be continuous across fluid-solid interfaces:

$$ \|v(x, t) \cdot n(x)\| = 0 \quad \forall x \in \Gamma_{FS}, \forall t \in I \qquad (27a) $$

$$ \sigma(x, t) n(x) = p(x, t) n(x) \quad \forall x \in \Gamma_{FS}, \forall t \in I. \qquad (27b) $$
<figure><img src="image_5.png" /><figcaption>Figure 2. Riemann problem in the vicinity of an interface point $x_{\Gamma} \in \Gamma$.</figcaption></figure>

We now introduce four fields, see Fig. 1, defined as: (1) the elastic hybrid velocity $\hat{v}$, defined on $\Gamma_{FS}$; (2) the hybrid normal component of the acoustic velocity $\hat{v}_n = \hat{v} \cdot n$ on $\Gamma_{FS}$; (3) the elastic traction $\hat{\sigma}\hat{n}^+$, defined on $\bar{\omega}^+ \cap \Gamma_{FS}$; (4) the acoustic pressure $\hat{p}^-$, defined on $\bar{\omega}^- \cap \Gamma_{FS}$.

Thus, by construction, the first fluid-solid boundary condition (27a) is satisfied while the second one (27b) can be enforced weakly as

$$ \int_{\Gamma_{FS}} \hat{p} \, n \cdot \mu \, d\Gamma = \int_{\Gamma_{FS}} \hat{\sigma}\hat{n} \cdot \mu \, d\Gamma \quad \text{for all } \mu \in H^{1/2}(\Gamma_{FS}). \qquad (28) $$

This conservativity condition (28) is the equivalent of (22) and (25) for acoustic-elastic interfaces.

The pressure $\hat{p}^-$ can be expressed as an implicit function of $\hat{v}_n$ using the acoustic DtN operator (26), and the tractions $\hat{\sigma}\hat{n}^+$ as an implicit function of $\hat{v}$ using the elastic DtN operator (23). When these DtN operators are invertible, enforcing (28) is, at the end, solving for $\hat{v}$ and then for $\hat{v}_n$. Once $\hat{v}$ and $\hat{v}_n$ are known, $\hat{p}^-$ is retrieved using (26), and $\hat{\sigma}\hat{n}^+$ using (23), making the local resolution independent on each subdomain $\omega^{\pm}$.

# 4 EXPLICIT CONSTRUCTION OF THE DTN OPERATORS USING RIEMANN SOLVER

The DtN operators introduced above can be built explicitly and locally through the exact resolution of the associated Riemann problem at each point of the interface, when identifying the hybrid variables $\hat{v}$, $\hat{v}_n$, $\hat{\sigma}\hat{n}$ and $\hat{p}$ with the state solution of the Riemann problem on the interface. Practically, this leads to a constructive method for the DtN operators that makes use of the Rankine-Hugoniot conditions on either side of the interface. As such the double-valued traction on $\Gamma$ are computed using information from only one side of the interface and from the single-valued hybrid velocity. The hybrid variable is therefore the only unknown coupling the subdomains, and can be implicitly determined using the conservativity conditions.

This leads to *hybridized* versions of the Riemann solvers previously presented by Wilcox et al. (2010) for heterogeneous acoustic-elastic interfaces. It is *hybridized* in the sense that the *hybrid* variable, that is, the hybrid velocity variable is explicitly computed first, and then used to compute independently the other variables, contrary to the direct computation of the numerical fluxes performed by classical Riemann solvers (LeVeque 2002; Wilcox et al. 2010; Käser & Dumbser 2006).

A more general study of hybridized upwind fluxes for a wide class of PDEs has been proposed by Bui-Thanh (2015). In some ways, this section is an extension of this last work for the elastic-acoustic wave equations.

## 4.1 Elastic-elastic interfaces

Let us consider an elastic-elastic interface as in Section 3.1. For each point of the interface $x \in \Gamma_{SS}$, we consider the evolution of quantities $q = (\epsilon, \rho v)^T$ along the normal direction of the interface (represented in Fig. 2).

On the interface, three *a priori* different quantities appear: the continuous prolongation $q|_{\omega^-}$ and $q|_{\omega^+}$ on $\Gamma$ denoted $q^-$ and $q^+$ respectively, and the fields $\hat{v}$ and $\hat{\sigma}\hat{n}^\pm$ defined only on the very interface.

Without any restriction, let us assume piecewise constant state variables ($q^-$, $q^+$) and Jacobian flux matrices ($\bar{A}_n^-$, $\bar{A}_n^+$) in the neighbourhood of the interface, as represented in Fig. 3. Then, finding the time evolution of this configuration leads to solve the associated Riemann problem for all points $x \in \Gamma_{SS}$.

Since the problem is linear, classical solution of the Riemann problem involves the diagonalization of the oriented Jacobian flux matrices $\bar{A}_n$. For isotropic linear elasticity, the eigenvalues are the elastic wave speeds

$$ \lambda_1 = -c_p, \quad \lambda_2 = -c_s, \quad \lambda_3 = 0, \quad \lambda_4 = c_s, \quad \lambda_5 = c_p, $$

where the compressional and shear wave speeds are explicit functions of Lamé parameters

$$ c_p = \sqrt{\frac{\lambda + 2\mu}{\rho}} \quad \text{and} \quad c_s = \sqrt{\frac{\lambda}{\rho}}. \qquad (29) $$
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Left: Riemann problem associated with discontinuous states of $q$ at the interface. Right: typical solution for a Riemann problem associated with isotropic elastodynamic equations: four waves are propagating away from the interface.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Negative characteristic lines used to build $(\sigma n)_2$.</figcaption></figure>

Each of these eigenvalues has a multiplicity of 1 in 2-D, while in 3-D, $\lambda_2$ and $\lambda_4$ have a multiplicity of 2, and $\lambda_3$ a multiplicity of 3. The diagonalization of $\bar{\bar{A}}_n$ in 3-D is therefore

$$ \bar{\bar{A}}_n = \bar{\bar{R}}\bar{\bar{\Lambda}}\bar{\bar{R}}^{-1} \quad \text{with} \quad \bar{\bar{\Lambda}} = \operatorname{diag}(-c_p, -c_s, -c_s, 0, 0, 0, c_s, c_s, c_p). \qquad (30) $$

These eigenvalues are associated with characteristic lines represented on Fig. 4. Matrices $\bar{\bar{A}}_n^-$ and $\bar{\bar{A}}_n^+$ may have different eigenvalues $c_p^-, c_s^-$ and $c_p^+, c_s^+$ since parameters $\lambda, \mu$ and $\rho$ may be discontinuous at the interface.

The solution of the Riemann problem, for $t \ge 0$ leads to six different states $\{q^-, q^1, q^2, q^3, q^4, q^+\}$. The Rankine–Hugoniot conditions state that each jump between two adjacent states is an eigenvector of either $\bar{\bar{A}}_n^-$ (if the jump propagates in $\omega^-$) or $\bar{\bar{A}}_n^+$ (if the jump propagates in $\omega^+$). This allows to determine explicitly all the intermediate states $q^1, \dots, q^4$. For an elastic–elastic interface, the Rankine–Hugoniot conditions are

$$ c_p^- (\mathbf{q}^- - \mathbf{q}_1) + \bar{\bar{A}}_n^- (\mathbf{q}^- - \mathbf{q}_1) = 0, \qquad (31a) $$

$$ c_s^- (\mathbf{q}_1 - \mathbf{q}_2) + \bar{\bar{A}}_n^- (\mathbf{q}_1 - \mathbf{q}_2) = 0, \qquad (31b) $$

$$ \bar{\bar{A}}_n^- \mathbf{q}_2 - \bar{\bar{A}}_n^+ \mathbf{q}_3 = 0, \qquad (31c) $$

$$ -c_s^+(\mathbf{q}_3 - \mathbf{q}_4) + \bar{\bar{A}}_n^+(\mathbf{q}_3 - \mathbf{q}_4) = 0, \qquad (31d) $$

$$ -c_p^+(\mathbf{q}_4 - \mathbf{q}^+) + \bar{\bar{A}}_n^+(\mathbf{q}_4 - \mathbf{q}^+) = 0. \qquad (31e) $$

The Rankine–Hugoniot condition (31c) associated with the stationary wave $\lambda_3 = 0$ insures that the Riemann solution leads to continuous velocity and traction states at the interface.

$$ \mathbf{v}_2 = \mathbf{v}_3, \qquad (32a) $$

$$ (\sigma n)_2 = (\sigma n)_3. \qquad (32b) $$

Let us identify $\hat{v}$ with the velocity state at the interface: $\hat{v} = \mathbf{v}_2 = \mathbf{v}_3$. The single-valued definition of $\hat{v}$ will insure by construction the first condition (32a). The pointwise second condition (32b) is now relaxed and will be satisfied only weakly using the conservativity condition (22) and identifying

$$ \widehat{\sigma \mathbf{n}^-} = (\sigma \mathbf{n})_2 \quad \text{and} \quad \widehat{\sigma \mathbf{n}^+} = (\sigma \mathbf{n})_3. \qquad (33) $$

Let us now express $\widehat{\sigma \mathbf{n}^-} = (\sigma \mathbf{n})_2$ as a function of the state $\mathbf{q}^-$ and of the single-valued hybrid velocity $\hat{v}$, using only the characteristics directions travelling towards the negative side of the interface (the ones on the left on Fig. 4).
It is worth to note here that the matrices $\tau^-$ and $\tau^+$ are symmetric positive definite.

## 4.2 Acoustic-acoustic interfaces

It is not possible to simply use the expressions (38) and (39) by setting $c_s^- = 0$ and $c_s^+ = 0$ in order to retrieve the DtN operators for an acoustic-acoustic interface, because these DtN would not be invertible. Indeed only the normal component of the velocity is continuous across such an interface. Therefore, we need to follow a development similar to the one presented for the elastic case, but now using the acoustic-acoustic Riemann problem, with the acoustic state variables $\mathbf{q}$ and acoustic Jacobian flux matrices $\bar{\mathbf{A}}_n$ provided by (16).

For an acoustic-acoustic interface, there are only three different eigenvalues after diagonalization of the matrices $\bar{\mathbf{A}}_n$. These eigenvalues are associated with the pressure wave velocities

$$ \lambda_1 = -c_p = -\sqrt{\frac{\kappa}{\rho}}; \quad \lambda_2 = 0; \quad \lambda_3 = c_p. $$

The second eigenvalue $\lambda_2 = 0$ has a multiplicity of 2 in 3-D and of 1 in 2-D. Therefore, there will be only two intermediary states ($\mathbf{q}_1, \mathbf{q}_2$) and three Rankine-Hugoniot conditions

$$ c_p^- (\mathbf{q}^- - \mathbf{q}_1) + \bar{\mathbf{A}}_n^- (\mathbf{q}^- - \mathbf{q}_1) = 0, \quad (40a) $$

$$ \bar{\mathbf{A}}_n^- \mathbf{q}_1 - \bar{\mathbf{A}}_n^+ \mathbf{q}_2 = 0, \quad (40b) $$

$$ -c_p^+(\mathbf{q}_2 - \mathbf{q}^+) + \bar{\mathbf{A}}_n^+(\mathbf{q}_2 - \mathbf{q}^+) = 0. \quad (40c) $$

Interface conditions (40b) are different from (32), since they express only the continuity of the pressure and of the normal component of the velocity,

$$ \mathbf{v}_1 \cdot \mathbf{n} = \mathbf{v}_2 \cdot \mathbf{n}, \quad (41a) $$

$$ \kappa^- \epsilon_1 = \kappa^+ \epsilon_2. \quad (41b) $$

Let us now identify the single-valued $\hat{v}_n$ with the continuous normal component of the velocity at the interface: $\hat{v}_n = \mathbf{v}_1 \cdot \mathbf{n} = \mathbf{v}_2 \cdot \mathbf{n}$. By construction the first acoustic continuity condition of (24) is then satisfied. The second condition will again be weakly imposed using the acoustic conservativity condition (28), identifying

$$ \hat{p}^- = -\kappa^- \epsilon_1 \quad \text{and} \quad \hat{p}^+ = -\kappa^+ \epsilon_2. \quad (42) $$

The construction of the acoustic DtN$^-$ operator, that is, the relation linking $\hat{p}^-$ to $\mathbf{q}^-$ and to $\hat{v}_n$, is straightforward using the second part of vectorial equality (40a)

$$ c_p^- \rho^- (\mathbf{v}^- - \mathbf{v}_1) - (\kappa^- \epsilon^- \mathbf{n} - \kappa^- \epsilon_1 \mathbf{n}) = 0. \quad (43) $$

Taking the scalar product of this last relation by $\mathbf{n}$, and using (42), the acoustic DtN$^-$ operator expression follows

$$ -\hat{p}^- = \mathrm{DtN}^-(\hat{v}_n) = (\kappa \epsilon)^- - \tau^- (\mathbf{v}^- \cdot \mathbf{n} - \hat{v}_n). \quad (44) $$

Here $\tau^-$ is the acoustic impedance, and is the scalar counterpart of the elastic $\tau^-$ matrix

$$ \tau^- = \rho^- c_p^- \quad (45) $$

Similarly, using the relation (40c), the DtN$^+$ operator is given as

$$ -\hat{p}^+ = \mathrm{DtN}^+(\hat{v}_n) = (\kappa \epsilon)^+ + \tau^+ (\mathbf{v}^+ \cdot \mathbf{n} - \hat{v}_n), \quad (46) $$

with $\tau^+ = \rho^+ c_p^+$ being the acoustic impedance on the positive side of the interface. It is worth to note here that these values of the acoustic $\tau^+$ are coherent with the ones proposed independently by Stanglmeier et al. (2016).

## 4.3 Acoustic-elastic interfaces

In the case of an acoustic-elastic interface, the Riemann problem deals with different set of equations on both sides of the interface. The diagonalization of the oriented Jacobian matrices $\bar{\mathbf{A}}_n$ leads to different number of eigenvalues on the acoustic and the elastic sides. Physically, pressure waves on the acoustic media are coupled to both compressional and shear waves of the elastic part. It means that there is one characteristic line in the acoustic side, and two characteristic lines in the elastic side (see Fig. 7). They correspond to the following four different eigenvalues:

$$ \lambda_1 = -c_p^-; \quad \lambda_2 = 0; \quad \lambda_3 = c_s^+; \quad \lambda_4 = c_p^+. \quad (47) $$
<figure><img src="image_8.png" /><figcaption>Figure 7. Characteristic lines for an acoustic-elastic interface.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Lamb's Test configuration and velocity norm fields at times $t = 0.75$ s (top) and $t = 2$ s (bottom). The vertical force near the free surface generates a P-wave (a) and an S-wave (b). Both waves are linked by the head wave (c). Surface Rayleigh wave (d) is travelling along the free surface with shallow depth penetration.</figcaption></figure>

Therefore, there will be three intermediate states ($\mathbf{q}_1, \mathbf{q}_2, \mathbf{q}_3$) and three Rankine-Hugoniot conditions, one for the acoustic side and two for the elastic side:

$$c_p^- (\mathbf{q}^- - \mathbf{q}_1) + \bar{\bar{A}}_n^- (\mathbf{q}^- - \mathbf{q}_1) = 0, \qquad (48a)$$

$$-c_s^+(\mathbf{q}_2 - \mathbf{q}_3) + \bar{\bar{A}}_n^+(\mathbf{q}_2 - \mathbf{q}_3) = 0, \qquad (48b)$$

$$-c_p^+(\mathbf{q}_3 - \mathbf{q}^+) + \bar{\bar{A}}_n^+(\mathbf{q}_3 - \mathbf{q}^+) = 0. \qquad (48c)$$

The states $\mathbf{q}$ and the Jacobian matrices $\bar{\bar{A}}_n$ in relations (48) shall be understood as (7) and (10), or as (15) and (16) according to the media in which the Rankine condition is applied. The system (48) has to be completed by the continuity condition (27) on the acoustic-elastic interface:

$$\mathbf{v}_1 \cdot \mathbf{n} = \mathbf{v}_2 \cdot \mathbf{n}, \qquad (49a)$$

$$\kappa^{-} \epsilon_{1} = - \hat{p} \mathbf{n} = (\sigma \mathbf{n})_{2}. \qquad (49b)$$
Although only the normal component $\hat{v}_n$ is continuous, we can still define a single-valued interface velocity $\hat{v}$ such that $\hat{v} = v_2$, and $\hat{v}_n = \hat{v} \cdot n$, which means that its tangential component is equal to the one of the elastic part. The last continuity condition will again be relaxed and only weakly imposed using the acoustic-elastic conservativity condition (28).

On the acoustic side, working with the condition (48a), leads to the relation (44). Now multiplying this last relation by the normal $n$ to retrieve $\hat{p} n$, one gets

$$- \hat{p} n = (\kappa^{-} \epsilon_{1}) n = (\kappa \epsilon)^{-} n - \tau (v \cdot n - \hat{v} \cdot n) n. \qquad (50)$$

This relation can be rewritten under the form

$$- \hat{p} n = (\kappa^{-} \epsilon_{1}) n = (\kappa \epsilon)^{-} n - \tau_{ac}^{-} (v - \hat{v}), \qquad (51)$$

using now the matrix $\tau_{ac}$, similar to the elastic one with $c_s^- = 0$:

$$\tau_{ac}^{-} = \begin{pmatrix} \rho^{-} c_{p}^{-} n_{x}^{2} & \rho^{-} c_{p}^{-} n_{x} n_{y} & \rho^{-} c_{p}^{-} n_{x} n_{z} \\ \rho^{-} c_{p}^{-} n_{x} n_{y} & \rho^{-} c_{p}^{-} n_{y}^{2} & \rho^{-} c_{s}^{-} n_{y} n_{z} \\ \rho^{-} c_{p}^{-} n_{x} n_{z} & \rho^{-} c_{p}^{-} n_{y} n_{z} & \rho^{-} c_{p}^{-} n_{z}^{2} \end{pmatrix}. \qquad (52)$$

On the elastic side, similar development used for the elastic-elastic interface allows to recover the DtN operator of equation (39)

$$\hat{\sigma} \hat{n} = (\sigma n)^{+} - \tau^{+} (v^{+} - \hat{v}), \qquad (53)$$

with the matrix $\tau^+$ is given by (37), taking material parameters at the positive side of the interface $(\rho^+, c_p^+, c_s^+)$.

## 4.4 Computing the velocities $\hat{v}$ and $\hat{v}_n$

When the conservativity conditions are satisfied pointwise, the single-valued hybrid unknowns $\hat{v}$ and $\hat{v}_n$ can also be computed pointwise using the previous results:

(i) for an elastic-elastic interface, applying condition (32b) and using (38) with (39) leads to

$$\hat{v} = (\tau^{+} + \tau^{-})^{-1} (\tau^{-} v^{-} + \tau^{+} v^{+} - [\sigma n]) \qquad (54)$$

(ii) for an acoustic-acoustic interface, applying condition (41b) and using (44) on both sides leads to

$$\hat{v}_n = \frac{1}{\tau^+ + \tau^-} (\tau^- v^- \cdot n - \tau^+ v^- \cdot n - [\kappa \epsilon]) \qquad (55)$$

(iii) for an acoustic-elastic interface, applying condition (49b) and using (51) with (53) leads to

$$\hat{v} = (\tau^{+} + \tau_{ac}^{-})^{-1} (\tau_{ac}^{-} v^{-} + \tau^{+} v^{+} - p^{-} n + (\sigma n)^{+}) \qquad (56)$$

Once the hybrid velocity are determined, the traction and pressure are computed using the DtN operators established in the previous sections, which completes the description of the Riemann solver.

## 4.5 Hybridized and classical Riemann solvers

It can be ascertained that the results given by this two-stage hybridized Riemann solver are equivalent, at least for the elastic-elastic case, to the numerical fluxes proposed by Wilcox et al. (2010) using a classical Riemann solver. Indeed the same initial states are considered, as well as the same Rankine-Hugoniot relations, although they are used in a different way. Both Riemann solvers end up producing the same results, under the form of numerical fluxes $\mathfrak{F}q^*$ (DG) or hybrid unknowns $\hat{v}$ and $\hat{\sigma}n$ (HDG). More precisely, the relation between them is given by

$$\mathfrak{F}q^* = \begin{pmatrix} \rho \bar{A}_{12} \hat{v} \\ - \hat{\sigma}n \end{pmatrix} \qquad (57)$$

using Wilcox notation for numerical fluxes.

Moreover, if Wilcox et al. (2010) had followed a scalar strain-velocity formulation, or a pressure-velocity formulation for the acoustic equations, their acoustic fluxes would also have been related to our hybrid unknowns $\hat{v}_n$ and $\hat{p}$ through a simple relation as well.

An important point is that our two-stage Riemann solver does not induce an extra computational complexity. When matrices $\tau^\pm$, $\tau_{ac}^-$, $(\tau^+ + \tau^-)^{-1}$ and $(\tau^+ + \tau_{ac}^-)^{-1}$ are pre-computed, our solver needs two matrix-vector products on each side of the interface, and another one at the interface to get all the hybrid unknowns, while classical solvers involve rather complex numerical fluxes expressions (Wilcox et al. 2010). Numerical results from Section 7.3 suggest that the two-stage Riemann solver is at least as fast as the classical one.
# 5 HDG DISCRETE FORMULATION

## 5.1 The HDG scheme

Solid and fluid domains are decomposed into a set of non-overlapping open subdomains. Let $\mathcal{T}_h^S$ be the partition associated with the solid domain $\Omega_S$ with $n_s$ non-overlapping subdomains and $\mathcal{T}_h^F$ the partition associated with the fluid region $\Omega_F$ with $n_f$ non-overlapping subdomains:

$$ \mathcal{T}_h^S = \{K_i\}_{i=1}^{n_s}, \quad \bar{\Omega}_S = \bigcup_{i=1}^{n_s} \bar{K}_i, \qquad (58a) $$

$$ \mathcal{T}_h^F = \{K_i\}_{i=n_s+1}^{n_s+n_f}, \quad \bar{\Omega}_F = \bigcup_{i=n_s+1}^{n_s+n_f} \bar{K}_i, \qquad (58b) $$

and the global partition is denoted $\mathcal{T}_h = \mathcal{T}_h^S \cup \mathcal{T}_h^F$.

The set of fluid and solid boundaries are defined respectively as $\partial \mathcal{T}_h^S = \{\partial K : K \in \mathcal{T}_h^S\}$ and $\partial \mathcal{T}_h^F = \{\partial K : K \in \mathcal{T}_h^F\}$. The solid-fluid set of faces (skeleton) is now defined as $\mathcal{E}_h^S = \{\partial K_k \cap \partial K_l : K_k \in \mathcal{T}_h^S, K_l \in \mathcal{T}_h^F\}$, while similarly $\mathcal{E}_h^S$ and $\mathcal{E}_h^F$ denotes respectively the solid and the fluid mesh skeleton. Finally the union of all interfaces, including the outer surface, is labelled $\mathcal{E}_h$ and is such that $\mathcal{E}_h : \mathcal{E}_h^S \cup \mathcal{E}_h^F \cup \mathcal{E}_h^{SF}$. It is worth to note here that in the set $\partial \mathcal{T}_h$ internal interfaces are counted twice and define the support of the double-valued interface fields, while $\mathcal{E}_h$ denotes the skeleton that defines the support of single-valued interface variables.

Let us consider the following broken polynomial spaces for both the approximate solution and test functions

$$ \mathcal{Q}_h = \{r \in L^2(\mathcal{T}_h^F) : r|_{K} \in \mathcal{P}_N(K), \forall K \in \mathcal{T}_h^F\}, \qquad (59a) $$

$$ \mathcal{V}_h = \{w \in (L^2(\mathcal{T}_h^S))^{n_d} : w|_K \in (\mathcal{P}_N(K))^{n_d}, \forall K \in \mathcal{T}_h^S\}, \qquad (59b) $$

$$ \mathcal{S}_h = \left\{ \mathbf{r} \in (L^2(\mathcal{T}_h^S))^{n_g} \text{ sym} : \mathbf{r}|_K \in (\mathcal{P}_N(K))^{n_g} \text{ sym}, \forall K \in \mathcal{T}_h^S \right\}, \qquad (59c) $$

$$ \hat{\mathcal{V}}_h(\mathbf{g}) = \{\mu \in L^2(\mathcal{E}_h^S) : \mu|_F \in (\mathcal{P}_N(F))^{n_d}, \forall F \in \mathcal{E}_h^S \text{ and } \mu = \mathbf{g} \text{ on } \bar{\Omega}_S \cap \Gamma^D\}, \qquad (59d) $$

$$ \hat{\mathcal{V}}_h(\mathbf{g}) = \{\mu \in L^2(\mathcal{E}_h^F) : \mu|_F \in \mathcal{P}_N(F), \forall F \in \mathcal{E}_h^F \text{ and } \mu = \mathbf{g} \text{ on } \bar{\Omega}_F \cap \Gamma^D\}. \qquad (59e) $$

With a slight abuse of notation, statements like $\mathbf{v}_h \in \mathcal{V}_h$ for time-dependent functions shall be understood as $\mathbf{v}_h(., t) \in \mathcal{V}_h, \forall t \in \mathbf{I}$.

As in Section 3, we define the following fields on the mesh faces:

(i) in the elastic part, a single-valued vector field $\hat{\mathbf{v}} \in \hat{\mathcal{V}}_h(\mathbf{v}^D)$ with support $\mathcal{E}_h^S$, and a double-valued vector field $\hat{\boldsymbol{\sigma}}\mathbf{n}$ with support $\partial \mathcal{T}_h^S$,

(ii) in the acoustic part, a single-valued scalar field $\hat{v}_h \in \hat{\mathcal{V}}_h(v_h^D)$ with support $\mathcal{E}_h^F$, and a double-valued scalar field $\hat{p}$ with support $\partial \mathcal{T}_h^F$.

Moreover, results from Section 4 allow to express $\hat{\boldsymbol{\sigma}}\mathbf{n}$ and $\hat{p}$ as

$$ \hat{\boldsymbol{\sigma}}\mathbf{n} = (\boldsymbol{\sigma}\mathbf{n}) - \boldsymbol{\tau}(\mathbf{v} - \hat{\mathbf{v}}) \quad \text{and} \quad \hat{p} = p + \boldsymbol{\tau}(\mathbf{v} \cdot \mathbf{n} - \hat{v}_h). \qquad (60) $$

The proposed HDG method is built considering the sum of the local weak forms of the elastic and acoustic waves equations on each element, enforcing the conservativity conditions (22), (25) and (28) on each internal face of the mesh using the DtN expressions (60).

We then search functions $(\epsilon_h, \mathbf{v}_h, \hat{v}_h) \in \mathcal{S}_h \times \mathcal{V}_h \times \hat{\mathcal{V}}_h(v^D)$ on the elastic domain, and functions $(\epsilon_h, \mathbf{v}_h, \hat{v}_h) \in \mathcal{Q}_h \times \mathcal{V}_h \times \hat{\mathcal{V}}_h(v^D)$ on the acoustic domain such that:

In the elastic domain, for all $(\mathbf{r}, \mathbf{w}, \boldsymbol{\mu}) \in \mathcal{S}_h \times \mathcal{V}_h \times \hat{\mathcal{V}}_h(0)$

$$ \int_{\mathcal{T}_h^F} \frac{\partial \epsilon_h}{\partial t} : \mathbf{r} \, d\Omega = - \int_{\mathcal{T}_h^S} \mathbf{v}_h \cdot (\nabla \cdot \mathbf{r}) \, d\Omega + \int_{\partial \mathcal{T}_h^S} \hat{\mathbf{v}}_h \cdot (\mathbf{r} \cdot \mathbf{n}) \, d\Gamma, \qquad (61a) $$

$$ \int_{\mathcal{T}_h^F} \rho \frac{\partial \mathbf{v}_h}{\partial t} \cdot \mathbf{w} \, d\Omega = \int_{\mathcal{T}_h^S} \nabla \cdot (\mathbf{C} : \boldsymbol{\epsilon}_h) \cdot \mathbf{w} \, d\Omega - \int_{\partial \mathcal{T}_h^S} \boldsymbol{\tau}(\mathbf{v}_h - \hat{\mathbf{v}}_h) \cdot \mathbf{w} \, d\Gamma + \int_{\mathcal{T}_h^S} \mathbf{f} \cdot \mathbf{w} \, d\Omega, \qquad (61b) $$

$$ \int_{\partial \mathcal{T}_h^S} [\hat{\boldsymbol{\sigma}}\mathbf{n}_h] \cdot \boldsymbol{\mu} \, d\Gamma = \int_{\Omega_S \cap \Gamma^N} \mathbf{t}^N \cdot \boldsymbol{\mu} \, d\Gamma. \qquad (61c) $$

In the acoustic domain, for all $(\mathbf{r}, \mathbf{w}, \boldsymbol{\mu}) \in \mathcal{Q}_h \times \mathcal{V}_h \times \hat{\mathcal{V}}_h(0)$

$$ \int_{\mathcal{T}_h^F} \frac{\partial \epsilon_h}{\partial t} r \, d\Omega = - \int_{\mathcal{T}_h^S} \mathbf{v}_h \cdot \nabla r \, d\Omega + \int_{\partial \mathcal{T}_h^F} \hat{v}_h r \, d\Gamma, \qquad (61d) $$

$$ \int_{\mathcal{T}_h^F} \rho \frac{\partial \mathbf{v}_h}{\partial t} \cdot \mathbf{w} \, d\Omega = \int_{\mathcal{T}_h^S} \nabla \cdot (\kappa \epsilon_h) \cdot \mathbf{w} \, d\Omega - \int_{\partial \mathcal{T}_h^F} \boldsymbol{\tau}(\mathbf{v}_h \cdot \mathbf{n} - \hat{v}_h) \mathbf{w} \cdot \mathbf{n} \, d\Gamma, \qquad (61e) $$

$$ \int_{\partial \mathcal{T}_h^F} [\hat{\mathbf{p}}_h] \cdot \boldsymbol{\mu} \, d\Gamma = \int_{\Omega_S \cap \Gamma^N} \mathbf{t}^N \cdot \boldsymbol{\mu} \, d\Gamma. \qquad (61f) $$
On the fluid-solid interface, for all $\mu \in \hat{\mathcal{V}}_h(0)$

$$ \int_{\varepsilon_h^{SF}} \widehat{\sigma n}_h \cdot \mu \, d\Gamma = \int_{\varepsilon_h^{SF}} \hat{p}_h n \cdot \mu \, d\Gamma. \qquad (61g) $$

The equation (61c) encompasses the conservativity condition (22) on all the interior elastic faces with the weakly imposed Neumann conditions on $\Omega_F \cap \Gamma^N$. Similarly, (61f) encompasses the conservativity condition (25) with the weakly imposed Neumann condition on $\Omega_F \cap \Gamma^N$.

## 5.2 Relations between HDG and DG methods

We emphasize that the HDG variational formulation (61) is different from any DG formulation, for two reasons. First, the conservativity conditions are applied weakly, while in DG a pointwise relation is enforced at each element boundary node. Second, with HDG an hybrid variable is explicitly built and stored at element boundaries, while in DG no such variable is defined, but numerical fluxes are built at element boundaries, and not necessarily stored from one time step to the other. HDG may coincide after discretization with a DG method for some specific discretization choices but this is not the case generally. With our HDG approach, the weak conservativity conditions are actually pointwise since all the involved variables are in the same polynomial space. Moreover, due to the Riemann solvers equivalence (see Section 4.5), and due to our common discretization choices (see Section 5.4), our HDG method—after discretization—coincides with the DG method of Wilcox et al. (2010) in the elastic domain. By coincide we mean the two methods are algebraically equivalent although their implementation differs.

The same cannot be said for the acoustic domain since the strain representations differ. Moreover, when $h/p$ non-conforming interfaces are considered: classical DG solvers would project all the flux variables on the appropriate space (Bui-Thanh & Ghattas 2012) while HDG would just need a projection of the hybrid velocity (Cockburn et al. 2009). Therefore DG and HDG are no longer equivalent when non-conforming interfaces are considered.

## 5.3 Conservativity consistency and stability

The solution of (61) is conservative by construction, in the sense of Arnold et al. (2002). This is a common property of the HDG methods.

Using the intuitive definition given by Arnold et al. (2002), a DG method is consistent when the numerical traces computed from the exact solution are the traces of the exact solution. The HDG method proposed in this paper is consistent. It is a trivial consequence of the design of the hybrid unknowns which are the solutions of Riemann problems associated with sides states $q^-$ and $q^+$. When these states are equal to the exact solution, then there is no discontinuity at the interface. And the Riemann solution is then simply equal to the common value of both sides, that is, equal to the trace of the exact solution.

We follow Wilcox et al. (2010) and Nguyen et al. (2011b) to prove the semi-discrete stability. For the sake of simplicity, Neumann conditions are assumed everywhere on $\partial\Omega$. Let us first introduce the total mechanical discrete energy $E_h$ of the approximate HDG solution on the global domain $\Omega$, that is, defined as the sum of the kinetic energy and the potential energy on the elastic and acoustic parts of the domain

$$ E_h(t) = \int_{\Omega_S} \frac{1}{2} (\rho \mathbf{v}_h \cdot \mathbf{v}_h + \epsilon_h : \mathbf{C} : \epsilon_h) \, d\Omega + \int_{\Omega_F} \frac{1}{2} (\rho \mathbf{v}_h \cdot \mathbf{v}_h + \kappa \epsilon_h^2) \, d\Omega. \qquad (62) $$

Now, let us consider the HDG scheme (61), with the following test functions: in the elastic part, $r = \mathbf{C} : \epsilon_h(t)$ and $w = \mathbf{v}_h(t)$ for all $t \in \mathbf{I}$; in the acoustic part, $r = \kappa \epsilon_h(t)$ and $w = \mathbf{v}_h(t)$ for all $t \in \mathbf{I}$.

Summing (61a)–(61d), while using the above test functions, the left hand side is the time derivative of the total discrete energy $E_h(t)$. On the right hand side, the gradient terms can be all simplified, and the surface integration terms rearranged using the DtN and the conservativity conditions (61c), (61f) and (61g), in order to obtain the relation

$$ \frac{d}{dt} E_h(t) = - \int_{\partial T_h^S} (\mathbf{v}_h - \hat{\mathbf{v}}) \cdot \boldsymbol{\tau} (\mathbf{v}_h - \hat{\mathbf{v}}_h) \, d\Gamma - \int_{\partial T_h^F} (\mathbf{v}_h \cdot \mathbf{n} - \hat{\mathbf{v}}_h) \boldsymbol{\tau} (\mathbf{v}_h \cdot \mathbf{n} - \hat{\mathbf{v}}_h) \, d\Gamma \\ + \int_{T_h^S} \mathbf{f} \cdot \mathbf{v}_h \, d\Omega + \int_{\Omega_S \cap \Gamma^N} \boldsymbol{t}^N \cdot \mathbf{v}_h \, d\Omega + \int_{\Omega_F \cap \Gamma^N} t^N \mathbf{n} \cdot \mathbf{v}_h \, d\Omega. \qquad (63) $$

The last terms are the power associated with the external body force $\mathbf{f}$ and to the prescribed Neumann condition on the external boundary of the domain. The first terms on the right hand side are denoted $-\Theta_h$,

$$ \Theta_h = \int_{\partial T_h^S} (\mathbf{v}_h - \hat{\mathbf{v}}_h) \cdot \boldsymbol{\tau} (\mathbf{v}_h - \hat{\mathbf{v}}_h) \, d\Gamma + \int_{\partial T_h^F} (\mathbf{v}_h \cdot \mathbf{n} - \hat{\mathbf{v}}_h) \boldsymbol{\tau} (\mathbf{v}_h \cdot \mathbf{n} - \hat{\mathbf{v}}_h) \, d\Gamma. \qquad (64) $$

This term quantifies the energy dissipation introduced by the HDG approximation. It is worth to note that in the elastic part, $\boldsymbol{\tau}$ is always a symmetric definite positive matrix; and in the acoustic part $\boldsymbol{\tau}$ is always a positive scalar. Therefore $\Theta_h > 0$, implying that when no energy is supplied in the computational domain through an external source (i.e. a generalized body force or external boundary conditions) the total energy decays with time and the HDG scheme is stable.
Another interesting point is that the dissipation introduced by the HDG method is related only to some weighted norm of the velocity jumps at the interelement interfaces. The dissipation does not explicitly depend on the accuracy of the approximation inside the elements. Of course, the more accurate is the HDG approximation, the smaller are the jumps—the exact solution has no velocity jump—and the less dissipative is the HDG solution.

## 5.4 High-order polynomial approximation

We assume that each element $K \in T_h$ is the image of a reference element $\square = [-1, 1]^{n_d}$ through a differentiable local mapping $F_K$. This mapping is an invertible transformation from the reference coordinates $\xi = (\xi, \eta, \zeta)$ to the physical coordinates $\mathbf{x} = (x, y, z)$:

$$F_K: \square \rightarrow K \ (\xi, \eta, \zeta) \rightarrow F_K(x(\xi, \eta, \zeta), y(\xi, \eta, \zeta), z(\xi, \eta, \zeta)). \quad (65)$$

The HDG polynomial approximation is built using the $n_d$-tensorization of the 1-D Lagrange interpolation polynomials associated with the GLL nodes. This polynomial representation ensures a spectral convergence in term of the polynomial degree $N$ when the solution is sufficiently smooth (Deville et al. 2002; Canuto et al. 2012). Moreover, the GLL-grid involves the nodal values on mesh boundaries. As a $(n_d - 1)$-tensorization of the same Lagrange polynomials are used for the surface integrals, the connexions between the volumetric and surface fields are straightforward.

We use a classical collocation approach, extensively described in Kopriva (2009), and used by Wilcox et al. (2010); Bui-Thanh & Ghattas (2012) among many others. The volumetric (resp. surface) integrals of (61) are expressed on the volume (resp. faces) of the reference element $\square$. In order to handle the curved elements, the fluxes are expressed on a contravariant basis, and the surface integrals have to be weighted by metric terms.

Integrals are then evaluated using the GLL-quadrature, which for $N + 1$ points integrates exactly polynomials of order $\le 2N - 1$ (Deville et al. 2002). This choice leads to a diagonal element mass matrix, and significant simplifications when both surface and volumetric terms are discretized. This choice is commonly used by classical DG and CG-SEM (e.g. Komatitsch & Vilotte 1998; Wilcox et al. 2010).

Let us store the nodal variables in a vectorial form, with the contribution for each element $K \in T_h$ and each face $F \in \mathcal{E}_h$

$$\langle \epsilon(\mathbf{x}_i, t), \epsilon(\mathbf{x}_j, t) \rangle \rightarrow (\mathcal{E}), \quad \nu(\mathbf{x}_i, t) \rightarrow (V), \quad (\hat{\nu}(\mathbf{x}_i, t), \hat{v}_h(\mathbf{x}_i, t)) \rightarrow (\Lambda), \quad (66)$$

where the $\mathbf{x}_i$ are the indexed GLL nodes. After performing the numerical quadrature in (61), the following semi-discrete equations are obtained:

$$\begin{pmatrix} A & 0 \\ 0 & M \end{pmatrix} \frac{d}{dt} \begin{pmatrix} \mathcal{E} \\ V \end{pmatrix} = \begin{pmatrix} 0 & -B_v \\ B_\epsilon & -E_v \end{pmatrix} \begin{pmatrix} \mathcal{E} \\ V \end{pmatrix} + \begin{pmatrix} C_\lambda \\ E_\lambda \end{pmatrix} \Lambda + F_{\text{ext}}, \quad (67a)$$

$$0 = C_{\epsilon} \mathcal{E} - B_{\lambda}^{T} V + G \Lambda. \quad (67b)$$

Both matrix $A$ and the mass matrix $M$ are diagonal by construction since the quadrature nodes are also the interpolation nodes. Matrices $B_v$ and $B_\epsilon$ correspond to the integrals with a gradient term. Due to the discontinuous nature of spaces $Q_h$, $V_h$ and $S_h$, the matrix $[0 - B_v; B_\epsilon - E_v]$ is block-diagonal. This matrix basically couples all the degrees of freedom that are in the same element. Matrices $C_\lambda$, $E_\lambda$ and $C_\epsilon$ couple the elements degrees of freedom with the hybrid ones located on the adjacent faces. Finally, the matrix $G$ is the discretization of terms involving $(\hat{\nu}, \hat{v}_h)$ in the conservativity conditions, that is,

$$\int_{\partial T_h^S} (\tau^+ + \tau^-) \hat{\nu} \cdot \mu \, d\Gamma + \int_{\partial T_h^F} (\tau^+ + \tau^-) \hat{v}_h \, \mu \, d\Gamma + \int_{\mathcal{E}_h^{SF}} (\tau^+ + \tau^-_{ac}) \hat{\nu} \cdot \mu \, d\Gamma \Rightarrow G\Lambda. \quad (68)$$

Thanks to the discontinuous nature of $\hat{V}_h$ and $\hat{V}_h$, the matrix $G$ is block-diagonal, with a block per mesh face. Moreover, thanks to the GLL quadrature rule, the degrees of freedom lying on different nodes of the same face are uncoupled. But the $n_d$ components of $\hat{v}$ located on a particular node are coupled together. Therefore, $G$ has a block-diagonal structure in the elastic domain, and is diagonal on the acoustic domain. The block sizes are at most $n_d \times n_d$, and they inherit their symmetric positive definite property from the matrix $\tau$. This point is important since an efficient $LL^T$ factorization of $G$ can be computed for its inversion.

The ODEs (67a) are the so-called *local solvers*: they provide an update of the strain and velocity fields in each element $K \in T_h$. Indeed when the hybrid variable $\Lambda$ are known, the element local solvers can be performed independently on each element. The algebraic equation (67b) is the discrete conservativity condition. It is the only *global* equation, coupling explicitly interelement contributions. In a typical HDG resolution strategy, global equation (67b) is solved first, and then local solvers (67a) are applied.

## 6 TIME INTEGRATION

Let us now introduce the following notations

$$\mathbf{y}(t) = \begin{pmatrix} \mathcal{E} \\ V \end{pmatrix}, \quad f(\mathbf{y}, \Lambda) = \begin{pmatrix} A & 0 \\ 0 & M \end{pmatrix}^{-1} \left[ \begin{pmatrix} 0 & -B_v \\ B_\epsilon & -E_v \end{pmatrix} \begin{pmatrix} \mathcal{E} \\ V \end{pmatrix} + \begin{pmatrix} C_\lambda \\ E_\lambda \end{pmatrix} \Lambda + F_{\text{ext}} \right], \quad (69)$$
together with the function g

$$g(y, \Lambda) = \mathbb{C}_e \mathcal{E} - \mathbb{E}_{\lambda}^{T} V + \mathbb{G} \Lambda. \qquad (70)$$

The system (67) can then be rewritten as

$$\frac{dy(t)}{dt} = f(y, \Lambda), \qquad (71a)$$

$$0 = g(y, \Lambda). \qquad (71b)$$

These semi-discrete equations thus define a set of ADEs, which states that the approximate solution should evolve in time staying on the manifold defined by the discrete conservativity condition (67b).

Following Hairer (2011) and Hairer et al. (2010), the ADE system (71) can be rewritten as an ODE system if $\frac{\partial g}{\partial \Lambda}$ is invertible. It is the case for elastic-acoustic wave problems as $\frac{\partial g}{\partial \Lambda} = \mathbb{G}$ is always invertible, being a block diagonal symmetric definite positive matrix.

Then the relation (71b) can be explicitly inverted for $\Lambda$

$$\Lambda = h(y) = \mathbb{G}^{-1} (\mathbb{E}_{\lambda}^{T} V - \mathbb{C}_e \mathcal{E}). \qquad (72)$$

This last relation can now be used to replace $\Lambda$ in equation (71a)

$$\frac{dy(t)}{dt} = f(y, h(y)), \qquad (73)$$

which defines a classical system of ODEs, and then any classical ODE integrator may be used, either implicit or explicit. This is an important practical point allowing to use the same numerical time integration scheme for both the HDG-SEM and the CG-SEM formulations. However, equation (73) should not be misleading in the sense that it is presented here just for the purpose of justifying the use of classical ODE integrators. Equation (71) remains the actually solved one. $\Lambda$ is computed and stored as such, and is not eliminated as would occur in the DG approach.

We choose to consider explicit time schemes for two reasons. First, we want to take advantage of the diagonal structure of the matrices $M$ and $\Lambda$ and of the fast inversion of $\mathbb{G}$. Second, the parallelization is more straightforward. In our analysis both a second-order midpoint and an explicit fourth-order Runge-Kutta have been used.

## 6.1 Second-order midpoint method

The midpoint integration scheme basically performs the force equilibrium at time $t^{n+1/2}$

$$\Lambda^{n+1/2} = h(y^{n+1/2}), \qquad (74a)$$

$$y^{n+1} = y^n + \Delta t f(y^{n+1/2}, \Lambda^{n+1/2}). \qquad (74b)$$

The mid-time step values are evaluated using a multi-corrector implementation, that is, with successive approximations ($\Lambda_{[i]}^{n+1/2}$, $y_{[i]}^{n+1/2}$) initialized by a first extrapolation step.

For most applications, only one iteration is actually needed. This integration scheme is second order accurate, and stable for CFL numbers less or equal to 0.8. It is almost energy conservative, in the sense that the energy dissipation comes from the HDG discretization itself (63) and not from the time scheme.

## 6.2 Fourth-order explicit Runge-Kutta method

A fourth-order Low Storage Explicit Runge-Kutta (LSERK4) method (Carpenter & Kennedy 1994; Hesthaven & Warburton 2008) has also been implemented in this study. Only one intermediary state needs to be stored in memory: $\Delta y_{[i]}$, but at the price of one extra computational stage. The LSERK4 method is implemented as algorithm 1.

<caption>**Algorithm 1 LSERK4 Time Integrator**</caption>

```
y_{[0]}^{n+1} = y^n

for i = 1..5 do
    $\Lambda_{[i-1]}$ = h(y_{[i-1]})
    $\Delta y_{[i]}$ = a_i $\Delta y_{[i-1]}$ + $\Delta t$ f(y_{[i-1]}^{n+1}, $\Lambda_{[i-1]}$), t^{n+c_i})
    y_{[i]}^{n+1} = y_{[i-1]}^{n+1} + b_i $\Delta y_{[i]}$

end for
$\Lambda^{n+1} = h(y^{n+1})$
```

The coefficients $a_i$, $b_i$, $c_i$ are given in Carpenter & Kennedy (1994). Note that $t^{n+c_i}$ means $t^n + c_i \Delta t$, with $0 \le c_i \le 1$.
## 6.3 Time scheme selection

To assess the relative performances of the two time schemes, let us introduce the Courant number as $n_C = \max[c\Delta t/\delta x]$ where the $\delta x$ is the minimum distance between two GLL nodes. Although the midpoint algorithm is stable for Courant number up to $n_C = 0.8$, a maximum value of $n_C \approx 0.2$ or $n_C \approx 0.25$ is actually imposed to minimize the dispersion error. LSERK4 allows $n_C$ up to 1, with usually good dispersion properties for high $n_C$. The LSERK4 is more accurate and significantly less expensive than the midpoint scheme. However the midpoint may be favoured when emphasis is put on the energy conservation.

# 7 NUMERICAL RESULTS

In this section, the accuracy of the method is illustrated by several numerical benchmarks. So far the HDG-SEM has been implemented only in 2-D, all the results presented thereafter are given for in-plane configurations. Except when LSERK4 is explicitly mentioned, the midpoint scheme is used for the numerical applications presented in this section.

Absorbing layers, based on the ADE-PML formulation from Gedney & Zhao (2010) have been adapted to the HDG framework following the stress-velocity formulation proposed by Martin et al. (2010). The auxiliary differential equations are spatially discretized with the HDG methodology, and integrated in time with the scheme used for the inner computational domain. Just one layer of PML elements is used. Their thickness in the direction of absorption is by default twice the size of the neighbouring inner elements.

## 7.1 Lamb's problem

The Lamb's problem (Lamb 1904) is a classical benchmark for elastic wave propagation codes, checking the accuracy of the surface waves modelling and whether the free surface condition is correctly modelled. In this problem, a vertical point force is applied at the surface of a linear isotropic elastic half-space. This force generates direct P-waves and S-waves together with a non-dispersive Rayleigh surface wave propagating along the free surface.

For this problem, a semi-infinite homogeneous isotropic elastic media is considered with a P-wave velocity $c_p = 3200 \text{ m s}^{-1}$, an S-wave velocity $c_s = 1700 \text{ m s}^{-1}$ and a mass density $\rho = 1700 \text{ kg m}^{-3}$. The simulation domain is truncated to be a rectangular box with dimensions $x \in [-5000, 5000] \text{ m}$, $y \in [-4000, 0] \text{ m}$, and the final time is $T = 5 \text{ s}$. The free surface is an horizontal line situated on the top $y_{\max} = 0 \text{ m}$. ADE-PMLs are set at the left, bottom and right sides of the domain. The source is a point force directed towards the bottom $f = (0, -1)^T f(t)$ and applied just beneath the surface at $x_s = 0$, $y_s = -20 \text{ m}$. The function source is defined as $f(t) = A G(t)$ where $G(t)$ is a Ricker wavelet, with a central frequency of 3 Hz and a time delay of 0.5 sec, and an amplification factor $A = 10^{10}$.

Accuracy of the numerical solution is assessed against the analytical solution based on a Cagniard de Hoop method, implemented in the EX2DDIR code (Berg et al. 1994; De Hoop 1960). The $L^2$-norm of the velocity error is computed at receivers along the surface for $0 \le t \le 2.8 \text{ s}$ such that the incomplete PML absorption does not have time to pollute the recorded signal. The convergence analysis is performed as follows: for a given polynomial order $N$, simulations are performed for a series of increasingly finer uniform quadrilateral meshes, with stretched elements in the absorbing PML. The investigated element sizes $h$ are: 400, 200, 100, 50, 25 and 12.5 m. Results are log-plotted against the inverse of the mesh size $1/h$, and are shown in Fig. 9 for a receiver located at position $x_r = 2000 \text{ m}$, $y_r = 0 \text{ m}$.

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Lamb's problem: errors and convergence orders observed. Left: $L^2$-errors for the numerical solution $v$ and the post-processed solution $v^*$ for different polynomial degrees $N$. Right: estimated orders of convergence for different $N$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Elastic-acoustic test configuration. The magnitude of the velocity fields are shown at $t = 1$ s (left) and $t = 1.8$ s (right). In the acoustic domain, the direct (a) and the reflected P-wave (b) are observed. In the elastic domain, the transmitted P-wave (c) and the converted P-to-S wave (d) are the most energetic waves. Less energy is conveyed into refracted waves in the elastic domain (e) and in the acoustic domain (f). The interface Scholte wave (g) propagates along the interface and is almost detached from the faster S-wave (d).</figcaption></figure>

Estimated orders of convergence shown in Fig. 9 are computed, for each element order $N$, from the three most precise points (except with $N = 10$ where only two points have been used). Note that the Courant number is adapted such that the error due to the time integration is never dominant. Thus $n_c$ varies between 0.0025 and 0.25.

The numerical velocity solution $v_h$ converges with the order $N + 1$ for all the polynomial orders considered in this study. A slightly higher degree of convergence of almost $N + 3/2$ can be observed in some cases. In particular, the convergence rate observed for $N = 2$ is surprisingly good. Those results are consistent with observations from Nguyen *et al.* (2011b).

Following the post-processing methodology presented in Appendix A, one can compute a post processed velocity $v_h^*$ that always exhibits higher accuracy and convergence rate. However no clear convergence order can be determined. While optimal superconvergence order $N + 2$ is observed for $N = 2$ and $N = 5$, the order is suboptimal for $N = 8$ and $N = 3$. These observations are consistent with Fu *et al.* (2015).

## 7.2 Elastic-acoustic interface

In this problem, we consider the wavefield generated by a point source located close to a horizontal interface, at $y = 0$ m, separating two homogeneous infinite half-spaces. The half-space below the interface is a homogeneous linear isotropic medium, while the half-space above the interface is a homogeneous acoustic medium.

The acoustic medium is characterized by a P-wave velocity $c_p = 1500 \text{ m s}^{-1}$, and a mass density $\rho = 1020 \text{ kg m}^{-3}$. The elastic medium is characterized by a P-wave velocity $c_p = 3600 \text{ m s}^{-1}$, an S-wave velocity $c_s = 2000 \text{ m s}^{-1}$ and a mass density $\rho = 2500 \text{ kg m}^{-3}$. The source is a pure explosion, that is, a spherical moment, and is located in the acoustic medium at the position $x_s = 0, y_s = 650 \text{ m}$. The source has an amplification coefficient $A = 1 \times 10^{11}$ and the time dependence is a Ricker wavelet with a central frequency of 5 Hz and a time delay of 0.5 s. The simulation domain is truncated as a finite rectangular box with dimensions $x \in [-4000, 4000] \text{ m}, y \in [-3500, 2500] \text{ m}$, and the duration of the simulation is $T = 4 \text{ s}$. Absorbing layers (ADE-PML) are set at the four end sides of the domain.

This test illustrates the accuracy of the HDG-SEM in simulating the conversion of the spherical pressure wave at the acoustic-elastic interface. When the source lies in the acoustic part, the pressure wave is partly reflected at the interface, while energy transmission at the interface generates both P-waves and S-waves in the elastic media together with an interface wave travelling along the acoustic-elastic interface.

The reference analytical solution, based on a Cagniard-de Hoop method (Cagniard 1962; De Hoop 1960), is implemented in the Gar6more-2D code provided by Diaz (2005). The convergence analysis is performed for different element sizes and element polynomial orders $N$. The mesh size $h$ is refined with a factor $1/\sqrt{2}$, and mesh sizes analysed here are: 400, 282.84, 200, 141.42, 100, 70.71, 50, 35.35, 25 and 12.5 m. The Courant number is adapted, and $n_c$ varies between 0.002 and 0.25. In Fig. 11, the $L^2$-norm of the velocity error recorded at a receiver in the acoustic domain ($x_r = 2000 \text{ m}, y_r = 550 \text{ m}$) is log-plotted against the inverse of the mesh size $1/h$, together with the observed order of convergence. Note that the $L^2$-norm of the velocity error is computed for $0 \le t \le 2.5 \text{ sec}$, such that the incomplete PML absorption does not affect the convergence results.

Regarding the velocity approximation $v_h$, the expected $N + 1$ order of convergence is observed for all the orders of element polynomial approximation considered in this study, in agreement with Cockburn & Quenneville-Bélair (2014), Nguyen *et al.* (2011b) and Stanglmeier *et al.* (2016). A higher order of convergence is surprisingly observed for the fourth order polynomial approximation.

The post-processed pressure, (see Appendix A), super-converges with an optimal $N + 2$ order as expected.
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Elastic-acoustic test case: the $L^2$-norm of the velocity error (left) and the estimated orders of convergence (right) are shown for different order of element polynomial approximation $N$.</figcaption></figure>

<table><thead><tr><th>Polynomial degree <i>N</i></th><th>Order ||<i>v</i> - <i>v</i><sub><i>h</i></sub>||</th><th>Order ||<i>p</i> - <i>p</i><sub><i>h</i></sub><sup>*</sup>||</th></tr></thead><tbody><tr><td>2</td><td>3.15</td><td>3.92</td></tr><tr><td>3</td><td>3.95</td><td>4.88</td></tr><tr><td>4</td><td>5.48</td><td>6.40</td></tr><tr><td>5</td><td>5.96</td><td>6.89</td></tr><tr><td>6</td><td>7.08</td><td>8.02</td></tr><tr><td>8</td><td>9.36</td><td>9.97</td></tr><tr><td>10</td><td>10.95</td><td>11.92</td></tr></tbody></table>

## 7.3 Elements of comparison between HDG, DG and CG-SEM

It is always a difficult task to compare the efficiency of different numerical methods. Such comparisons are often problem-dependent and always implementation-dependent. We however want to give here an idea about the relative computational cost of our HDG method when compared with the CG-SEM and the DG method of Wilcox et al. (2010).

From the previous convergence studies, let us consider both the Lamb and the elastic-acoustic configurations with $h = 100$ m and $N = 5$. For these settings, all methods show approximately the same accuracy, and their performances are summarized in Table 1. All the schemes have been time integrated with the LSERK4 algorithm setting $n_c = 0.2$ small enough such that the error in space is clearly dominant. We use a displacement-velocity formulation for the CG method in the elastic domain, but unfortunately we do not have a velocity-potential based CG method available for the acoustic domain. A single core Intel Xeon CPU E5-2660, 2.60 GHz has been used for all the simulations.

For the Lamb's problem, it can be noted that the CG-SEM is around 31 per cent faster than HDG, which is expected due to the extra computational burden of the hybrid unknowns. DG are around 16 per cent slower than HDG. This may be due to some higher complexity of the DG flux computation (see Section 4.5) or to a suboptimal implementation of that DG method. It is worth to note that DG and HDG give here almost exactly the same velocity outputs (with a relative difference around 1.e-12) confirming the equivalence DG-HDG mentioned previously. Finally, although HDG is slower than CG, its accuracy may be significantly improved thanks to a low-cost post-processing—performed only for the element where the receiver lies.

For the elastic-acoustic problem, both DG and HDG approximate the velocity and the pressure with the same accuracy. This time the DG approach is much slower since it uses a full second order $\epsilon$ tensor in the acoustic domain, while HDG makes use of a scalar $\epsilon$ instead. However the DG approach could be certainly more efficient if reformulated with a pressure-velocity acoustic formulation. Note that for the DG simulation, the pressure is retrieved using $p = -\kappa(\epsilon_{11} + \epsilon_{22})$.

The authors warn that none of the implementations used here have been seriously optimized, so the CPU-time results should not be overstated. Moreover, the relative overcosts of DG and HDG decrease with the polynomial order.

## 7.4 Curved interfaces

The possibility of using curve-shaped elements is important to capture interface geometric features. Even though arbitrary high-order local geometrical transformations can be used, as long as conditions mentioned by Kopriva (2006) are satisfied, quadratic local geometrical transformations are often used in practice.

<table><caption>**Table 1.** Comparisons of performances and accuracy for various numerical methods, on both Lamb and elastic-acoustic problems. The symbol (*) means that post-processing has been performed.</caption><thead><tr><th></th><th>CG-SEM</th><th>DG-SEM</th><th>HDG-SEM</th><th>HDG-SEM (*)</th></tr></thead><tbody><tr><td>Lamb</td><td></td><td></td><td></td><td></td></tr><tr><td>CPU time (s)</td><td>1226.8</td><td>2074.5</td><td>1786.4</td><td>1788.6</td></tr><tr><td>L<sub>2</sub>-error on <i>v</i></td><td>4.01e-4</td><td>3.60e-4</td><td>3.60e-4</td><td>2.76e-4</td></tr><tr><td>Acou-Elas</td><td></td><td></td><td></td><td></td></tr><tr><td>CPU time (s)</td><td>–</td><td>2114.6</td><td>1422.2</td><td>1423.2</td></tr><tr><td>L<sub>2</sub>-error on <i>v</i></td><td>–</td><td>4.28e-3</td><td>4.28e-3</td><td>–</td></tr><tr><td>L<sub>2</sub>-error on <i>p</i></td><td>–</td><td>6.55e3</td><td>6.55e3</td><td>4.09e3</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** On the left, configuration of the two-layer computational domain, with the localizations of the explosive source and of the line of receivers. On the right, typical snapshot of norm of the velocity at $t = 0.8$ s, showing the direct P-wave (a), and the S-wave (b) converted at the free surface. The Rayleigh surface wave (c) is propagating along the curved free surface with strong Rayleigh-to-body-wave conversion (d) due to the geometry of the free surface.</figcaption></figure>

The test configuration is inspired from Komatitsch & Vilotte (1998). The computational domain is a two-layer linear elastic half-space, with strong contrast in both wave speeds and Poisson's ratios. The free surface and the material interfaces have a smooth topography, that is, the material interface has a typical dome shape, while the free surface has an uncorrelated topography. The top layer is characterized by $\rho = 1000 \text{ kg m}^{-3}$, $c_p = 2000 \text{ m s}^{-1}$, $c_s = 1300 \text{ m s}^{-1}$, while the lower layer is characterized by $\rho = 1500 \text{ kg m}^{-3}$, $c_p = 2800 \text{ m s}^{-1}$, $c_s = 1473 \text{ m s}^{-1}$.

The computational domain is a box with $x \in [0, 2500]$ m, and $y \in [0, \sim 1700]$ m, spatially discretized with a conforming mesh of 48×40 elements, and an element polynomial order of 9. Thus, the number of points per minimal wavelength is at least 5. The explosive point source, that is, spherical moment tensor, is located at $x_s = 620 \text{ m}$, $y_s = 1630 \text{ m}$, with a Ricker source time function of central frequency 15 Hz. An horizontal line of 100 receivers is located within the top layer, with a spacing of 14 m, from position (700, 1300) m to (2100, 1300) m, as shown on Fig. 12. The Courant number is set to $n_c = 0.2$ and the duration of the simulation is $t = 2.5$ s.

In the snapshot shown on Fig. 12, a converted surface-to-body-wave can be observed, which is generated during the propagation of the Rayleigh wave along the free surface topography, and which is partially superimposed on the S-wave generated by the conversion of the direct P-wave at the free surface. This phenomenon is consistent with the results of Komatitsch & Vilotte (1998) using CG-SEM.

A comparison between CG-SEM and HDG-SEM is shown on Fig. 13, for the same computational domain, mesh discretization, element polynomial approximation, and time step. For HDG-SEM, the solution has been post-processed according to A. CG-SEM uses a second-order Leapfrog time integration scheme, while HDG-SEM uses a second-order midpoint scheme. Both numerical methods produce very similar results.

## 7.5 Coupling CG-SEM with HDG

As shown by Cockburn et al. (2007, 2009), CG methods can also be hybridized. In particular, for the CG-SEM approximation of the elastic wave equations, this leads to the classical formulation of the elastic wave equation as a velocity-displacement first-order hyperbolic system (Festa & Vilotte 2005). Hybridization results of CG methods provides a natural framework for coupling CG-SEM and DG-SEM in a common computational domain in a very simple way. Let us consider an interface $\Gamma_{CD}$ separating a domain $\Omega_{CG}$ where the problem is approximated by CG-SEM from a domain $\Omega_{DG}$ where the problem is approximated by HDG-SEM. The coupling relations can therefore be simply formulated as in Cockburn et al. (2009); Perugia & Schötzau (2001):

$$ \hat{v}_h = v|_{\Omega_{CG}}, \qquad (75a) $$

$$ \widehat{\sigma \bar{n}}_h = -(\sigma n)|_{\Omega_{DG}} + \tau (v|_{\Omega_{DG}} - \hat{v}_h). \qquad (75b) $$

One must pay attention to the fact that from the CG side, the hybrid field needs to be continuous on the mesh skeleton.

A simple test case configuration is considered here, with a linear isotropic heterogeneous elastic medium of dimensions $x \in [-2500, 2500]$ m and $y \in [-2500, 2500]$ m. A planar interface located at $y = -1030 \text{ m}$ is separating the domain into an upper domain where the solution is approximated using CG-SEM and a lower domain where the solution is approximated by HDG-SEM. Reflection boundary conditions are prescribed on the external boundary. The upper elastic domain is characterized by $\rho = 1700 \text{ kg m}^{-3}$, $c_p = 3200 \text{ m s}^{-1}$ and $c_s = 1700 \text{ m s}^{-1}$. The lower elastic domain is characterized by $\rho = 800 \text{ kg m}^{-3}$, $c_p = 2400 \text{ m s}^{-1}$ and $c_s = 1200 \text{ m s}^{-1}$.

The domain is discretized using a uniform quadrangular mesh, with 51 × 51 elements, and an element polynomial order of 4. The source is a horizontal point force located at the centre of the domain, with an amplification factor equal to $10^{10}$ and a Ricker source time function of central frequency 3 Hz and time delay 0.5 s. For both the CG-SEM and HDG-SEM the explicit LSERK4 time scheme is used, with the same time step. The simulation duration is 10 s, with 3780 time iterations and $n_c = 0.5$, allowing waves to be transmitted across the interface
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Top: time response of the horizontal and vertical components of the displacement at the horizontal line of receivers located within the top layer as obtained with the HDG-SEM method. The direct *P*-wave (a) can be clearly distinguished, as well as a strong Rayleigh-to-body-wave (d) generated by the free surface curvature. This wave travels at the *S*-wave speed and cannot always be clearly distinguished from the main *S*-wave in this figure. Bottom: horizontal and vertical velocities recorded at receiver 50 for the HDG-SEM method, together with the error residual when compared with the CG-SEM solution.</figcaption></figure>

several times after multiple reflections at the domain boundaries. The coupling condition (75) are applied along the planar interface, and are solved here explicitly using the LSERK4 scheme.

Particle velocities are recorded in the HDG-SEM part of the computational domain at a receiver located at $x_r = 1500$ m, $y_r = -2000$ m. Results are shown on Fig. 14. Results are assessed by comparing with a reference solution obtained using CG-SEM in the whole domain. As observed on Fig. 14, the coupling between CG-SEM and HDG-SEM is quite accurate during the whole simulation.

## 7.6 Infrasonic waves emitted by a seismic event

Seismic events can generate infrasonic waves that may be ducted by atmospheric waveguides over regional and even global scales. Infrasound signals have been recorded and related to seismic sources and surface topography features such as mountains (e.g. Walker *et al.* 2013, for the Tohoku earthquake), or cliffs (Green *et al.* 2009), even for low magnitude earthquakes.
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Upper left: test configuration for the CG-SEM and HDG-SEM coupling. Upper right: snapshot of the velocity norm at time $t = 1.5$ s. Lower left and right: velocity components recorded at the receiver located in the HDG domain. The amplified difference between the numerical solution and the CG-SEM reference solution is shown in blue. The difference remains small even after several wave transmissions across the coupling interface.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Left: numerical domain for the seismic-infrasonic simulation. Right: vertical evolution of the sound speed $c_p$ at the origin of the atmospheric waveguide.</figcaption></figure>

Seismic-infrasonic coupling involves a variety of wave phenomena, such as elastic-to-acoustic wave conversion associated to surface waves propagating along the Earth's surface, wave diffraction by small topographic features and infrasound channelization in atmospheric waveguides for example.

A simple test case configuration is considered here. A homogeneous linear elastic half-space is coupled with an acoustic half-space by a piecewise horizontal interface that includes a 50 m high cliff. The lower elastic domain is characterized by $c_p = 3200 \text{ m s}^{-1}$, $c_s = 1700 \text{ m s}^{-1}$ and $\rho = 1700 \text{ kg m}^{-3}$. In the upper acoustic domain, a vertical variation of the sound speed is introduced as a simplified representation of the observed stratospheric structure, see Fig. 15, with near the ground a sound speed $c_p^0 = 330 \text{ m s}^{-1}$ and a mass density $\rho^0 = 1.3 \text{ kg m}^{-3}$. This variation is imposed through an inhomogeneous $\kappa(x)$ in our code.
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Snapshots of velocity norm for a 75 m source depth, at times $t = 2.5$, $t = 8.5$, and $t = 29$ s. Elastic P-wave (a), and Rayleigh waves (b1) and (b2) generate infrasonic waves (c), (d1), (d2), (e), (f) and (g). Only (e), (f) and (g) are ducted in the atmospheric waveguide, and their multiple refractions/reflections can be clearly seen at $t = 29$ s.</figcaption></figure>

The simulation domain is truncated as a rectangular box of dimensions $x \in [0, 20000]$ m, $y \in [-2500, 2500]$ m. The Earth's surface is at $y = 50$ m on the left side of the cliff, and at $y = 0$ m on the right side. The 50 m cliff is located at $x = 5000$ m. Absorbing layers (ADE-PML) are used at the four sides of the domain. A very shallow source, at 75 m depth, is considered in the elastic domain located at 2500 m from the left side of the domain. The explosive source time function is a Ricker wavelet with a 2 Hz central frequency.

The computational domain is spatially discretized using 400 × 100 quadrangular elements, with an element polynomial order equal to 4, such that at least 6 points sample the minimal wavelength in the atmosphere. The ADE-PML layer is 100 m thick and its element polynomial order is set to 8 along the absorbing directions. The minimum propagating wavelength is $\lambda_{min} = 227$ m, and therefore the depth of the seismic source is much smaller than the minimum propagating wave length. Snapshots of the velocity norm are recorded at $t = 2.5$ s, $t = 8.5$ s and $t = 29$ s, and are reported on Fig. 16.

In this configuration, a complex wavefield is generated. In the solid domain, a direct P-wave (a), a converted S-wave and a surface Rayleigh wave (b1) are generated around the shallow source. A second Rayleigh wave (b2) results from the diffraction of the direct P wave at the cliff. In the acoustic domain, infrasonic waves can be classified into three families. First, refracted conical waves (c), (d1) and (d2) that are generated by the direct P-wave (a) and Rayleigh surface waves (b1) and (b2), just like head waves can be. Due to the strong contrasts at the interface between seismic and acoustic velocities, wave fronts of (c), (d1) and (d2) are almost horizontal. Therefore these last three waves are not ducted in the atmospheric waveguide. They are tightly related to seismic wave fronts at the solid-air interface, and do not carry information about the atmospheric structure. The second family of infrasonic waves is composed of waves (e) and (f) respectively generated by diffraction of the P-wave and the Rayleigh wave by the cliff. Finally, the (g) wave is a pressure wave directly transmitted at the Earth's surface on top of the seismic source. In contrast with the first family, waves (e) (f) and (g) have all a spherical-shaped wave front, and a significant part of their energy can therefore be ducted in the atmospheric waveguide. All these waves clearly appear when the source is shallow enough. Given the sound speed gradient in the atmosphere, infrasonic phases (e), (f) and (g) are refracted enough to partially impinge back on the solid/air interface, where they are reflected, and so on. The $t = 29$ s snapshot on Fig. 16 is a typical illustration of this phenomenon.

This simple example illustrates some capabilities of the HDG-SEM method in dealing with elastic-acoustic coupling in heterogeneous geological configurations.

## 7.7 Sismo-hydroacoustic coupling

Finally, a sismo-hydroacoustic coupling test case is considered derived from a real world application. The computational domain is a simplified 2-D cross-section of the French Antilles region that is 350 km long and 60 km deep. A part of the top section is filled with the oceanic domain,
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Top: schematic representation of the Antilles section studied. Below: three successive snapshots showing the increasing complexity of the wave field.</figcaption></figure>

with a maximum depth of 4 km. The oceanic domain is on top of an Earth's crust, which is approximately 16 km deep. Finally, the Earth's mantle occupies the lower part of the domain, between 16 and 60 km deep.

Realistic acoustic and elastic parameters have been prescribed. In the ocean, the SOFAR channel is modelled with a variation of the acoustic wave speed between 1500 and 1600 m s⁻¹. In the crustal layer, the mass density is ρ = 2700 kg m⁻³, and varying (λ(x), μ(x)) are applied so that the elastic velocities increases linearly with the depth. cₚ varies from 4500 at the top to 7000 m s⁻¹ at the bottom, and cₛ = cₚ/√₃. For the mantle, ρ = 3300 kg m⁻³, cₚ = 8000 m s⁻¹ and cₛ = 4600 m s⁻¹. A point source is located in the crust, with a maximum frequency around 5 Hz.

The domain is spatially discretized with 1.4 × 10⁵ quadrangular elements, using fourth-order polynomial approximation, such that at least 5 or 6 points sample the minimal wavelength of the hydroacoustic waves. This represents almost 20 × 10⁶ degrees of freedom. The duration of the simulation is T = 50 s. The LSERK4 time scheme is used here with n<sub>c</sub> = 1, leading to 17 000 time steps for the whole simulation.

The complexity of the generated wave field can clearly be observed on the velocity snapshots of Fig. 17. Among all the phases, the transmitted *P* wave (a) to the mantle, the fastest wave, can be observed. Its conversion at the mantle-crust interface is associated with the arrival of a precursor *P*-refracted wave (a1) that would be recorded first on a remote seismogram on shore. The shape of the wave front (b) results from the effect of the velocity gradient in the crust, deviating the energy upwards. The wave front (c) is associated with the reflected *S*-wave at the seabed. It is quite energetic due to the low energy transmission of an almost normal incident *S*-wave at the elastic-acoustic interface. Finally, (e) is associated with the train of acoustic *T*-waves ducted in the ocean by multiple reflections at the free surface and the seafloor.

This simulation is rather demanding, that is, the highest frequencies propagate over 250 wavelengths, but can still be run on a single core of an Intel Xeon E5-2650, 2.0 GHz processor.
## 8 CONCLUSION

We have presented an HDG-SEM method for solving elastic-acoustic wave propagation in heterogeneous media. This HDG-SEM makes use of a global conservativeness condition through local DtN operators that can be explicitly constructed using an hybridization of exact Riemann solvers at the element interfaces. The resulting semi-discrete system of algebraic differential equations can be reformulated as a classical system of ODEs allowing efficient use of classical explicit and implicit time schemes. The computational cost associated with this method, seems intermediate between those of classical CG and DG methods.

The HDG-SEM inherits attractive properties from spectral element methods. It has an arbitrary high order of accuracy in space, and its computational efficiency increases with the polynomial order. It also inherits the interesting superconvergence property of HDG methods, and it provides a natural framework for coupling efficiently CG-SEM and HDG-SEM in an unique computational domain. As shown by Cockburn et al. (2009), the HDG methods add new elegant *h/p*-adaptive flexibility properties that still need to be further investigated for geophysical wave propagation applications.

The HDG-SEM formulation has attractive properties for efficient parallelization in new massively parallel hybrid architectures. We thus consider the HDG-SEM quite attractive for real world geophysical wave applications involving large scale problems with complex geometries and heterogeneities.

## ACKNOWLEDGEMENTS

The authors want to thank the CEA for funding this study. The numerical computations were partly performed on the S-CAPAD IPGP platform, with the support of G. Moguilny and F. Houssen.

## REFERENCES

Ainsworth, M., Monk, P. & Muniz, W., 2006. Dispersive and dissipative properties of discontinuous Galerkin finite element methods for the second-order wave equation, *J. Sci. Comput.*, **27**(1–3), 5–40.

Arnold, D.N. & Brezzi, F., 1985. Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates, *RAIRO-Modélisation mathématique et analyse numérique*, **19**(1), 7–32.

Arnold, D.N., Brezzi, F., Cockburn, B. & Marini, L.D., 2002. Unified analysis of discontinuous Galerkin methods for elliptic problems, *SIAM J. Numer. Anal.*, **39**(5), 1749–1779.

Berg, P. et al., 1994. Exploration oriented seismic modeling and inversion, in *Modeling the Earth for Oil Exploration*, pp. 417–524, ed. Helbig, K., Pergamon, Amsterdam.

Bonnasse-Gahot, M., 2015. High order discontinuous Galerkin methods for time-harmonic elastodynamics, *PhD thesis*, Université Nice Sophia Antipolis.

Bottero, A., Cristini, P., Komatitsch, D. & Asch, M., 2016. An axisymmetric time-domain spectral element method for full-wave simulations: application to ocean acoustics, *J acoust. Soc. Am.*, **140**, 3520.

Brezzi, F. & Fortin, M., 1991. *Mixed and Hybrid Finite Element Methods*, Springer-Verlag.

Brezzi, F., Douglas, J. & Marini, L., 1985. Two families of mixed finite elements for second order elliptic problems, *Numer. Math.*, **47**(2), 217–235.

Bui-Thanh, T., 2015. From Godunov to a unified hybridized discontinuous Galerkin framework for partial differential equations, *J. Comp. Phys.*, **295**(C), 114–146.

Bui-Thanh, T. & Ghattas, O., 2012. Analysis of an hp-nonconforming discontinuous Galerkin spectral element method for wave propagation, *SIAM J. Numer. Anal.*, **50**(3), 1801–1826.

Cagniard, L., 1962. *Reflection and Refraction of Progressive Seismic Waves*, McGraw Hill Book Company.

Canuto, C., Hussaini, M.Y., Quarteroni, A.M. & Thomas, A. Jr, 2012. *Spectral Methods in Fluid Dynamics*, Springer Science & Business Media.

Carpenter, M.H. & Kennedy, C.A., 1994. Fourth-order 2n-storage Runge-Kutta schemes, Nasa Report, TM 109112.

Chabaud, B. & Cockburn, B., 2012. Uniform-in-time superconvergence of hdg methods for the heat equation, *Math. Comput.*, **81**(277), 107–129.

Chaljub, E. & Valette, B., 2004. Spectral element modelling of three-dimensional wave propagation in a self-gravitating earth with an arbitrarily stratified outer core, *Geophys. J. Int.*, **158**(1), 131–141.

Chaljub, E., Capdeville, Y. & Vilotte, J.-P., 2003. Solving elastodynamics in a fluid-solid heterogeneous sphere: a parallel spectral ele-

ment approximation on non-conforming grids, *J. Comp. Phys.*, **187**(2), 457–491.

Chaljub, E., Komatitsch, D., Vilotte, J.-P., Capdeville, Y., Valette, B. & Festa, G., 2007. Spectral Element analysis in seismology, in *Advances in Wave Propagation in Heterogeneous Media*, Vol. 72 of Advances in Geophysics series, pp. 105–108, eds Wu, R.-S., Maupin, V. & Dmowvska, R., Elsevier.

Chapman, C., 2004. *Fundamentals of Seismic Wave Propagation*, Cambridge Univ. Press.

Chung, E. & Engquist, B., 2006. Optimal discontinuous Galerkin methods for wave propagation, *SIAM J. Numer. Anal.*, **44**(5), 2131–2158.

Cockburn, B. & Gopalakrishnan, J., 2004. A characterization of hybridized mixed methods for second order elliptic problems, *SIAM J. Numer. Anal.*, **42**(1), 283–301.

Cockburn, B. & Quenneville-Bélair, V., 2014. Uniform-in-time superconvergence of the HDG methods for the acoustic wave equation, *Math. Comput.*, **83**(285), 65–85.

Cockburn, B. & Shi, K., 2013. Superconvergent HDG methods for linear elasticity with weakly symmetric stresses, *IMA J. Numer. Anal.*, **33**(3), 747–770.

Cockburn, B. & Stolarski, H.K., 2014. Analysis of an HDG method for linear elasticity, *Int. J. Numer. Methods Eng.*, **102**, 551–575.

Cockburn, B., Gopalakrishnan, J. & Wang, H., 2007. Locally conservative fluxes for the continuous Galerkin method, *SIAM J. Numer. Anal.*, **45**(4), 1742–1776.

Cockburn, B., Gopalakrishnan, J. & Lazarov, R., 2009. Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems, *SIAM J. Numer. Anal.*, **47**(2), 1319–1365.

Cockburn, B., Nguyen, N.C. & Peraire, J., 2010. A Comparison of hdg Methods for Stokes flow, *J. Sci. Comput.*, **45**(1), 215–237.

Cockburn, B., Gopalakrishnan, J., Nguyen, N., Peraire, J. & Sayas, F.-J., 2011. Analysis of HDG methods for Stokes flow, *Math. Comput.*, **80**(274), 723–760.

Cohen, G., 2002. *Higher-order Numerical Methods for Transient Wave Equations*, Scientific Computation, Springer-Verlag.

Cohen, G., Joly, P., Roberts, J. & Tordjman, N., 2001. Higher order triangular finite elements with mass lumping for the wave equation, *SIAM J. Numer. Anal.*, **38**(6), 2047–2078.

Cupillard, P., Delavaud, E., Burgos, G., Festa, G., Vilotte, J.-P., Capdeville, Y. & Montagner, J.-P., 2012. RegSEM: a versatile code based on the spectral element method to compute seismic wave propagation at the regional scale, *Geophys. J. Int.*, **188**(3), 1203–1220.
de Basabe, J., Sen, J. & Wheeler, M., 2008. The interior penalty discontinuous Galerkin method for elastic wave propagation: grid dispersion, *Geophys. J. Int.*, **175**(1), 83–93.

De Hoop, A., 1960. A modification of Cagniard's method for solving seismic pulse problems, *Appl. Sci. Res. B*, **8**(1), 349–356.

de la Puente, J., Käser, M., Dumbser, M. & Igel, H., 2007. An arbitrary high order discontinuous Galerkin method for elastic waves on unstructured meshes – IV. anisotropy, *Geophys. J. Int.*, **169**(3), 1210–1228.

Delcourte, S. & Glinsky, N., 2015. Analysis of a high-order space and time discontinuous Galerkin method for elastodynamic equations. application to 3D wave propagation, *ESAIM: Math. Modelling Numer. Anal.*, **49**(4), 1085–1126.

Delcourte, S., Fezoui, L. & Glinsky-Olivier, N., 2009. A high-order discontinuous Galerkin method for the seismic wave propagation, *ESAIM Proc.*, **27**, 70–89.

De Martin, F., 2011. Verification of a spectral-element method code for the Southern California EarthquakeCenter LOH.3 viscoelastic case, *Bull. seism. Soc. Am.*, **101**(6), 2855–2865.

Deville, M.O., Fischer, P.F. & Mund, E.H., 2002. *High-order Methods for Incompressible Fluid Flow*, Vol. 9, Cambridge Univ. Press.

Diaz, J., 2005. Approches analytiques et numériques de problèmes de transmission en propagation d'ondes en régime transitoire. Application au couplage fluide-structure et aux méthodes de couches parfaitement adaptées, *PhD thesis*, ENSTA ParisTech.

Diaz, J. & Joly, P., 2005. Robust high order non-conforming finite element formulation for time domain fluid-structure interaction, *J. Comput. Acoust.*, **13**(03), 403–431.

Dumbser, M., Käser, M. & Toro, E., 2007. An arbitrary high order discontinuous Galerkin method for elastic waves on unstructured meshes – V. Local time stepping and p-adaptivity, *Geophys. J. Int.*, **171**(2), 695–717.

Etienne, V., Chaljub, E., Virieux, J. & Glinsky, N., 2010. An hp-adaptive discontinuous Galerkin finite-element method for 3D elastic wave modelling, *Geophys. J. Int.*, **183**(2), 941–962.

Faccioli, E., Maggio, F., Paolucci, R. & Quarteroni, A., 1997. 2D and 3D elastic wave propagation by a pseudo-spectral domain decomposition method, *J. Seismol.*, **1**(3), 237–251.

Festa, G. & Vilotte, J.-P., 2005. The Newmark scheme as a velocity-stress time-staggering: an efficient PML implementation for spectral element simulations of elastodynamics, *Geophys. J. Int.*, **161**(3), 789–812.

Fraejis de Veubeque, B., 1965. Displacement and equilibrium models in the finite element method, in *Stress Analysis*, eds Zienkiewicz, O. & Holister, G., Wiley.

Fu, G., Cockburn, B. & Stolarski, H., 2015. Analysis of an HDG method for linear elasticity, *Int. J. Numer. Methods Eng.*, **102**(3-4), 551–575.

Gedney, S.D. & Zhao, B., 2010. An auxiliary differential equation formulation for the complex-frequency shifted PML, *IEEE Trans. Antennas Propag.*, **58**(3), 838–847.

Giorgiani, G., Fernández-Méndez, S. & Huerta, A., 2013. Hybridizable discontinuous Galerkin p-adaptivity for wave propagation problems, *Int. J. Numer. Methods Fluids*, **72**(12), 1244–1262.

Green, D.N., Guilbert, J., Le Pichon, A., Sebe, O. & Bowers, D., 2009. Modelling ground-to-air coupling for the shallow ML 4.3 Folkestone, United Kingdom, earthquake of 28 April 2007, *Bull. seism. Soc. Am.*, **99**(4), 2541–2551.

Grote, M., Schneebeli, A. & Schötzau, D., 2006. Discontinuous Galerkin finite element method for the wave equation, *SIAM J. Numer. Anal.*, **44**(6), 2408–2431.

Hairer, E., 2011. *Solving Differential Equations on Manifolds*, Lecture Notes, Université de Geneve.

Hairer, E., Nrsett, S.P. & Wanner, G., 2010. *Solving Ordinary Differential Equations: Nonstiff Problems. v. 2: Stiff and Differential-algebraic Problems*, Springer Verlag.

Hermann, V., Käser, M. & Castro, C., 2011. Non-conforming hybrid meshes for efficient 2-D wave propagation using the discontinuous Galerkin method, *Geophys. J. Int.*, **184**(2), 746–758.

Hesthaven, J. & Warburton, T., 2008. *Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications*, Springer.

Hungria, A., Prada, D. & Sayas, F.-J., 2017. Hdg methods for elastodynamics, *Comput. Math. Appl.*, **74**(11), 2671–2690.

Käser, M. & Dumbser, M., 2006. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes - i. the two-dimensional isotropic case with external source terms, *Geophys. J. Int.*, **166**(2), 855–877.

Komatitsch, D. & Vilotte, J.-P., 1998. The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures, *Bull. seism. Soc. Am.*, **88**(2), 368–392.

Komatitsch, D., Vilotte, J.-P., Vai, R., Castillo-Covarrubias, J.M. & Sánchez-Sesma, F.J., 1999. The spectral element method for elastic wave equations – application to 2-D and 3-D seismic problems, *Int. J. Numer. Methods Eng.*, **45**(9), 1139–1164.

Komatitsch, D., Barnes, C. & Tromp, J., 2000. Wave propagation near a fluid-solid interface: A spectral-element approach, *Geophysics*, **65**(2), 623–631.

Komatitsch, D., Martin, R., Tromp, J., Taylor, M. & Wingate, B., 2001. Wave propagation in 2-D elastic media using a spectral element method with triangles and quadrangles, *J. Comput. Acoust.*, **9**(2), 703–718.

Kopriva, D.A., 2006. Metric identities and the discontinuous spectral element method on curvilinear meshes, *J. Sci. Comput.*, **26**(3), 301–327.

Kopriva, D.A., 2009. *Implementing Spectral Methods for Partial Differential Equations: Algorithms for Scientists and Engineers*, Springer Science & Business Media.

Lamb, H., 1904. On the propagation of tremors over the surface of an elastic solid, *Phil. Trans. R. Soc. A*, **203**, 1–42.

LeVeque, R.J., 2002. *Finite Volume Methods for Hyperbolic Problems*, Cambridge Univ. Press.

Liu, T., Sen, M.K., Hu, T., De Basabe, J.D. & Li, L., 2012. Dispersion analysis of the spectral element method using a triangular mesh, *Wave Motion*, **49**(4), 474–483.

Maday, Y. & Patera, A.T., 1989. Spectral element methods for the incompressible Navier-Stokes equations, in *State of the Art Survey in Computational Mechanics*, pp. 71–143, eds Noor, A.K. & Oden, J.T., ASME.

Martin, R., Komatitsch, D., Gedney, S.D. & Bruthiaux, E., 2010. A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using Auxiliary Differential Equations (ADE-PML), *Comput. Modelling Eng. Sci. (CMES)*, **56**(1), 17.

Maupin, V. & Wu, R.-S., 2007. Preface, in *Advances in Geophysics: Advances in Wave Propagation in Heterogenous Earth*, Vol. 48, pp. v–viii, eds Wu, R.-S., Maupin, V. & Dmowska, R., Elsevier.

Mazzieri, I. & Rapetti, F., 2012. Dispersion analysis of triangle-based spectral element methods for elastic wave propagation, *Numer. Algorithms*, **60**(4), 631–650.

Mercerat, E., Vilotte, J.-P. & Sánchez-Sesma, F., 2006. Triangular spectral element simulation of 2D elastic wave propagation using unstructured triangular grids, *Geophys. J. Int.*, **166**(2), 679–698.

Mercerat, E.D. & Glinsky, N., 2015. A nodal high-order discontinuous Galerkin method for elastic wave propagation in arbitrary heterogeneous media, *Geophys. J. Int.*, **201**(2), 1101–1118.

Meyers, R.J., Tautges, T.J. & Tuchinsky, P.M., 1998. The “hex-tet” hex-dominant meshing algorithm as implemented in cubit., in *IMR*, pp. 151–158.

Minisini, S., Zebel, E., Konokov, A. & Mulder, W., 2013. Local time stepping with the discontinuous Galerkin method for wave propagation in 3D heterogeneous media, *Geophysics*, **78**(3), 67–77.

Moczo, P., Kristek, J. & Gális, M., 2014. *The Finite-Difference Modelling of Earthquake Motions: Waves and Ruptures*, Cambridge Univ. Press.

Nguyen, N., Peraire, J. & Cockburn, B., 2011a. Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell's equations, *J. Comput. Phys.*, **230**(19), 7151–7175.

Nguyen, N.C., Peraire, J. & Cockburn, B., 2011b. High-order implicit hybridizable discontinuous Galerkin methods for acoustics and elastodynamics, *J. Comput. Phys.*, **230**(10), 3695–3718.

Patera, A.T., 1984. A spectral element method for fluid dynamics: laminar flow in a channel expansion, *J. Comput. Phys.*, **54**(3), 468–488.
Peraire, J., Nguyen, N. & Cockburn, B., 2010. A hybridizable discontinuous Galerkin method for the compressible Euler and Navier-Stokes equations, in 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, p. 363.

Perugia, I. & Schötzau, D., 2001. On the coupling of local discontinuous Galerkin and conforming finite element methods, *J. Sci. Comput.*, **16**(4), 411–433.

Peter, D. et al., 2011. Forward and adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes, *Geophys. J. Int.*, **186**(2), 721–739.

Raviart, P. & Thomas, J., 1977. A mixed finite element method for second order elliptic problems, in *Mathematical Aspects of Finite Element Method, Lecture Notes in Mathematics, Vol. 606*, pp. 292–315, eds Galleni, I. & Magenes, E., Springer Verlag.

Seriani, G. & Priolo, E., 1994. Spectral element method for acoustic wave simulation in heterogeneous media, *Finite Elem. Anal. Des.*, **16**(34), 337–348.

Soon, S.-C., Cockburn, B. & Stolarski, H.K., 2009. A hybridizable discontinuous Galerkin method for linear elasticity, *Int. J. Numer. Methods Eng.*, **80**(8), 1058–1092.

Stanglmeier, M., Nguyen, N., Peraire, J. & Cockburn, B., 2016. An explicit hybridizable discontinuous Galerkin method for the acoustic wave equation, *Comput. Methods Appl. Mech. Eng.*, **300**, 748–769.

Tsuboi, S., Ando, K., Takayuki, M., Peter, D., Komatitsch, D. & Tromp, J., 2016. A 1.8 trillion degree of freedom, 1.24 Petaflops global seismic wave simulation on the K computer, *Int. J. High Perform. Comput. Appl.*, **30**(4), 411–422.

Virieux, J., Etienne, V. & Cruz-Atienza, V., 2012. Modelling seismic wave propagation for geophysical imaging, in *Seismic Waves – Research and Analysis*, Chap. 13, pp. 253–304, ed. Kanao, M., InTech.

Walker, K.T., Pichon, A.L., Kim, T.S., Groot-Hedlin, C., Che, I.-Y. & Garcés, M., 2013. An analysis of ground shaking and transmission loss from in-frasound generated by the 2011 Tohoku earthquake, *J. geophys. Res.*, **118**(23), 12–831.

Wenk, S., Pelties, C., Igel, H. & Käser, M., 2013. Regional wave propagation using the discontinuous Galerkin method, *Solid Earth*, **4**, 43–57.

Wilcox, L.C., Stadler, G., Burstedde, C. & Ghattas, O., 2010. A high-order discontinuous Galerkin method for wave propagation through coupled elastic-acoustic media, *J. Comput. Phys.*, **229**(24), 9373–9396.

Yamakawa, S. & Shimada, K., 2003. Fully-automated hex-dominant mesh generation with directionality control via packing rectangular solid cells, *Int. J. Numer. Methods Eng.*, **57**(15), 2099–2129.

Ye, R., de Hoop, M.V., Petrovitch, C.L., Pyrak-Nolte, L.J. & Wilcox, L.C., 2016. A discontinuous Galerkin method with a modified penalty flux for the propagation and scattering of acousto-elastic waves, *Geophys. J. Int.*, **205**(2), 1267–1289.

## APPENDIX A: LOCAL POST-PROCESSING

The post-processing is an elementwise operation that provides a more accurate local solution. Due to its local nature, the post-processing adds only a marginal cost to the whole simulation.

When polynomials of degree $N$ are used, the post-processing accuracy relies on the fact that when both the primary field, i.e the velocity (resp. pressure), and the gradient field, that is, the strain (resp. velocity), in the elastic (resp. acoustic) domain converge with order $N + 1$, then the primary field can be post-processed to obtain a superconvergence of order $N + 2$.

For the acoustic equations, Nguyen et al. (2011b) proposed a post-processing of the pressure that numerically achieves a convergence with order $N + 2$, a result that has been theoretically proven by Cockburn & Quenneville-Bélair (2014). In this study, the same post-processing is considered. At the end of each time step, a new post-processed pressure $p_h^* \in \mathcal{P}_{N+1}(K)$, is computed in two stages. The first stage is to get an approximation of the gradient of the pressure $q_h \in (\mathcal{P}_N(K))^{nd}$ solving

$$ \int_K q_h \cdot w \, d\Omega = \int_K (\kappa \epsilon_h) \nabla \cdot w \, d\Omega + \int_{\partial K} \hat{p}_h w \cdot n \, d\Gamma \quad \forall w \in (\mathcal{P}_N(K))^{nd}. \quad (A1) $$

Then, $p_h^*$ is computed from $q_h$ solving the following system

$$ \int_K \nabla p_h^* \cdot \nabla w \, d\Omega = \int_K q_h \cdot \nabla w \, d\Omega \quad \forall w \in \mathcal{P}_{N+1}(K), \quad (A2a) $$

$$ \int_K p_h^* \, d\Omega = - \int_K \kappa \epsilon_h \, d\Omega. \quad (A2b) $$

The first step is straightforward since it only involves the inversion of a diagonal matrix of size $(n_d(N+1)) \times (n_d(N+1))$. The second step is more expensive since it involves the resolution of a rectangular system of size $(N+3) \times (N+2)$.

In the elastic domain, the superconvergence property may not be obtained, if the symmetric strain tensor does not converge with the optimal order $N + 1$, which in turn depends on the space used to build the approximated strain, and how its symmetry is enforced. Here is a brief review of the post-processing results obtained with different HDG approaches in elastostatics and elastodynamics.

It has been shown by Cockburn & Shi (2013) that for a weakly symmetric stresses HDG formulations, the superconvergence can be achieved. For the displacement gradient-velocity-pressure HDG formulation used by Nguyen et al. (2011b), the superconvergence has been both observed numerically and proven theoretically by Cockburn et al. (2011).

The HDG method proposed by Soon et al. (2009) is the elastostatic equivalent to our formulation when C is constant. This formulation has numerically achieved an optimal superconvergence in some circumstances, although Fu et al. (2015) demonstrated that it was a matter of chance. Indeed the a priori error analysis ensures orders of convergence $N + 1$ for the displacement, and only $N + 1/2$ for the symmetric strain tensor. Thus no superconvergence property is ensured with this formulation, although it may be sometimes observed. Our numerical experiments agree with these conclusions. Indeed our results suggest that the post-processed velocity is always more precise, and converges with a higher order when compared with the regular velocity approximation, although the new order of convergence is less than $N + 2$ most of the time.
The proposed two-stage post-processing is inspired from Soon et al. (2009). First, the strain rate $d_h \in (P_N(K))_{\text{sym}}^{n_d^2}$ is computed solving

$$ \int_K d_h : r \, d\Omega = - \int_K v_h \cdot \nabla r \, d\Omega + \int_{\partial K} \hat{v}_h \cdot r \cdot n \, d\Gamma \quad \forall r \in (P_N(K))_{\text{sym}}^{n_d^2}. \quad (A3) $$

Then, the post-processed velocity $v_h^* \in (P_{N+1}(K))^{n_d}$ is computed solving the system

$$ \int_K \nabla^s v_h^* : \nabla^s w \, d\Omega = \int_K d_h : \nabla^s w \, d\Omega, \quad \forall w \in (P_{N+1}(K))^{n_d}, \quad (A4a) $$

$$ \int_K v_h^* \, d\Omega = \int_K v_h \, d\Omega, \quad (A4b) $$

$$ \int_K \mathrm{rot} v_h^* \, d\Omega = \int_K \mathrm{rot} v_h \, d\Omega. \quad (A4c) $$

The equation (A4a) determines completely the symmetric gradient of the post-processed velocities. To obtain the velocity field $v_h^*$ from it, the rigid motion needs to be prescribed via (A4b) for its translation components and via (A4c) for its rotation components.

The first stage involves the inversion of a diagonal matrix whose size is the dimension of $(P_N(K))_{\text{sym}}^{n_d^2}$. The second stage involves the resolution of a rectangular system whose size is $(2(N+2)+3) \times (2(N+2))$ in 2-D, and $(3(N+2)+6) \times (3(N+2))$ in 3-D.

## APPENDIX B: NOTES FOR 2-D IMPLEMENTATION

Here is some information regarding the 2-D implementation of the hybrid Riemann solver. In 2-D, the elastodynamic variables become $q = (\epsilon_{11}, \epsilon_{22}, \epsilon_{13}, \rho v_1, \rho v_2)^T$, and the acoustic ones become $q = (\epsilon, \rho v_1, \rho v_2)^T$. The elastic flux matrix (10) and the acoustic flux matrix (16) become respectively

$$ \bar{A}_n = - \begin{pmatrix} 0 & 0 & 0 & \frac{n_x}{\rho} & 0 \\ 0 & 0 & 0 & 0 & \frac{n_y}{\rho} \\ 0 & 0 & 0 & \frac{n_y}{2\rho} & \frac{n_x}{2\rho} \\ (\lambda + 2\mu)n_x & \lambda n_x & 2\mu n_y & 0 & 0 \\ \lambda n_y & (\lambda + 2\mu)n_y & 2\mu n_x & 0 & 0 \end{pmatrix} \quad \text{and} \quad \bar{A}_n = - \begin{pmatrix} 0 & \frac{n_x}{\rho} & \frac{n_y}{\rho} \\ \kappa n_x & 0 & 0 \\ \kappa n_y & 0 & 0 \end{pmatrix}. \quad (B1) $$

The elastic $\tau$ matrix (37) and the acoustic-elastic $\tau_{ac}^{-}$ matrix (52) become respectively

$$ \tau^{-} = \begin{pmatrix} \rho(c_p n_x^2 + c_s n_y^2) & \rho(c_p - c_s)n_x n_y \\ \rho(c_p - c_s)n_x n_y & \rho(c_s n_x^2 + c_p n_y^2) \end{pmatrix}^{-} \quad \text{and} \quad \tau_{ac}^{-} = \begin{pmatrix} \rho^{-} c_p^{-} n_x^2 & \rho^{-} c_p^{-} n_x n_y \\ \rho^{-} c_p^{-} n_x n_y & \rho^{-} c_p^{-} n_y^2 \end{pmatrix}. \quad (B2) $$