DG method for acousto-elastic waves ] A Discontinuous Galerkin method with a modified penalty flux for the propagation and scattering of acousto-elastic waves R. Ye et al.] Ruichao Ye, Maarten V. de Hoop, Christopher L. Petrovitch, Laura J. Pyrak-Nolte, Lucas C. Wilcox
Department of Earth Science, Rice University, Houston TX, USA
Simons Chair in Computational and Applied Mathematics and Earth Science, Rice University, Houston TX, USA
Applied Research Associates, Inc., Raleigh-Durham NC, USA
Department of Physics, Purdue University, West Lafayette IN, USA
Department of Applied Mathematics, Naval Postgraduate School, Monterey CA, USA


We develop an approach for simulating acousto-elastic wave phenomena, including scattering from fluid-solid boundaries, where the solid is allowed to be anisotropic, with the Discontinuous Galerkin method. We use a coupled first-order elastic strain-velocity, acoustic velocity-pressure formulation, and append penalty terms based on interior boundary continuity conditions to the numerical (central) flux so that the consistency condition holds for the discretized Discontinuous Galerkin weak formulation. We incorporate the fluid-solid boundaries through these penalty terms and obtain a stable algorithm. Our approach avoids the diagonalization into polarized wave constituents such as in the approach based on solving elementwise Riemann problems.

iscontinuous Galerkin method – penalty flux – fluid-solid boundaries – anisotropy

1 Introduction

The accurate computation of waves in realistic three-dimensional Earth models represents an ongoing challenge in local, regional, and global seismology. Here, we focus on simulating coupled acousto-elastic wave phenomena including scattering from fluid-solid boundaries, where the solid is allowed to be anisotropic, with the Discontinuous Galerkin method. Of particular interest are applications in geophysics, namely, marine seismic exploration and global Earth inverse problems using earthquake-generated seismic waves as the probing field. In the first application, we are concerned with the presence of the ocean bottom and in the second one with the core-mantle-boundary (CMB) and inner-core-boundary (ICB). Our formulation closely follows the analysis of existence of (weak) solutions of hyperbolic first-order systems of equations by [Blazek et al.(2013)Blazek, Stolk, & Symes]. We use an unstructured tetrahedral mesh with local refinement to accommodate highly heterogeneous media and complex geometries, which is also an underlying motivation for employing the Discontinuous Galerkin method from a computational point of view.

In the past three decades, a wide variety of numerical techniques has been employed in the development of computational methods for simulating seismic waves. The most widely used one is based on the finite difference method [e.g., [Madariaga(1976)] and [Virieux(1986)]]. This method has been applied to computing the wavefield in three-dimensional local and regional models [e.g., [Graves(1996)] and [Ohminato & Chouet(1997)]]. The use of optimal or compact finite-difference operators has provided a certain improvement [e.g., [Zingg et al.(1996)Zingg, Lomax, & Jurgens] and [Zingg(2000)]]. Methods that resort to spectral and pseudospectral techniques based on global gridding of the model have also been used both in regional [e.g., [Carcione(1994)]] and global [e.g., [Tessmer et al.(1992)Tessmer, Kessler, Kosloff, & Behle] and [Igel(1999)]] seismic wave propagation and scattering problems. However, because of the use of global basis functions (polynomial: Chebyshev or Legendre, or harmonic: Fourier), these techniques are limited to coefficients which are (piecewise) sufficiently smooth. The finite difference method suffers from a limited accuracy in the presence of a free surface or surface discontinuities with topography within the model [e.g., [Robertsson(1996)] and [Symes & Vdovina(2009)]]. A procedure for the stable imposition of free-surface boundary conditions for a second-order formulation can be found in [Appelö & Petersson(2009)]. Another approach, belonging to a broader family of interface methods, handles both free surfaces [e.g., [Lombard et al.(2008)Lombard, Piraux, Gélis, & Virieux]] and fluid-solid interfaces [e.g., [Lombard & Piraux(2004)]] in such a way, conjectured by the authors, that enables higher-order accuracy to be obtained. [Kozdon et al.(2013)Kozdon, Dunham, & Nordstr枚m] use summation-by-parts finite difference operators along with a weak enforcement of boundary conditions to develop a multi-block finite difference scheme which achieves higher-order accuracy for complex geometries.

A key development in the computation of seismic waves has been based on the spectral element method (SEM) [[Komatitsch & Tromp(2002)]]. In its original formulation, in terms of displacement [[Komatitsch & Vilotte(1998)]], continuity of displacement and velocity is enforced everywhere within the model. In the case of a boundary between an inviscid fluid and a solid, however, the kinematic boundary condition is perfect slip; therefore, only the normal component of velocity is continuous across such a boundary, and thus this formulation is not applicable. Some classical finite-element methods (FEMs) alternatively introduce coupling conditions on fluid-solid interfaces between displacement in the solid and pressure in the fluid [e.g. [Zienkiewicz & Bettess(1978), Bermúdez et al.(1999)Bermúdez, Hervella-Nieto, & Rodrıguez]].

The FEM and SEM are commonly (but not exclusively) based on the second-order form of the system of equations describing acousto-elastic waves. In this case, the acousto-elastic interaction is effected by coupling the respective wave equations through appropriate interface conditions. To resolve the coupling, a predictor-multicorrector iteration at each time step has been used [[Komatitsch et al.(2000)Komatitsch, Barnes, & Tromp], [Chaljub et al.(2003)Chaljub, Capdeville, & Vilotte]]. A computationally more efficient time stepping method for global seismic wave propagation accommodating the effects of fluid-solid boundaries, as well as transverse isotropy with a radial symmetry axis and radial models of attenuation, was proposed in [Komatitsch et al.(2005)Komatitsch, Tsuboi, & Tromp]. It uses a velocity potential formulation and a second-order accurate Newmark time integration, in which a time step is first performed in the acoustic fluid and then in the elastic solid using interface values based on the fluid solution. Currently the SEM is used in a variety of implementations in global and regional seismic simulation, with the effects of variations in elastic parameters, density, ellipticity, topography and bathymetry, fluid-solid interfaces, anisotropy, and self-gravitation included [e.g. [Carrington et al.(2008)Carrington, Komatitsch, Laurenzano, Tikir, Michéa, Le Goff, Snavely, & Tromp]].

In contrast to classical finite element discretizations, the Discontinous Galerkin (DG) method imposes continuity of approximate solutions between elements only weakly through a numerical flux. The Discontinuous Galerkin method has been employed for solving second-order wave equations in both the acoustic and elastodynamic settings [e.g. [Rivière & Wheeler(2000)], [Grote et al.(2006)Grote, Schneebeli, & Schötzau], [Chung & Engquist(2006)] and [De Basabe et al.(2008)De Basabe, Sen, & Wheeler]]. [Etienne et al.(2010)Etienne, Chaljub, Virieux, & Glinsky] employ a central numerical flux in a DG scheme combined with a leap-frog time integration for the velocity-stress elastic-wave formulation. [Dumbser & Käser(2006), Käser & Dumbser(2008)] developed a non-conservative formulation with an upwind numerical flux using only the material properties from the side of the interface that is opposite to the outer normal direction. [Wilcox et al.(2010)Wilcox, Stadler, Burstedde, & Ghattas] derived an upwind numerical flux by solving the exact Riemann problem on interior boundaries of each element with material discontinuities based on a velocity-strain formulation of the coupled acousto-elastic equations.

In this work, we essentially extend the upwind flux, given by [Warburton(2013)] for hyperbolic systems, to a penalty flux based on the boundary continuity condition for general fluid-solid interfaces. The novelties of our approach are the following: we

  1. use a coupled first-order elastic strain-velocity, acoustic velocity-pressure formulation,

  2. obtain a self-consistent Discontinuous Galerkin weak formulation without diagonalization into polarized wave constituents,

  3. append penalty terms, derived from interior boundary continuity conditions, with an appropriate weight to the numerical (central) flux so that the consistency condition holds for the discretized Discontinuous Galerkin weak formulation,

  4. incorporate fluid-solid boundaries through the mentioned penalty terms.

We note that the DG method is naturally adapted to well-posedness, in the sense that it makes use of coercivity of the operator defining the part of the system containing the spatial derivatives separately in the solid and fluid regions.

2 The system of equations describing acousto-elastic waves

We consider a bounded domain which is divided into solid and fluid regions, and , respectively. The interior boundaries include solid-solid interface , fluid-fluid interface , and fluid-solid interface , (where we distinguish whether the fluid or solid is on a particular side). We present the weak form of the coupled acousto-elastic system of equations.

Hooke’s law in an elastodynamical system is expressed by relating stress, , and strain, . Assuming small deformations gives a linear relationship, that is, , where is the stiffness tensor. Through the relevant symmetries, this tensor only contains 21 independent components. We use the Voigt notation which simplifies the writing of tensors while introducing and . In this notation the stiffness tensor takes the form of a 6 by 6 matrix, C, defined by,


The isotropic case is obtained by setting all of the components to zero except for , , , , and ; are the Lamé parameters. Furthermore, denotes the density. The anisotropic elastodynamical equations are written in terms of the strain, E, and the particle velocity, v,


in . In fluid regions, , we use the pressure-velocity formulation,


Here, is the pressure, while we use to distinguish acoustic field quantities and material parameters from the elastic ones. In the above, denotes a volume source density of injection and f denotes a volume source density of force.

The solid-solid, fluid-solid and fluid-fluid boundary conditions are given by


The convention is determined by the direction of the interface normal, n. The outer normal vector points in the direction of the “” side of the interface.

We introduce test functions (tensors) in the solid regions and in the fluid regions, which are assumed to be contained in the same spaces and satisfy the same boundary conditions as and . Using (2) and (3), we find that


Assuming an outer traction-free boundary condition in (5b) and an outer pressure-free boundary condition in (5c), and applying an integration by parts, we obtain


We use the fluid-solid boundary conditions (4b), replacing the fluid-solid surface integrals in (6a) and (6b) by taking the average of both sides consistent with a central flux scheme, and obtain


This form of the equations is analogous to the one used in the spectral element method, see [Chaljub & Valette(2004)]. Applying an integration by parts, again, in (7), we recover the coupled strong formulation,


We use this system of equations together with (5a) and (5d) to develop our Discontinuous Galerkin method based approach.

3 Discontinuous Galerkin method with fluid-solid boundaries

The domain is partitioned into elements, . We distinguish elements, , in the solid regions from elements, , in the fluid regions. Correspondingly, we distinguish fluid-fluid (), solid-solid () and fluid-solid () faces for each element; thus the interior boundaries are decomposed as

and so are the elements’ boundaries: and . The mesh size, , is defined as the maximum radius of each tetrahedral’s inscribed sphere.

We introduce the broken polynomial space where the local space is defined elementwise as , with a set of polynomial basis further discussed in Section 3.2. The subscript “” indicates the refinement of with decrease in mesh size. The semi-discrete time-domain, discontinuous Galerkin formulation using a central flux yields: Find , with each component for each one of them in such that




hold for each element or , for all test functions . The notations and indicate polynomial approximation of f and . Here,


in the solid regions, while


in the fluid regions, using interior boundary continuity conditions. A similar formulation for Maxwell’s equations, using the central flux, can be found in [Hesthaven & Warburton(2007), Chapter 10, Page 434].

3.1 Energy function of central flux

We consider a time-dependent energy function comprising both the solid and fluid regions, , with


The functions in (13) define a norm both in the solid and in the fluid regions. Taking the time derivative and noting that C is symmetric, we have


Starting from (9) and (10) and carrying out the summation over all the elements yields


This property is obtained as follows:

In (9) and (10) we let , and obtain elementwise


and similarily


From (9), (11), (14) and (17),


In the above,


The surface integration terms cancel out when summed from both sides of the solid-solid interfaces because of the continuity condition (4a) and the opposite outer normal directions. We are left with the contributions from solid-fluid inner faces, ,


A similar result in the fluid region obtained from (10), (12), (15) and (18) yields


and the surface integration terms on the solid-fluid and fluid-solid interfaces in (20) and (21) cancel out due to (4b). Therefore (16) is obtained. We note that the surface integration along solid-fluid interfaces and are essential to guarantee energy conservation.

3.2 Nodal basis functions

The discretized solution follows an expansion, componentwise, into nodal trial basis functions of order , as is in [Hesthaven & Warburton(2007)],


and similarly for the other fields, . The superscript, , indicates a local expansion within element . In the above, is a set of three-dimensional Lagrange polynomials associated with the nodal points, (see Figure 1), with each polynomial defined as

We use the warp & blend method [[Warburton(2006)]] to determine the coordinates of nodal points in the tetrahedron by numerically minimizing the Lebesgue constant of interpolation. For an order interpolation there are nodal points.

= 1 = 3 = 8
Figure 1: Warp & blend tetrahedral nodal point distribution for = 1, 3, 8. For clarity only facial nodes are illustrated.

The medium coefficients are expanded in a likewise manner


and similarly for . When refining a mesh, we expect an increase in number of elements with decreased size.

3.3 The system of equations in matrix form

To simplify the notation in the further development of a numerical scheme, we introduce a joint matrix form of the system of equations. We map the components of and to and matrices, respectively,


and, correspondingly, the components of body forces f and to the matrix

Equations (2) and (3) attain the form




that is,