A weight-adjusted discontinuous Galerkin method for the poroelastic wave equation: penalty fluxes and micro-heterogeneities

TL;DR

This paper presents a high-order weight-adjusted discontinuous Galerkin method for 3D poroelastic wave simulations, emphasizing energy stability, micro-heterogeneity handling, and avoiding complex Jacobian diagonalization.

Contribution

The paper introduces a novel energy-stable penalty flux and a weight-adjusted mass matrix approach for micro-heterogeneity in a high-order DG scheme for poroelastic waves.

Findings

Proven convergence of the numerical scheme.

High-order accuracy maintained near material discontinuities.

Comparable or improved performance against spectral element methods.

Abstract

We introduce a high-order weight-adjusted discontinuous Galerkin (WADG) scheme for the numerical solution of three-dimensional (3D) wave propagation problems in anisotropic porous media. We use a coupled first-order symmetric stress-velocity formulation. Careful attention is directed at (a) the derivation of an energy-stable penalty-based numerical flux, which offers high-order accuracy in presence of material discontinuities, and (b) proper treatment of micro-heterogeneities (sub-element variations) in the numerical scheme. The use of a penalty-based numerical flux avoids the diagonalization of Jacobian matrices into polarized wave constituents necessary when solving element-wise Riemann problems. Micro-heterogeneities are accurately and stably incorporated in the numerical scheme using easily-invertible weight-adjusted mass matrices. The convergence of the proposed numerical scheme is proven and verified by using convergence studies against analytical plane wave solutions. The proposed method is also compared against an existing implementation using the spectral element method to solve the poroelastic wave equation.

Tables

Table 1. Table 1: Material properties for several poroelastic media used in the examples [ 9 ]

Properties	Sandstone	Epoxy-glass	Sandstone	Shale
	(Orthotropic)	(Orthotropic)	(Isotropic)	(Isotropic)
$K_s$(GPa)	80	40	40	7.6
${\rho }_s$ $($kg/m3$)$	2500	1815	2500	2210
$c_{11}$(GPa)	71.8	39.4	36	11.9
$c_{12}$(GPa)	3.2	1.2	12	3.96
$c_{13}$(GPa)	1.2	1.2	12	3.96
$c_{33}$(GPa)	53.4	13.1	36	11.9
$c_{55}$(GPa)	26.1	3	12	3.96
$\varphi$	0.2	0.2	0.2	0.16
${\kappa }_1\left({10}^{-15}{m}^2\right)$	600	600	600	100
${\kappa }_3\left({10}^{-15}{m}^2\right)$	100	100	600	100
$T_1$	2	2	2	2
$T_3$	3.6	3.6	2	2
$K_f$(GPa)	2.5	2.5	2.5	2.5
${\rho }_f$(Kg/m3)	1040	1040	1040	1040
$\eta$$({10}^{-3}$ Kg/m.s)	1	1	1	1

Table 2. Table 2: Material properties used for comparison of solution with SPECFEM implementations [ 4 ]

Properties	Medium I	Medium II	Medium III
${\rho }_s$ $($kg/m3$)$	2650	2200	2650
${\rho }_f$ $($kg/m3$)$	880	950	750
$\varphi$	0.1	0.4	0.2
$T_1$	2	2	2
$T_3$	2	2	2
$K_s$(GPa)	12.2	6.9	6.9
$K_f$(GPa)	1.985	2.0	2.0
${\kappa }_{fr}$(GPa)∗	9.6	6.7	6.7
${\mu }_{fr}^{\ast \ast }$(GPa)∗∗	5.1	3.0	3.0
$\eta$$({10}^{-3}$ Kg/m.s)	0	0	0
∗: Frame bulk modulus, ∗∗: Frame shear modulus.

Equations

$${\tau }_{ij}$$

$$p$$

$$\rho \frac{\partial \mathbf{v}}{\partial t}+{\rho }_f\frac{\partial \mathbf{q}}{\partial t}=\nabla \cdot {\mathbf{\tau }}_s,$$

$${\rho }_f\frac{\partial \mathbf{v}}{\partial t}+{\psi }_i(t)\ast \frac{\partial \mathbf{q}}{\partial t}=-\nabla p,$$

$${\psi }_i(t)=m_i\delta (t)+\left(\eta /{\kappa }_i\right)H(t),$$

$${\rho }_f\frac{\partial \mathbf{v}}{\partial t}+m_i\frac{\partial \mathbf{q}}{\partial t}+\frac{\eta }{{\kappa }_i}\mathbf{q}=-\nabla p,$$

$$\frac{\partial \mathbf{Q}}{\partial t}+\mathbf{A}\frac{\partial \mathbf{Q}}{\partial x_1}+\mathbf{B}\frac{\partial \mathbf{Q}}{\partial x_2}+\mathbf{C}\frac{\partial \mathbf{Q}}{\partial x_3}=\mathbf{DQ}+\mathbf{f},$$

$$\displaystyle\bm{Q}=\left[\begin{array}[]{ccccccccccccc}\tau_{11},&\tau_{22},&\tau_{33},&\tau_{23},&\tau_{13},&\tau_{12},&p,&v_{1},&v_{2},&v_{3},&q_{1},&q_{2},&q_{3}\end{array}\right]^{T},~{}\bm{A}=\left[\begin{array}[]{c|c}\bm{0}_{7\times 7}&\bm{A}_{11}\\ \hline\cr\bm{A}_{21}&\bm{0}_{6\times 6}\end{array}\right],$$

$$\displaystyle\bm{B}=\left[\begin{array}[]{c|c}\bm{0}_{7\times 7}&\bm{B}_{11}\\ \hline\cr\bm{B}_{21}&\bm{0}_{6\times 6}\end{array}\right],~{}\bm{C}=\left[\begin{array}[]{c|c}\bm{0}_{7\times 7}&\bm{C}_{11}\\ \hline\cr\bm{C}_{21}&\bm{0}_{6\times 6}\end{array}\right],~{}\text{and}~{}\bm{D}=\left[\begin{array}[]{c|c}\bm{0}_{7\times 7}&\bm{0}_{6\times 7}\\ \hline\cr\bm{0}_{6\times 7}&\bm{D}_{22}\end{array}\right],~{}\text{with}$$

$$\displaystyle\bm{A}_{11}=-\left[\begin{array}[]{cccccc}c_{11}^{u}&0&0&\alpha_{1}M&0&0\\ c_{12}^{u}&0&0&\alpha_{2}M&0&0\\ c_{13}^{u}&0&0&\alpha_{3}M&0&0\\ 0&0&0&0&0&0\\ 0&0&c_{55}&0&0&0\\ 0&c_{66}&0&0&0&0\\ -M\alpha_{1}&0&0&-M&0&0\end{array}\right],~{}\bm{A}_{21}=-\left[\begin{array}[]{ccccccc}\displaystyle\frac{m_{1}}{\beta_{1}}&0&0&0&0&0&\displaystyle\frac{\rho_{f}}{\beta_{1}}\\ 0&0&0&0&0&\displaystyle\frac{m_{2}}{\beta_{2}}&0\\ 0&0&0&0&\displaystyle\frac{m_{3}}{\beta_{3}}&0&0\\ -\displaystyle\frac{\rho_{f}}{\beta_{1}}&0&0&0&0&0&-\displaystyle\frac{\rho}{\beta_{1}}\\ 0&0&0&0&0&-\displaystyle\frac{\rho_{f}}{\beta_{2}}&0\\ 0&0&0&0&-\displaystyle\frac{\rho_{f}}{\beta_{3}}&0&0\end{array}\right],$$

$$\displaystyle\bm{B}_{11}=-\left[\begin{array}[]{cccccc}0&c_{12}^{u}&0&0&\alpha_{1}M&0\\ 0&c_{22}^{u}&0&0&\alpha_{2}M&0\\ 0&c_{33}^{u}&0&0&\alpha_{3}M&0\\ 0&0&c_{44}&0&0&0\\ 0&0&0&0&0&0\\ c_{66}&0&0&0&0&0\\ 0&-M\alpha_{2}&0&0&-M&0\end{array}\right],~{}\bm{B}_{21}=-\left[\begin{array}[]{ccccccc}0&0&0&0&0&\displaystyle\frac{m_{1}}{\beta_{1}}&0\\ 0&\displaystyle\frac{m_{2}}{\beta_{2}}&0&0&0&0&\displaystyle\frac{\rho_{f}}{\beta_{2}}\\ 0&0&0&\displaystyle\frac{m_{3}}{\beta_{3}}&0&0&0\\ 0&0&0&0&0&-\displaystyle\frac{\rho_{f}}{\beta_{1}}&0\\ 0&-\displaystyle\frac{\rho_{f}}{\beta_{2}}&0&0&0&0&0\\ 0&0&0&-\displaystyle\frac{\rho_{f}}{\beta_{3}}&0&0&0\end{array}\right],$$

$$\displaystyle\bm{C}_{11}=-\left[\begin{array}[]{cccccc}0&0&c_{13}^{u}&0&0&\alpha_{1}M\\ 0&0&c_{23}^{u}&0&0&\alpha_{2}M\\ 0&0&c_{33}^{u}&0&0&\alpha_{3}M\\ 0&0&c_{44}&0&0&0\\ c_{55}&0&0&0&0&0\\ 0&0&0&0&0&0\\ 0&0&-M\alpha_{3}&0&0&-M\end{array}\right],~{}\bm{C}_{21}=-\left[\begin{array}[]{ccccccc}0&0&0&0&\displaystyle\frac{m_{1}}{\beta_{1}}&0&0\\ 0&0&0&\displaystyle\frac{m_{2}}{\beta_{2}}&0&0&0\\ 0&0&\displaystyle\frac{m_{3}}{\beta_{3}}&0&0&0&\displaystyle\frac{\rho_{f}}{\beta_{3}}\\ 0&0&0&0&-\displaystyle\frac{\rho_{f}}{\beta_{1}}&0&0\\ 0&0&0&-\displaystyle\frac{\rho_{f}}{\beta_{2}}&0&0&0\\ 0&0&-\displaystyle\frac{\rho_{f}}{\beta_{3}}&0&0&0&-\displaystyle\frac{\rho}{\beta_{3}}\end{array}\right],~{}\text{and}$$

$$\displaystyle\bm{D}_{22}=\left[\begin{array}[]{cccccc}0&0&0&\displaystyle\frac{\rho_{f}\eta}{\beta_{1}\kappa_{1}}&0&0\\ 0&0&0&0&\displaystyle\frac{\rho_{f}\eta}{\beta_{2}\kappa_{2}}&0\\ 0&0&0&0&0&\displaystyle\frac{\rho_{f}\eta}{\beta_{3}\kappa_{3}}\\ 0&0&0&-\displaystyle\frac{\rho\eta}{\beta_{1}\kappa_{1}}&0&0\\ 0&0&0&0&-\displaystyle\frac{\rho\eta}{\beta_{2}\kappa_{2}}&0\\ 0&0&0&0&0&-\displaystyle\frac{\rho\eta}{\beta_{3}\kappa_{3}}\end{array}\right],$$

$${\mathbf{Q}}_s\frac{\partial \mathbf{\tau }}{\partial t}=\sum _{i=1}^d{\mathbf{A}}_i\frac{\partial \mathbf{v}}{\partial {\mathbf{x}}_i},$$

$${\mathbf{Q}}_v\frac{\partial \mathbf{v}}{\partial t}=\sum _{i=1}^d{\mathbf{A}}_i^T\frac{\partial \mathbf{\tau }}{\partial {\mathbf{x}}_i}+\mathbf{D}\mathbf{v}+\mathbf{f},$$

$$\displaystyle\bm{Q}_{s}=\left[\begin{array}[]{c|c}\bm{S}&\bm{S\alpha}\\ \hline\cr\bm{\alpha^{T}}\bm{S}&\displaystyle\frac{1}{M}+\bm{\alpha}^{T}S\bm{\alpha}\\ \end{array}\right],~{}~{}\bm{Q}_{v}=\left[\begin{array}[]{cccccc}\rho&0&0&\rho_{f}&0&0\\ 0&\rho&0&0&\rho_{f}&0\\ 0&0&\rho&0&0&\rho_{f}\\ \rho_{f}&0&0&m_{1}&0&0\\ 0&\rho_{f}&0&0&m_{2}&0\\ 0&0&\rho_{f}&0&0&m_{3}\\ \end{array}\right],$$

$$\mathbf{S}$$

$${\mathbf{A}}_1$$

$${\mathbf{A}}_3$$

$$\bm{Q}_{s}^{-1}=\left[\begin{array}[]{ccccccc}c_{11}+M\alpha_{1}^{2}&c_{12}+M\alpha_{1}\alpha_{2}&c_{13}+M\alpha_{1}\alpha_{3}&0&0&0&-M\alpha_{1}\\ c_{12}+M\alpha_{1}\alpha_{2}&c_{11}+M\alpha_{2}^{2}&c_{13}+M\alpha_{2}\alpha_{3}&0&0&0&-M\alpha_{2}\\ c_{13}+M\alpha_{1}\alpha_{3}&c_{13}+M\alpha_{2}\alpha_{3}&c_{33}+M\alpha_{3}^{2}&0&0&0&-M\alpha_{3}\\ 0&0&0&c_{55}&0&0&0\\ 0&0&0&0&c_{55}&0&0\\ 0&0&0&0&0&\displaystyle\frac{c_{11}-c_{12}}{2}&0\\ -M\alpha_{1}&-M\alpha_{2}&-M\alpha_{3}&0&0&0&M\\ \end{array}\right],$$

$$\bm{Q}_{v}^{-1}=\left[\begin{array}[]{ccccccc}\displaystyle\frac{m_{1}}{\rho m_{1}-\rho_{f}^{2}}&0&0&\displaystyle\frac{-\rho_{f}}{\rho m_{1}-\rho_{f}^{2}}&0&0\\ 0&\displaystyle\frac{m_{2}}{\rho m_{2}-\rho_{f}^{2}}&0&0&\displaystyle\frac{-\rho_{f}}{\rho m_{2}-\rho_{f}^{2}}&0\\ 0&0&\displaystyle\frac{m_{3}}{\rho m_{3}-\rho_{f}^{2}}&0&0&\displaystyle\frac{-\rho_{f}}{\rho m_{3}-\rho_{f}^{2}}\\ \displaystyle\frac{-\rho_{f}}{\rho m_{1}-\rho_{f}^{2}}&0&0&\displaystyle\frac{\rho}{\rho m_{1}-\rho_{f}^{2}}&0&0\\ 0&\displaystyle\frac{-\rho_{f}}{\rho m_{2}-\rho_{f}^{2}}&0&0&\displaystyle\frac{\rho}{\rho m_{2}-\rho_{f}^{2}}&0\\ 0&0&\displaystyle\frac{-\rho_{f}}{\rho m_{3}-\rho_{f}^{2}}&0&0&\displaystyle\frac{\rho}{\rho m_{3}-\rho_{f}^{2}}\\ \end{array}\right].$$

$$0<s_{\text{min}}\le {\mathbf{u}}^T{\mathbf{Q}}_s(\mathbf{x})\mathbf{u}\le s_{\text{max}}<\infty$$

$$0<{\tilde{s}}_{\text{min}}\le {\mathbf{u}}^T{\mathbf{Q}}_s^{-1}(\mathbf{x})\mathbf{u}\le {\tilde{s}}_{\text{max}}<\infty$$

$$0<v_{\text{min}}\le {\mathbf{u}}^T{\mathbf{Q}}_v(\mathbf{x})\mathbf{u}\le v_{\text{max}}<\infty$$

$$0<{\tilde{v}}_{\text{min}}\le {\mathbf{u}}^T{\mathbf{Q}}_v^{-1}(\mathbf{x})\mathbf{u}\le {\tilde{v}}_{\text{max}}<\infty$$

$$K=\frac{1}{2}{\mathbf{\tau }}^T{\mathbf{Q}}_s\mathbf{\tau },V=\frac{1}{2}{\mathbf{v}}^T{\mathbf{Q}}_v\mathbf{v}.$$

$${\mathbf{x}}^k={\Phi }^k\hat{\mathbf{x}},$$

$$V_h(D^k)={\Phi }^k\circ V_h(\hat{D}).$$

$$V_h({\Omega }_h)=\underset{k}{⨁}{V}_h(D^k).$$

$$P^N(\hat{D})=\{{\hat{x}}^i{\hat{y}}^j,0\le i+j\le N\},$$

$$P^N(\hat{D})=\{{\hat{x}}^i{\hat{y}}^j{\hat{x}}^k,0\le i+j+k\le N\}.$$

$$\left(\mathbf{g},\mathbf{h}\right)={\int }_{D^k}\mathbf{g}\cdot \mathbf{h}d\mathbf{x}={\int }_{\hat{D}}\mathbf{g}\cdot \mathbf{h}Jd\hat{\mathbf{x}},||\mathbf{g}|{|}_{L^2\Omega }^2=(\mathbf{g},\mathbf{g}{)}_{L^2(D^k)},$$

Full text

1 Introduction

Theories of poroelasticity deal with the physics of elastic wave propagation in porous media. These physics are applicable where pore-filling materials are of interest, such as oil and gas exploration, gas hydrate detection, hydrogeology, seismic monitoring of CO2 storage, and medical imaging. The most popular theory of poroelastic wave propagation was pioneered by Maurice Biot [5, 6]. Biot obtained a dynamical system of equations describing wave propagation in a porous medium, saturated with a single phase fluid with an ability to flow through the pore networks. Biot’s system of equations is divided in two categories, based on the material’s characteristic frequency $(\omega_{c})$ [6, 7, 8], which describes the nature of fluid flow (laminar or non-laminar) induced by wave motion. Biot’s theory of poroelasticity predicts two compressional waves, the fast P wave and slow P wave, and one shear wave.

Biot’s theory of poroelasticity predicts two compressional waves (a fast P wave and slow P wave) and one shear wave. In the fast P wave, the fluid and matrix are locked together and move in-phase. This locking arises through the action of viscous or inertial forces which result from wave-induced fluid flow. Viscous forces arise due to friction between the layers of the viscous fluid (laminar flow), whereas inertial forces result from the momentum of the fluid flow (non-laminar flow). The slow P wave results from relative motion between the fluid in the pores and the matrix. Biot’s system of equations is divided in two categories based on the material’s characteristic frequency $(\omega_{c})$ [6, 7, 8], which describes the nature of fluid flow (laminar or non-laminar) induced by wave motion. Due to the relative motion between the solid and the fluid, a pore boundary layer is formed. The pore boundary layer is related to the viscous-skin depth, which is the distance where the fluid particle velocity attains a certain percentage of its maximum value.

A major cause of energy dissipation in porous media is a dissipative (or drag) “memory” force resulting from the relative motion between the solid skeleton and pore fluids. The dissipative force is incorporated in the equations of motion through a time dependent viscodynamic operator $\psi(t)$. In the low-frequency regime, the dissipative force is linearly proportional to the relative velocity between solid and fluid and $\psi(t)$ is compactly supported. This occurs because the viscous skin depth in the low-frequency regime is greater than the pore size, and thus does not form any pore boundary layer. However, in the high frequency regime, the dominance of inertial forces over viscous forces causes the formation of a pore boundary viscous layer. Thus, in principle, the relaxation mechanism of the medium due to drag force is expressed as a memory kernel and $\psi(t)$ is incorporated into the equation of motion through a convolution operator [9]. Our current work focuses on poroelasticity in the low-frequency regime $(\omega_{c})$. In subsequent work, we will present numerical methods for the broad-band Biot system based on the Johnson-Koplik-Dashen (JKD) model of dynamic permeability [10], which will unify both the high and low frequency regime.

Numerical simulations of the poroelastic wave equation in the low-frequency regime have been previously explored in the literature. The first numerical simulations were performed using the reflectivity method for flat layers [11, 12] and cylindrical structures [13]. Numerical simulation of the poroelastic wave equation using direct or grid based methods, such as finite difference (FD) and pseudo-spectral (PS) methods, dates back to the 1970s [14]. Most of the methods presented in [14] regard pseudospectral [7, 9], staggered pseudospectral [8] and finite-difference methods [15, 16] and are based on 2D structured meshes. Santos and Oreña [17] used a finite-element method to solve the poroelastic wave equation (second-order form) using quadrilateral meshes for spatial discretization. Recent work on the numerical solution of orthotropic poroelasticity is reported by Lemoine et al. [1], using a finite-volume method on structured meshes.

In the present study, we introduce a high-order numerical scheme based on the discontinuous Galerkin method to solve the 3D poroelastic wave equation on unstructured tetrahedral meshes. High order methods provide one avenue towards improving fidelity in numerical simulations while maintaining reasonable computational costs, and methods which can accommodate unstructured meshes are desirable for problems with complex geometries. Among such methods, high order discontinuous Galerkin (DG) methods are particularly well-suited to the solution of time-dependent hyperbolic problems on modern computing architectures [18, 19]. The accuracy of high order methods can be attributed in part to their low numerical dissipation and dispersion compared to low order schemes [20]. This accuracy has made them advantageous for the simulation of electro-magnetic and elastic wave propagation [18, 21].

Since the time-domain wave propagation is described by a hyperbolic system of partial differential equations, an explicit time integration can efficiently be applied. The finite-element methods, when coupled with an explicit time integrator, require the inversion of a global mass matrix, unless special techniques, such as diagonal mass lumping, are applied. Spectral element methods avoid the inversion of the global mass matrix for hexahedral elements by choosing nodal basis functions, which are discretely orthogonal with respect to an under-integrated $L^{2}$ inner product and result in a diagonal mass matrix. In contrast, high order DG methods produce block diagonal mass matrices, which are locally invertible. High order DG methods are often used for seismic simulation (elastic approximation) through the use of simplicial meshes [22, 23, 24].

DG methods impose inter-element continuity of approximate solutions between elements weakly through a numerical flux, of which several choices are common. de la Puente et al. [25] solved the poroelastic wave equation using a local space-time DG method with a Rusanov type flux, which tends to be dissipative. In another study, Ward et al. [26], derived an upwind flux by solving the exact Riemann problem on inter-element boundaries for the 2D isotropic strain-velocity form of the poroelastic wave equations. Recently, Zhan et al. [27] also solved the poroelastic wave equation in 3D (formulated in stress-velocity form) using an upwind flux by solving the exact Riemann problem. The solution of Riemann problem requires diagonalization Jacobian matrices into polarized waves constituents, which is a computationally intensive process for the poroelastic system and does not extend naturally to anisotropic materials. Zhan et al. [27] addresses the complexity of the diagonalization by reducing the rank of Jacobian matrices by using the method of generalized wave impedances. Ye at al. [24] completely avoid the process of diagonalization for the coupled acoustic-elastic wave equation by using a penalty flux based on natural boundary conditions. In this study, we use a similar approach and derive an energy-stable penalty flux for the poroelastic wave equations.

A poroelastic material is defined by the static acoustic properties of fluid, solid and frame exclusively, e.g. fifteen physical properties are required to define an orthotropic medium. The spatial scale of variability of these properties can range from macro (piece-wise constant approximations) to micro-heterogeneities (sub-element variations) in the medium. Most high order methods for wave equations on simplical meshes assume that material coefficients are constant over each element. However, if the media is such that material gradients are non-zero in interior of the an element, piecewise constant approximations can yield inaccurate solutions [28]. This can be circumvented by incorporating sub-element heterogeneities into weighted mass matrices, which recovers a highly accurate and energy-stable DG method [29]. On tetrahedral meshes, this approach requires precomputation and storage of inverses for each local mass matrix, which increases the storage cost and data movement at high orders of approximations. Chan et al. [3] circumvented these storage costs for elastic wave propagation by approximating weighted mass matrices with easily invertible “weight-adjusted” approximations. We extend the same approach to the poroelastic wave equations, where matrix-valued weight functions arise after symmetrization of the system.

In brief, the novelties of our approach are the following:

We obtain a consistent DG weak formulation for the poroelastic wave equations without the diagonalization into polarized wave constituents required for upwind fluxes. 2. 2.

We introduce upwind-like dissipation through simple penalty terms in the numerical flux. 3. 3.

We implement an efficient “weight-adjusted” approximate mass matrix to address micro-heterogeneities at the sub-element level present in poroelastic materials. 4. 4.

We prove stability and high order accuracy of the proposed DG method.

The method is presented for tetrahedral meshes but extends naturally to quadrilateral and hexahedral meshes as well. The outline of the paper is as follows: Section 2 will discuss the system of equations describing the poroelastic wave equation along with interface conditions. Section 3 presents an energy stable formulation with simple penalty fluxes for the symmetric hyperbolic form of the poroelastic wave equations. Section 3 also discusses issues pertaining to storage and inversion of local mass matrices for material coefficients with micro-heterogeneities (sub-element variations). Section 4 incorporates weight-adjusted approximations of weighted $L^{2}$ inner products and mass matrices with matrix-valued weights into DG discretizations of the poroelastic wave equations. Finally, numerical results in Section 5 demonstrate the accuracy of this method for several problems in linear poroelasticity.

2 System of equations describing porelastic waves

Biot’s theory describes wave propagation in a saturated porous medium i.e., a medium made up of a solid matrix (the skeleton or frame) fully saturated with single phase fluid. Biot also assumed that an infinitesimal transformation relates the reference and current states of deformation; thus the displacement, strains, and particle velocity are small. The concept of infinitesimal deformation allows the use of the Lagrangian and Hamilton’s principle to derive the equations governing the propagation of waves in such a medium. In the following two subsections, we will review the constitutive equations, and equations of motions for a poroelastic medium. Readers are advised to refer to Biot’s original papers [5, 30, 6] and [31] for further detail.

2.1 Constitutive equations

The constitutive equations for an inhomogeneous and anisotropic poroelastic medium are expressed as [5, 30, 6] as

[TABLE]

where $c^{u}_{ijkl}$ is the undrained stiffness constant, which can be expressed as $c^{u}_{ijkl}=c_{ijkl}+M\alpha_{i}\alpha_{j}$. Here, $c_{ijkl}$ represents the elements of the stiffness tensor of the solid matrix, with $i,~{}j,~{}k,~{}l=1...6$ denoting indices of stress and strain tensors. The Biot effective stress coefficient $\alpha$ and Biot’s parameter $M$ are derived experimentally from the bulk moduli of the solid and fluid [31]. Finally, $\smash[t]{\epsilon_{kl}^{(m)}}$ signifies the strain with $\zeta$ denoting the variation of fluid content.

2.2 Dynamical equations and Darcy’s law

The dynamic equations describing the wave propagation in anisotropic heterogeneous porous media are given by [5, 30, 6]

[TABLE]

where $\rho=(1-\phi)\rho_{s}+\phi\rho_{f}$ is the bulk density, and $\rho_{s}$ and $\rho_{f}$ are the solid and fluid density, respectively. The variables $\bm{v}$ and $\bm{q}$ are the solid and fluid (relative to solid) particle velocities, respectively, and $\bm{\tau}_{s}$ is total stress (total) tensor.

The generalized dynamic Darcy’s law governing the fluid flow in an anisotropic porous media, is expressed as

[TABLE]

where ${\psi}_{i},~{}i=1...3$ are viscodynamic operators. Here, $``*"$ represents the convolution operation in time. In the low frequency range, i.e., for frequencies lower than $\omega_{c}=\min\left(\displaystyle\frac{\eta\phi}{\rho_{f}T_{i}\kappa_{i}}\right)$, $\psi_{i}$ can be expressed as

[TABLE]

where $m_{i}=T_{i}\rho_{f}/\phi$, with $T_{i}$ being the tortuosity, $\eta$ the fluid viscosity, and $\kappa_{i}$s principal components of the global permeability tensor, while $\delta(t)$ is Dirac’s function and $H(t)$ the Heaviside step function. Substituting (5) into (4), we get

[TABLE]

The convolution operation in (4) represents the dissipation of energy due to fluid-flow induced by the wave motion. As this work focuses on the low-frequency regime, the convolution is approximated by viscosity-dependent damping terms, shown in (6).

2.3 The system of equations in matrix form

To simplify notation, we introduce a matrix form of the system of equations while accounting for anisotropy and heterogeneity. Combining (1), (2), (4) and (6), we get

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\beta_{i}=\rho m_{i}-\rho_{f}^{2}$ and $\bm{f}$ is the vector of forcing terms.

2.4 Symmetric form of system of poroelastic equations

To ensure the stability and convergence of the numerical scheme, we utilize a symmetric form of (7). The symmetric form of (7) for velocity and stress is expressed as

[TABLE]

where

[TABLE]

with

[TABLE]

and $\bm{\tau}=[\tau_{11},~{}\tau_{22},~{}\tau_{33},~{}\tau_{23},~{}\tau_{13},~{}\tau_{12},~{}p]^{T}$, $\bm{v}=[v_{1},~{}v_{2},~{}v_{3},~{}q_{1},~{}q_{2},~{}q_{3}]^{T}$, with

[TABLE]

Explicit expressions for $\bm{Q}_{s}^{-1}$ and $\bm{Q}_{v}^{-1}$ are given by

[TABLE]

[TABLE]

The matrices $\bm{Q}_{s}$ and $\bm{Q}_{v}$ are Hessians of the potential and kinetic energy, respectively, computed with respect to state variables in $\bm{Q}$. Since the energy of the poroelastic system is quadratic positive-definite [9], the Hessians ($\bm{Q_{s}}$, $\bm{Q}_{v}$) are symmetric positive-definite. We also assume the Hessians are bounded as follows

[TABLE]

[TABLE]

[TABLE]

[TABLE]

for $\forall$ $\bm{x}\in\mathbb{R}^{d}$ and $\forall$ $\bm{u}\in\mathbb{R}^{N_{d}}$. Note that in (9), the dissipation matrix $\bm{D}$ is negative-definite.

Using the expressions for $\bm{Q}_{s}$ and $\bm{Q}_{v}$, the kinetic $(K)$ and potential energy $(V)$ of the poroelastic system are expressed as

[TABLE]

3 An energy stable discontinuous Galerkin formulation for poroelastic wave propgation

Energy stable discontinuous Galerkin methods for elastic wave propagation have been constructed based on symmetric formulations of the elastodynamic equations [28], and it is straightforward to extend such discontinuous Galerkin formulations to the symmetric poroelastic system. We assume that the domain $\Omega$ is exactly triangulated by a mesh $\Omega_{h}$ which consists of elements $D^{k}$ which are images of a reference element $\hat{D}$ under the local affine mapping.

[TABLE]

where $\bm{x}^{k}=\{x^{k},y^{k}\}$ for $d=2$ and $\bm{x}^{k}=\{x^{k},y^{k},z^{k}\}$ for $d=3$ denote the physical coordinates on $D^{k}$ and $\hat{\bm{x}}=\{\hat{x},\hat{y}\}$ for $d=2$ and $\widehat{\bm{x}}=\{\widehat{x},\widehat{y},\widehat{z}\}$ for $d=3$ denote coordinates on the reference element. We denote the determinant of the Jacobian of $\Phi^{k}$ as $J$.

Solutions over each element $D^{k}$ are approximated from a local approximation space $V_{h}(D^{k})$, which is defined as composition of mapping $\Phi^{k}$ and reference approximation space $V_{h}(\widehat{D})$

[TABLE]

Subsequently, the global approximation space $V_{h}(\Omega_{h})$ is defined as

[TABLE]

In this work, we will take $V_{h}(\widehat{D})=P^{N}(\widehat{D})$, with $P^{N}(\widehat{D})$ being the space of polynomials of total degree $N$ on the reference simplex. In two dimensions, $P^{N}$ on a triangle is

[TABLE]

and in three dimensions, $P^{N}$ on a tetrahedron is

[TABLE]

The $L^{2}$ inner product and norm over $D^{k}$ is represented as

[TABLE]

where $\bm{g}$ and $\bm{h}$ are real vector-valued functions. Global $L^{2}$ inner products and squared norms are defined as the sum of local $L^{2}$ inner products and squared norms over each elements. The $L^{2}$ inner product and norm over the boundary $\partial D^{k}$ of an element are similarly defined as

[TABLE]

where $J^{f}$ is the Jacobian of the mapping from a reference face $\hat{f}$ to a physical face $f$ of an element.

Let $f$ be a face of an element $D^{k}$ with neighboring element $D^{k,+}$ and unit outward normal $\bm{n}$. Let $u$ be a function with discontinuities across element interfaces. We define the interior value $u^{-}$ and exterior value $u^{+}$ on face $f$ of $D^{k}$

[TABLE]

The jump and average of a scalar function $u\in V_{h}(\Omega_{h})$ over $f$ are then defined as

[TABLE]

Jumps and averages of vector-valued functions $\bm{u}\in\mathbb{R}^{m}$ and and matrix-valued functions $\tilde{\bm{S}}\in\mathbb{R}^{m\times n}$ are defined component-wise.

[TABLE]

We can now specify a DG formulation for poroelastic wave equation (7). The symmetric hyperbolic system in (7) readily admits a DG formulation based on a penalty flux [32, 28]. For the symmetric first order poroelastic wave equation, the DG formulation in strong form is given as

[TABLE]

for all $\bm{h},~{}\bm{g}\in V_{h}(\Omega_{h})$. Here, $\bm{A}_{n}$ is the normal matrix defined on a face $f$ of an element

[TABLE]

The factors $\alpha_{\bm{\tau}}$ and $\alpha_{\bm{v}}$ are penalty parameters and defined on element interfaces. We assume that $\alpha_{\bm{\tau}},~{}\alpha_{\bm{v}}\geq 0$ and are piecewise constant over each shared face between two elements. These penalty constants can be taken to be zero, which results in a non-dissipative central flux, while $\alpha_{\bm{\tau}},~{}\alpha_{\bm{v}}>0$ results in energy dissipation similar to the upwind flux [18]. The stability of DG formulations are independent of the magnitude of these penalty parameters. However, a naive choice of these parameters will result in a stiffer semi-discrete system of ODEs and necessitates a smaller time under explicit time integration schemes. In this work, we take $\alpha_{\tau},\alpha_{v}=O(1)$ unless stated otherwise.

Remark**.**

For several problems in Section 5, $\bm{f}$ is taken as a scaled point source perturbation or Dirac delta function $\bm{f(x)}=\bm{\beta}(\bm{x})\Delta(\bm{x-x_{0}})$ (with $\bm{\beta(x)}\in\mathbb{R}$). In this setting, $\bm{f}$ is not $L^{2}$ integrable and thus $\left(\bm{f},\bm{g}\right)_{L^{2}(D^{k})}$ may not be well-defined. In such cases, we commit a variation crime and evaluate its contribution as

[TABLE]

3.1 Physical interpretation of central flux terms

Apart from providing the global solution of the system, the central flux terms corresponding to the velocity $(\bm{A}_{n}\llbracket\bm{v}\rrbracket)$ and the stress $(\bm{A}_{n}^{T}\llbracket\bm{\tau}\rrbracket)$ also enforce natural boundary conditions at the interface between two poroelastic media. These conditions are related to the phenomena of reflection, refraction, and diffraction of waves in the presence of inhomogeneities and interfaces. The component-wise expressions of the central flux terms $\bm{A}_{n}^{T}\llbracket\bm{\tau}\rrbracket$ and $\bm{A}_{n}\llbracket\bm{v}\rrbracket$ correspond to

[TABLE]

where the “$\pm$” convention is determined by the outward interface normal, $\bm{n}$. In this work, the outer normal vector is assumed point in the direction of the “$+$” side of the interface.

The physical significance of the above conditions are:

(12) describes the continuity of traction across the interface. 2. 2.

(13) states that the medium across the interface stays intact. 3. 3.

(14) implies that all fluid entering the interface should exit the other side. 4. 4.

(15) assumes perfect hydraulic conductivity across the interface, providing continuity of pore-fluid pressure.

These interface conditions are also consistent with the open-pore boundary conditions of [33, 34] between two different poroelastic media.

3.2 Boundary conditions

In many applications, the boundary condition for the top surface of a domain is assumed to be a free surface (stress free), with the remaining surfaces taken to be absorbing boundaries. We impose boundary conditions on the DG formulation by choosing appropriate exterior values which result in modified boundary numerical fluxes. Boundary conditions on the normal component of the stress can be imposed by noting that that the numerical flux contain the term $\llbracket\bm{A}_{n}^{T}\bm{\tau}\rrbracket=\llbracket\widetilde{\bm{S}}\bm{n}\rrbracket$ with

[TABLE]

For a face which lies on the top surface of the domain, the free surface boundary or zero traction boundary conditions can be imposed by setting

[TABLE]

For problems which require the truncation of infinite or large domains, basic absorbing boundary conditions can be imposed by setting

[TABLE]

In addition to the above boundary conditions, more accurate absorbing boundary conditions can be also imposed using perfectly matching layers (PML) [35] or high order absorbing boundary conditions (HABC) [36, 37]. However, the implementation of such boundary conditions results in an augmented system of PDEs which can become very expensive. For example, the implementation of PML for the poroelastic wave equations results in a system of thirty PDEs [35].

In all cases, the boundary conditions are imposed by computing the numerical fluxes based on the modified jumps. This imposition guarantees energy stability for free surface and absorbing boundary conditions.

3.3 Energy stability

One can show that the DG formulation in (11) is energy stable in the absence of external forces $(\bm{f}=0)$, and free-surface and absorbing boundary conditions. Integrating by parts the velocity equation in (11) gives

[TABLE]

Taking $\left(\bm{h},\bm{g}\right)=(\bm{\tau},\bm{v})$ and adding both equations together yields

[TABLE]

where the term

[TABLE]

is the total energy of the system. Let $\Gamma_{h}$ denote the set of unique faces in $\Omega_{h}$ and let $\Gamma_{\bm{\tau}}$, $\Gamma_{\text{abc}}$ denote boundaries where free-surface and absorbing boundary conditions are imposed, respectively. We separate surface terms into contributions from interior shared faces and from boundary faces. On interior shared faces, we sum the contributions from the two adjacent elements to yield

[TABLE]

where $\bm{v}^{T}\bm{D}\bm{v}<0$, since $\bm{D}$ is a negative semi-definite matrix. For faces which lie on $\Gamma_{\bm{\tau}}$ , $\bm{A}_{n}^{T}=-2\bm{A}_{n}^{T}\bm{\tau}^{-},~{}\bm{A}_{n}^{T}\left\{\!\{\tau\}\!\right\}=0$ and $\llbracket\bm{v}\rrbracket=0$ yielding

[TABLE]

Finally, for faces in $\Gamma_{\text{abc}}$, we have $\bm{A}_{n}^{T}\left\{\!\{\tau\}\!\right\}=\displaystyle\frac{1}{2}\bm{A}_{n}^{T}\tau{{}^{-}}$, $\bm{A}_{n}^{T}\llbracket\tau\rrbracket=-\bm{A}_{n}^{T}\bm{\tau}^{-}$ and $\llbracket\bm{v}\rrbracket=-\bm{v^{-}}$, yielding

[TABLE]

Combining contributions from all faces and dissipation in the system yields the following result:

Theorem 1**.**

The DG formulation in (11) is energy stable for $\alpha_{\bm{\tau}},\alpha_{\bm{v}}\geq 0$ such that

[TABLE]

Since $\bm{Q}_{s}$ and $\bm{Q}_{v}$ are positive definite, the left hand side of (1) is an $L^{2}$-equivalent norm on $(\bm{\tau},\bm{v})$ and Theorem 1 implies that magnitude of the DG solution is non-increasing in time. This also shows that dissipation is present for positive penalization parameters, i.e. ${\alpha_{\bm{\tau}}},{\alpha_{\bm{v}}}>0$.

3.4 The semi-discrete matrix system for DG

Let ${\{\phi_{i}\}}_{i=1}^{N_{p}}$ be a basis for $P^{N}\left(\widehat{D}\right)$. In our implementation, we use nodal basis functions located at Warp and Blend interpolation points [18], which are defined implicitly using an orthogonal polynomial basis on the reference simplex. We define the reference mass matrix $\widehat{\bm{M}}$ and the physical mass matrix $\bm{M}$ for an element $D^{k}$ as

[TABLE]

For affine mappings, $J$ is constant and $\bm{M}=J\bm{\widehat{M}}$. We also define weak differentiation matrices $\bm{S}_{k}$ and face mass matrices $\bm{M}_{f}$ such that

[TABLE]

where $J^{f}$ is the Jacobian of the mapping from the reference face $\widehat{f}$ to $f$. For affinely mapped simplices, $J^{f}$ is also constant and $M_{f}=J^{f}\widehat{\bm{M}}_{f}$, where the definition of the reference face mass matrix $\widehat{\bm{M}}_{f}$ is analogous to the definition of the reference mass matrix $\widehat{\bm{M}}$.

Finally, we introduce weighted mass matrices. Let $w(\bm{x})\in\mathbb{R}$ and $\bm{W(x)}\in\mathbb{R}^{m\times n}$. Then, scalar and matrix-weighted mass matrices $\bm{M}_{w}$ and $\bm{M_{W}}$ are defined as

[TABLE]

where $\bm{M_{W_{i,j}}}$ is the scalar weighted mass matrix weighted by the $(i,j)^{\text{th}}$ element of $\bm{W}$. Note that $\bm{M}_{w},\bm{M}_{\bm{W}}$ are positive definite if $w(x),\bm{W}$ are pointwise positive definite.

Local contributions to the DG variational form may be evaluated in a quadrature-free manner using matrix-weighted mass matrices as defined above. Let $\bm{\Sigma_{i}},~{}\bm{V}_{i}$ denote vectors containing degrees of freedom for solutions components $\bm{\tau}_{i}$ and $\bm{v}_{i}$, such that

[TABLE]

Then, the local DG formulation can be written as a block system of ordinary differential equations by concatenating $\bm{\Sigma}_{i},~{}\bm{V}_{i}$ into single vectors $\bm{\Sigma}$ and $\bm{V}$ and using the Kronecker product $\otimes$

[TABLE]

where $\bm{F}_{v}$ and $\bm{F}_{\tau}$ denote the degrees of freedom for the velocity and stress numerical fluxes.

In order to apply a time integrator, we must invert $\bm{M}_{\bm{Q}_{s}}$ and $\bm{M}_{\bm{Q}_{v}}$. While the inversion of $\bm{M}_{\bm{Q}_{s}}$ and $\bm{M}_{\bm{Q}_{v}}$ can be parallelized from element to element, doing so typically requires either the precomputation and storage of the dense matrix inverses or on-the-fly construction and solution of a large dense matrix system at every time step. The former option requires a large amount of storage, while the latter option is computationally expensive and difficult to parallelize among degrees of freedom. This cost can be avoided when ${\bm{Q}_{s}}$ and $\bm{Q}_{v}$ are constant over an element $D^{k}$. In this case, $\bm{M}_{\bm{Q}_{s}}$ reduces to

[TABLE]

Similarly ${\bm{M}^{-1}_{\bm{Q}_{v}}}$ and ${\bm{M}^{-1}_{\bm{D}}}$ can be expressed as $\bm{M}^{-1}_{\bm{Q}_{v}}={\bm{Q}_{v}}^{-1}\otimes\left(\displaystyle\frac{1}{J}\bm{\widehat{M}^{-1}}\right)$, and $\bm{M}^{-1}_{\bm{D}}={\bm{D}}^{-1}\otimes\left(\displaystyle\frac{1}{J}\bm{\widehat{M}^{-1}}\right)$ respectively. Applying these observations to (19) and (20) yields the following sets of local ODEs over each element

[TABLE]

where we have introduced the differentiation matrix $\bm{D}_{i}=\bm{M}{{}^{-1}}\bm{S}_{i}$ and lift matrix $\bm{L}_{f}=\bm{M}^{-1}\bm{M}_{f}$. For affine elements, both derivative and lift matrices are applied using products of geometric factors and reference derivative and lift matrices.

Unfortunately, if $\bm{Q}_{s}$ and $\bm{Q}_{v}$ varies spatially with in the element, then above approach can no longer be used to invert $\bm{Q}_{s}$ and $\bm{Q}_{v}$. Here, we follow the approach of [3], where $\bm{M}_{\bm{Q}_{s}},\bm{M}_{\bm{Q}_{v}}$ are replaced with weight-adjusted approximations. These approximations are low storage, simple to invert, and yield an energy stable and high order accurate DG method to approximate the matrix-weighted $L^{2}$ inner product (and corresponding matrix-weighted mass-matrices $\bm{Q}_{s}$ and $\bm{Q}_{v}$).

We also note that, material coefficients $\bm{Q_{s}},~{}\bm{Q}_{v}$ appear only the left hand side of (7). The right hand side of (7) is equivalent to the discretization of a constant coefficient system. This provides additional advantages in that the right hand side can be evaluated using efficient techniques for DG discretizations of constant-coefficient problems [38, 39].

4 Weight-adjusted discontinuous Galerkin (WADG) formulation for poroelastic wave propagation

We wish to apply weight-adjusted approximations to avoid the inversion of $\bm{M}_{\bm{Q}{{}_{s}}}$ and $\bm{M}_{\bm{Q}{{}_{v}}}$. Weight-adjusted inner products are high order approximations of weighted $L^{2}$ inner products. These weight-adjusted inner products result in weight-adjusted mass matrices whose inverses approximate the inverses of weighted $L^{2}$ mass matrices. We briefly review scalar and matrix-valued weight-adjusted mass matrices in the following section. Matrix weight-adjusted inner products and weight-adjusted approximations with matrix weights are discussed in more detail in A and B.

4.1 Approximation of weighted mass matrix inverses

The advantage of weighted-adjusted inner products is that the corresponding weight-adjusted mass matrices are straightforward to invert. For scalar weights, weight-adjusted mass matrices are given as follows [3]

[TABLE]

By evaluating $\bm{M}_{1/w}$ in a matrix-free fashion using a sufficiently accurate quadrature rule, the inverse of the weight-adjusted mass matrix $\bm{M}^{-1}\bm{M}_{1/w}\bm{M}^{-1}$ yields a low storage implementation. Let $\bm{\hat{x}}_{i},\bm{\hat{w}}_{i}$ denote quadrature points and weights on the reference element, and let $\bm{V}_{q}$ denote the generalized Vandermonde matrix

[TABLE]

whose columns correspond to evaluations of basis functions at quadrature points. Then, for affine elements, $M=J\widehat{\bm{M}}=J\bm{V}_{q}^{T}\text{diag}(\hat{\bm{w}}_{i})\bm{V}_{q}$, where $\widehat{\bm{M}}$ is the reference mass matrix and $J$ is the determinant of the Jacobian of the reference-to-physical mapping, which is constant for affine mappings. Additionally,

[TABLE]

where $w(\hat{\bm{x}}_{i})$ denotes the evaluation of the weight function $w(\bm{x})$ at quadrature points. Thus, for a vector $\bm{u}$, the inverse of the weight-adjusted mass matrix can be applied as follows

[TABLE]

where $\bm{P}_{q}=\bm{\widehat{M}}^{-1}\bm{V}_{q}^{T}\text{diag}(\widehat{\bm{w}}_{i})$ is the quadrature-based $L^{2}$ projection operator on the reference element.

For weight-adjusted inner products with matrix-valued weights, the corresponding weight-adjusted mass matrices approximate weighted $L^{2}$ mass matrices in a similar fashion

[TABLE]

where $\bm{M_{W}}$ is the matrix-weighted mass matrix defined in (18). We note that $\bm{M_{W}}$ can be applied in a quadrature based fashion component by component.

In the context of DG with explicit time-stepping, the factor of $\widehat{\bm{M}}^{-1}$ can be pre-multiplied into the right hand side (the spatial discretization). The application of the weight-adjusted mass matrices then requires only two reference matrices $\bm{V}_{q}$ and $\bm{P}_{q}$ and the values of the weight function at quadrature points $w(\bm{\widehat{x}}_{i})$. In this study the number of quadrature points is $O(N^{3})$ [40]. The overall storage cost for applying weight-adjusted mass matrices using the above implementation is $O(N^{3})$ per element, while the pre-computation and storage of DG operators involving inverses of weighted mass matrices require $O(N^{6})$ storage per element.

We derive a weight-adjusted DG method by replacing the $L^{2}$ inner products in the left hand side of the DG formulation (11) with weight-adjusted approximations. The right hand side of the WADG formulation is similar to the right hand side of the DG formulation (11), and preserves a variant of the energy stability in Theorem 1 i.e.,

[TABLE]

where $T_{{\bm{Q}_{s}}^{-1}}$ are $T_{{\bm{Q}_{v}}^{-1}}$ are weighting operators defined in B.

We replace the weighted $L^{2}$ mass matrices in (19) and (20) by their weight-adjusted approximations. Inverting these weight-adjusted mass matrices yields the following local system of ODEs for $\bm{\Sigma}$ and $\bm{V}$:

[TABLE]

Matrices $\left(\bm{I}\otimes\bm{M}^{-1}\right)\bm{M}_{\bm{Q}_{s}^{-1}}$ and $\left(\bm{I}\otimes\bm{M}^{-1}\right)\bm{M}_{\bm{Q}_{v}^{-1}}$ are applied in a matrix-free fashion using reference element matrices and values of $\bm{Q}_{s}^{-1}$ and $\bm{Q}_{v}^{-1}$ at quadrature points in (8) and (9). The convergence analysis of the numerical scheme is performed in C.

5 Numerical experiments

In this section, we present several numerical experiments validating the stability and accuracy of proposed method in two and three dimensions. The convergence of the new DG formulation in piecewise constant isotropic poroelastic media is confirmed. Finally, the method is applied to problems with anisotropy and micro-heterogeneities (sub-element variations).

In all experiments, we follow [41] and compute the application of weight-adjusted mass matrices using a quadrature which is exact for polynomials of degree $(2N+1)$. Time integration is performed using the low-storage 4${}^{\text{th}}$ order five-stage Runge-Kutta scheme of Carpenter and Kennedy [42], and the time step is chosen based on the global estimate

[TABLE]

where $C_{N}=O(N^{2})$ is the order-dependent constant in the surface polynomial trace inequality [43] and $C_{\rm CFL}$ is a tunable global CFL constant. This estimate is derived by bounding the eigenvalues of the spatial DG discretization matrix appearing in the semi-discrete system of ODEs. This choice of $dt$ is very conservative as it is derived based on an upper bound on the spectral radius.

5.1 Spectra and choice of penalty parameter

We first verify the energy stability of proposed DG formulation. Let $\bm{A}_{h}$ denote the matrix induced by the global semi-discrete DG formulation, such that time evolution of the solution $\bm{\tau},\bm{v}$ is governed by

[TABLE]

where $\bm{Q}$ denotes a vector of degrees of freedom for $(\bm{\tau},\bm{v})$. We show in Figure 1 eigenvalues of $\bm{A}_{h}$ for $\alpha_{\bm{\tau}}=\alpha_{\bm{v}}=0$ and $\alpha_{\bm{\tau}}=\alpha_{\bm{v}}=1$ with material parameters of isotropic sandstone (given in Table 1). The discretization parameters are $N=3$ and $h=1/2$. In both cases, the largest real part of any eigenvalues is $O(10^{-14})$, which suggests that the semi-discrete scheme is indeed energy stable.

For practical simulations, the choice of $\alpha_{\bm{\tau}},\alpha_{\bm{v}}$ remains to be specified. Taking $\alpha_{\bm{\tau}},\alpha_{\bm{v}}>0$ results in damping of of under-resolved spurious components of the solutions. However, a naive selection of $\alpha_{\bm{\tau}},\alpha_{\bm{v}}$ can result in an overly restrictive time-step restriction for stability. A guiding principle for determining appropriate values of the penalty parameters $\alpha_{\bm{\tau}},\alpha_{\bm{v}}$ is to ensure that the spectral radius is the same magnitude as the case when $\alpha_{\bm{\tau}}=\alpha_{\bm{v}}=0$. For example, the spectral radius of $\bm{A}_{h}$, $\rho(\bm{A}_{h})$ is $18.0034$ for $\alpha_{\bm{\tau}},\alpha_{\bm{v}}=0$ which is $O(N^{2}/h)$. The spectral radius $\rho(\bm{A}_{h})$ is $19.1217$ for $\alpha_{\bm{\tau}},\alpha_{\bm{v}}=0.5$, while the spectral radius for $\alpha_{\bm{\tau}},\alpha_{\bm{v}}=1$ is $\rho(\bm{A}_{h})=44.44$. Since the maximum stable timestep is proportional to the spectral radius, taking $\alpha_{\bm{\tau}},\alpha_{\bm{v}}=1$ in this case results in a more restrictive CFL condition. This phenomena is related to observations in [32] that large penalty parameters result in extremal eigenvalues of $\bm{A}_{h}$ with very large negative real parts. The optimal choice of scaling also depends on the media heterogeneities on elements adjacent to an interface.

In all following experiments, we use ${\alpha}_{\bm{\tau}}={\alpha}_{\bm{v}}=1$ unless specified otherwise.

5.2 Analytic solution

Next, we study the accuracy and convergence of the our DG method for a plane wave propagating in an isotropic porous sandstone with material properties given in Table 1 (Column 4). Unless otherwise stated, we report relative $L^{2}$ errors for all components of the solution $\bm{U}=(\bm{\tau},\bm{v})$

[TABLE]

5.3 Plane wave in a poroelastic medium

The analytical solution to (7) for a plane wave is given as

[TABLE]

where $\bm{Q}_{n}^{0}$ is the initial amplitude vector of stress and velocity components; $\omega$ are wave frequencies; $\bm{k}=(k_{x},k_{y},k_{z})$ is the wave-number vector. To achieve realistic poroelastic behavior, we superimpose three plane waves, of the form given by (24), corresponding to a fast P-wave, an S-wave and a slow P-wave.

Now, we briefly describe how we determine the the wave frequencies $\omega$. Substituting (24) into (7) yields

[TABLE]

Solving the three eigenvalues problem for in (24) for each wave mode $l$ yields in matrix of right eigenvectors $(R^{(l)}_{mn})$ and eigenvalue $(\omega_{l})$. Following [44, 25], the solution of (7) can be constructed as

[TABLE]

where $\gamma_{n}^{(l)}$ is a amplitude coefficient with $\gamma_{1}^{(1)}=\gamma_{2}^{(2)}=\gamma_{4}^{(3)}=100$.

For completeness, we perform convergence analyses for inviscid and viscid media separately.

5.3.1 Inviscid case ($\eta=0$)

The error is computed for an inviscid brine-filled ($\eta=0)$ isotropic sandstone. The inviscid nature of the fluid implies that $\bm{D}=\bm{0}$. We show in Figure 2 the $L^{2}$ errors computed at T=1 and CFL=1, using uniform triangular meshes constructed by bisecting an uniform mesh of quadrilaterals along the diagonal. Figure 2a and 2b show the convergence plot using the central flux $(\alpha_{\bm{v}}=\alpha_{\bm{\tau}}=0)$ and penalty flux $(\alpha_{\bm{v}}=\alpha_{\bm{\tau}}=1)$ respectively. For $N=1,...,5$, $O(h^{N+1})$, rates of convergence are observed. We note that for $N=4$ and $N=5$, we observe results for both fluxes which are better than the $4^{\text{th}}$ order accuracy of our time-stepping scheme. This is most likely due to the benign nature of the solution in time and the choice of time step (23) which scales as $O(h/N^{2})$. For $N=4,5$, the results of Figure 2 suggest that the resulting time step is small enough such that temporal errors of $O(dt^{4})$ are small relative to spatial discretization errors of $O(h^{N+1})$.

5.3.2 Viscid case ($\eta\neq 0$)

The poroelastic formulations used by Carcione [9] and de la Puente et al. [25] produce a stiff system of ODEs due to the presence of the dissipative matrix $\bm{D}$. Consequently, Carcione used a Strang’s 2${}^{\text{nd}}$ order operator splitting approach [45] to avoid small $\Delta t$, while de la Puente et al. [25] circumvented the effect of stiffness using a local implicit time-stepping approach to achieve the convergence rate of $O(h^{N+1})$. The poroelastic formulation and numerical discretization used in this study do not appear to produce equally stiff systems of ODEs, and Strang splitting is not required. However, for a complete analysis, the convergence of the proposed DG formulation in viscid case is studied using both a unified DG scheme (where we incorporate the dissipative matrix $\bm{D}$ directly into the explicit time-stepping scheme) as well as a second order Strang operator splitting [45].

Operator splitting separates the dissipative term from the conservative term at each time step. We solve the dissipative part $(\text{stiff})$ part of the system analytically and solve the conservative part by using the proposed DG formulation. We rewrite the system with non-zero forcing function as follows

[TABLE]

The formal solution of (27) is given as

[TABLE]

where $\exp(\bm{A}_{h}t)$ is an evolution operator.

Using the operator splitting approach, the propagation matrix can be partitioned as

[TABLE]

where the subscript $c$ indicates the matrix representing the conservative part of the system, and the subscript $d$ indicates the diffusive matrix, representing low order terms coupled with the viscosity $(\eta)$. Subsequently, the evolution operator can be expressed as

[TABLE]

By using the product formula, (30) is expressed as

[TABLE]

It can be shown that (30) is second order accurate in $dt$. Thus, (31) allows us to solve the stiff part separately. The solution of the stiff part of the system can be derived analytically and is shown in Appendix A.

We show in Figure 3 the $L^{2}$ errors for viscid case, computed at T=1 and CFL=1 using uniform triangular meshes. Figure 3a and 3b show convergence plots using the central flux $(\alpha_{\bm{v}}=\alpha_{\bm{\tau}}=0)$ and penalty flux $(\alpha_{\bm{v}}=\alpha_{\bm{\tau}}=1)$, respectively. For $N=1,...,5$, $O(h^{N+1})$ rate of convergence are observed when using the penalty flux $\alpha_{\bm{v}}=\alpha_{\bm{\tau}}=1$. When using a central flux $(\alpha_{\bm{v}}=\alpha_{\bm{\tau}}=0)$, we observe a so-called “even-odd” pattern [18, 29], with convergence rate $O(h^{N})$ for odd $N$ and between $O(h^{N+1/2})$ and $O(h^{N+1})$ for even $N$.

Similar to the inviscid case, for $N=4$ and $N=5$, we observe results for both fluxes which are better than the $4^{\text{th}}$ order accuracy of our time-stepping, due to the fact that temporal error are small relative to spatial discretization error. Furthermore, the spectral radius $\rho(\bm{A}_{h})$ for $N=4$ and $h=1/2$ is 25.91 which is of $O(N^{2}/h)$. Since ${N^{2}}/{h}>||\bm{D}=1.667||$, the dissipative term does not increase the stiffness of the high order DG scheme.

We show in Figure 4 the $L^{2}$ errors for the viscid case using Strang splitting, computed at T=1 and CFL=1. Figure 4a and 4b show convergence plots using the central flux $(\alpha_{\bm{v}}=\alpha_{\bm{\tau}}=0)$ and penalty flux $(\alpha_{\bm{v}}=\alpha_{\bm{\tau}}=1)$, respectively. Observations on convergence are similar to those in the unified DG scheme. Here, it is worth noting that for $N=5$ the convergence rate is $5.60$, which does not match the optimal theoretical rate of $O(h^{N+1})$. This is may be due to the second order accuracy of the splitting scheme resulting in larger temporal errors relative to the spatial discretization errors $O(h^{N+1})$.

5.4 Application examples

We next demonstrate the accuracy and flexibility of the proposed DG for several application-based problems in linear poroelasticity with micro heterogeneities and anisotropy. All computations are done using penalty parameters $\alpha_{\bm{v}}=\alpha_{\bm{\tau}}=1$ unless specified otherwise. In the subsequent sections, the field $\bm{b}$ represents the center of mass particle velocity vector [46], which is expressed as

[TABLE]

In subsequent simulations, the forcing is applied to both the $z-$ component of stress $\tau_{zz}$ and the fluid pressure $p$ by a Ricker wavelet point source

[TABLE]

where $\bm{x}_{0}$ is the position of the point source and $f_{0}$ is the central frequency.

In the following simulations, three types of poroelastic waves are observed: the fast P wave, S wave, and slow P wave. The fast P wave has solid and fluid motion in phase, while the slow P wave (Biot’s mode) has the solid and fluid motion out of phase with one another. At low frequencies, the slow P wave is diffusive in character as viscous forces dominate over inertial forces. However, at high frequencies, the dominance of inertial forces over viscous forces results in the propagation of the slow P wave.

5.4.1 Orthotropic sandstone

To illustrate the effect of anisotropy on poroelastic wave propagation, we perform a computational experiment in orthotropic sandstone with material properties given in Table 1. The size of the computational domain is $18.25~{}\text{m}\times 18.25~{}\text{m}$. The domain is discretized with uniform triangular element with a minimum edge length of $5~{}\text{cm}$. Figures 5(a)-(d) represent the $x-$ and $z-$ components of the center of mass particle velocity of the orthotropic sandstone, where (a) and (b) correspond to the inviscid case $(\eta=0)$, and (c) and (d) to the viscid case $(\eta\neq 0)$. The central frequency of the forcing function is $f_{0}=3730~{}\text{Hz}$, and polynomials of degree $N=4$ are used. The propagation time is $1.56~{}\text{ms}$. Three events can be observed: the fast P mode ($\text{P}_{\text{f}}$, outer wavefront), the shear wave (S, middle wavefront), and the slow P mode ($\text{P}_{\text{s}}$, inner wavefront). In the viscid case, the slow mode diffuses faster and the medium behaves almost as a single phase medium.

5.4.2 Epoxy-glass material

Snapshots of the $x-$ and $z-$ components of the center of mass particle velocity in the epoxy-glass porous medium are shown in Figure 6. Figures 6(a) and (b) correspond to the inviscid case $(\eta=0)$, and (c) and (d) to the viscid case $(\eta\neq 0)$. The central frequency is $f_{c}=3135~{}\text{Hz}$ and the propagation time is $1.8~{}\text{ms}$. It is worth noting the cuspidal triangles of S and $\text{P}_{\text{s}}$, which are phenomena typical in anisotropic materials. At $45^{o}$, the polarization of the $\text{P}_{\text{s}}$ mode wave is almost horizontal, which confirms the results shown in Figure 3(b) of [9].

5.4.3 Isotropic-isotropic and isotropic-anisotropic layered model

In this example, we illustrate the effect of an interface between two layers of porous media. We constructed two models, each with two layers. In the first model, the top and bottom layer are isotropic and made of shale and sandstone. In second model, the sandstone is replaced by orthotropic sand stone. The detailed material properties are given in Table 1. These models are considered to be filled with brine $(\eta=0)$. The size of the computational domain is $1400~{}\text{m}\times 1500~{}\text{m}$ in the $x$ and $z$ directions, respectively. The minimum edge size of the triangular elements used to mesh the domain is $8~{}\text{m}$. The point source is located at $(750~{}\text{m},900~{}\text{m})$ with a Ricker wavelet of frequency 45 Hz. The propagation time for isotropic and anisotropic model are $0.25~{}\text{s}$ and $0.22~{}\text{s}$, respectively. The simulation is performed using polynomials of degree $N=4$. Snapshots of the $z$ component of the center of mass particle velocity are shown in Figures 7a and 7b for inviscid $(\eta=0)$ and viscid fluid $(\eta\neq 0)$, respectively. Figure 7 clearly shows the direct, reflected, and transmitted wavefront, corresponding to all three modes. The slow P wave is more prominent in the shale. The effect of anisotropy on all three modes is clearly seen as wavefronts moves with different phase velocities.

Here, we illustrate the effect of viscosity of fluid on a slow P wave using an isotropic two layer model with the same discretization parameters as those used in Figure 7a. We compute numerical solutions for two cases. In first case, both layers are assumed to be brine filled with $(\eta=0)$ and are represented by inviscid-inviscid model. However, in the second case we assume that the bottom layer is filled with a viscous fluid and represented as inviscid-viscid model. Figure 8a shows the snapshot of the $z$ component of the center of mass particle velocity for the inviscid-inviscid model, which shows a direct, reflected, and transmitted slow P wave. On the other hand, Figure 8b shows a snapshot of the $z$ component of the center of mass particle velocity for the inviscid-viscid model, which shows only a direct and reflected slow P wave. The transmitted slow P wave is missing in the inviscid-viscid model, which confirms Biot’s observation on slow P waves for the low frequency regime [5].

5.4.4 Piecewise constant vs sub-cell heterogeneities

To address the effect of micro-heterogenities on poroelastic wave propagation, we illustrate the efficacy of WADG method over a traditional DG method using piecewise constant coefficients on each element. We modulate the material properties contained in $\bm{Q}^{-1}_{s}$ and $\bm{Q}^{-1}_{v}$ with

[TABLE]

where $\bm{x}$ and $\bm{y}$ are coordinates of of cubature nodes. We construct piecewise constant approximations to $\bm{Q}^{-1}_{s}$ and $\bm{Q}^{-1}_{v}$ by taking the local average over each element.

Figure 9a and 9b show piecewise smooth and high order heterogeneous approximations of bulk density, respectively. To demonstrate the effect of subelement variations on wave propagation, we first compute the solutions approximated by polynomials of degree $N=2$ with $K_{\rm 1D}=64$, where $K_{\rm 1D}$ denotes the number of elements in one direction. These solutions are shown in Figures 10a (piecewise constant) and 10b (WADG). The solutions are similar, as the lower resolution due to the piecewise constant approximation of media is offset by the larger number of elements. Next, we compute the solutions for $N=4$ and $K_{1D}=32$. The effects of subelement variations are now clearly marked by spurious reflections in Figures 10c and 10d. These effects are more prominent for higher orders as shown in Figures 10e and f, which display the piecewise constant and WADG solutions for $N=8$ and $K_{1D}=16$.

5.4.5 Validation of numerical solution

Here, we show comparisons between numerical solutions of the poroelastic wave equation obtained from the proposed DG method and from the spectral element method [4]. Morency and Tromp [4] used the spectral element method to solve the poroelastic wave equation in second order displacement form for 2D isotropic porous materials. First, we perform a validation study in a homogeneous porous medium. The domain dimensions are $1000~{}\text{m}\times 1000~{}\text{m}$. The domain is discretized with triangular and quadrilateral elements with minimum edge length $5~{}\text{m}$ for the DG and spectral element method, respectively.

Mornecy and Tromp implemented a pure dilatational explosive source as an external forcing function. To match the assumption of the dilatational explosive source, we have set the forcing functions corresponding to $\tau_{11}$ and $\tau_{22}$ to be equal, while the forcing function corresponding to $\tau_{12}$ is set to zero since a pure dilatational source does not radiate shear waves. It was also assumed in [4] that the fluid and the solid contain equal energy, and thus the forcing function corresponding to the pressure variable $p$ is the same as for $\tau_{11}$ and $\tau_{12}$.

The time domain response of the forcing function is a Ricker wavelet with a central frequency of $f_{c}=30~{}\text{Hz}$. The material properties of the medium are given in column 2 of Table 2. In [3], the material properties were represented using Lame constants $\kappa_{fr}$, $\mu_{fr}$. These can be converted to compliance coefficients as follows:

[TABLE]

Figure 11a shows a snapshot of the solid particle velocity $v_{3}$, computed at $t=0.14~{}\text{s}$. $\text{P}_{\text{f}}$ and $\text{P}_{\text{s}}$ represent the modes corresponding to fast and slow P waves, respectively. The point source is located at $\bm{x}_{0}=(500~{}\text{m},~{}300~{}\text{m})$, marked by a white star in Figure 11a. Figure 11b shows a comparison between the normalized solid vertical particle velocity $(v_{3})$ obtained from both the DG method and spectral-element method. The time history of the solution is recorded at $(600~{}\text{m},~{}400~{}\text{m})$, denoted by the red diamond in Figure 11a. The agreement between the solutions is good. The root mean square error (RMSE) between the solutions is $8.68\text{e}-03$. A comparison between the solution obtained from the DG method and a pseudo-spectral method [9] is also shown in Figure 11c. The comparison shows a good agreement between the solutions except at the onset time $0.06~{}\text{s}$. The RMS misfit between the solution is $3.18\text{e}-02$.

Next, we validate numerical solutions obtained from the DG method with a spectral-element method for a heterogenous medium. A two layer medium is constructed with material properties given in column 2 (top layer) and 3 (bottom layer) of Table 2. The size of the computational domain is $(4800~{}\text{m}\times 4800~{}\text{m})$ and discretized with the triangular and quadrilateral elements for DG and spectral-element method, respectively. The minimum size of the edge of the element is taken as $200~{}\text{m}$. Figure 12a represents the snapshot of solid particle velocity at $t=0.73~{}\text{s}$. The pure dilatational point source is located at $\bm{x}_{0}=(1600~{}\text{m},2900~{}\text{m})$, denoted by white star in Figure 12a, and implemented in same way as discussed for homogeneous medium. Figure 12b shows a comparison between solutions obtained from the DG and spectral-element methods. The time history of the solutions is recorded at $(2000~{}\text{m},1867~{}\text{m})$, marked by a red diamond in Figure 12a. The comparison shows a good overall agreement, though some discrepancies are present around time $T=1$. The RMS misfit between the solutions is $6.65\text{e}-02$.

5.4.6 3D models and computational results

Finally, we present numerical solutions of the poroelastic wave equation on 3D domains. First, we compute the solution in an homogeneous medium and then extend the computation to a realistic synthetic reservoir model constructed from a varying layer with an undulated topography.

5.5 Epoxy-glass model

First, we compute the solution for a 3D cube made of epoxy-glass with material properties given in column 3 of Table 1. The size of the computational domain is $2~{}\text{km}\times 2~{}\text{km}\times 2~{}\text{km}$ in $x$, $y$ and $z$ directions, respectively. The domain is discretized with tetrahedral elements with a minimum edge length of $45~{}\text{m}$. Figure 13(a)-(b) represents the $x-$ and $z-$ components of the solid particle velocity in the epoxy-glass medium at $0.25~{}\text{s}$. The central frequency of the forcing function is $f_{0}=35~{}\text{Hz}$, and polynomials of degree $N=5$ are used.

Three events can be observed: the fast P mode ($\text{P}_{\text{f}}$, outer wavefront), the shear wave (S, middle wavefront), and the slow P mode (inner wavefront). In an orthotropic material, the effect of anisotropy on wave propagation is observed in two orthogonal planes, whereas isotropic effects are observed only in one plane. Figure 13a clearly shows the effect of anisotropy on wave propagation with an isotropic effect in the plane parallel to the $y$ axis.

5.6 Reservoir model

We construct a 3D reservoir model characterized by rock layers, discontinuity, and a surface with undulated topography. The discretized model along with the density distribution is shown in Figures 14a and 14b, respectively. The dimension of the model is $22.8~{}\text{km}\times 17.4~{}\text{km}\times 7.0~{}\text{km}$ in $x$, $y$ and $z$ directions, respectively. The domain is discretized with tetrahedral elements with a minimum edge length of $125~{}\text{m}$. The top surface of the model is perturbed so that the effects of the topography, assumed as a free surface, could be incorporated into numerical simulations. Figure 15(a)-(b) represent the $x-$ and $z-$ components of the solid particle velocity at $3.5~{}\text{s}$. The central frequency of the forcing function is $f_{0}=10~{}\text{Hz}$, and polynomials of degree $N=3$ are used. The various modes of transmissions, reflections and scattering can be clearly seen in Figures 15a and 15b. Figure 16(a)-(b) show the $x-$ and $z-$ components of the fluid particle velocity at $3.5~{}\text{s}$.

6 Conclusions

This work presents a weight-adjusted discontinuous Galerkin (WADG) method for the linear poroelastic wave equations with arbitrary heterogeneous media. The method is energy stable and high order accurate for arbitrary stiffness tensors. The penalty numerical fluxes for this formulation are simple to derive and implement, and their lack of dependence on the stiffness tensors allows for a unified treatment of isotropic and anisotropic media with micro-heterogeneities.

We confirm the high-order accuracy of the numerical method using an analytic plane wave solution in a poroelastic media. Results obtained using this method also show good agreement with existing results in the literature for both problems involving both homogeneous and heterogeneous media. Finally, we provide computational results for a large 3D model consisting the surface topography.

We note that the implementation of this method reduces to the application of the weight-adjusted mass matrix inverse and the evaluation of constant-coefficient terms in the DG formulation. The cost of the latter step can be reduced (especially at high orders of approximation) by using fast methods based on Bernstein-Bezier polynomials for the application of derivative and lift matrices for constant-coefficient terms [38]. Future work will also involve an extension of the proposed formulation for the broad-band Biot’s model and coupled acousto-poroelastic media.

7 Acknowledgments

The authors gratefully thank sponsors of the Geo-Mathematical Imaging Group at Rice University for providing the resources to carry out this work. Jesse Chan gratefully acknowledges support from the NSF under grants DMS-1719818 and DMS-1712639. MVdH gratefully acknowledges support from the Simons Foundation under the MATH + X program and the NSF under grant DMS-1559587. The authors gratefully thank Ruichao Ye, Vitaly Katsnelson, Peter Caday, Christopher Wong and Sundeep Sharma for helpful and informative discussions. We also acknowledge Dr. Octavio Castillo from Barcelona Supercomputing Center (BSC), Barcelona, Spain for his help in 3D meshing of reservoir model.

Appendix A Matrix weight-adjusted inner products

Let $D^{\alpha}\bm{v}$ denote the component-wise differentiation of $\bm{v}$ with respect to a $d-$dimensional multi-index $\alpha$. Then, vector $L^{p}$ Soblev norms for $\bm{v(x)}\in\mathbb{R}$ can be defined as

[TABLE]

[TABLE]

where $W^{k,p}$ and $W^{k,\infty}$ are the Sobolev spaces [47]. Let $\bm{W}(\bm{x})$ be a matrix-valued weight function which is pointwise symmetric positive definite

[TABLE]

The $k^{th}$ order Sobolev norm for $\bm{W}(\bm{x})$ is defined as

[TABLE]

Appendix B Weight-adjusted approximations with matrix weights

Let $\Pi_{N}\bm{u}$ defined as $L^{2}$ projection applied to each component of the vector-valued function $\bm{u}$. We then define the operator $T_{\bm{W}}$ as

[TABLE]

The inverse operator $T_{W}^{-1}$ is defined implicitly via

[TABLE]

Lemma 1**.**

Let $\Pi_{N}$ denote the component-wise $L^{2}$ projection and $\bm{W}\in{(L^{\infty})}^{m\times m}$. Then $T_{\bm{W}}$ satisfies the following properties:

$T_{\bm{W}}^{-1}T_{\bm{W}}=\Pi_{N}$** 2. 2.

$\Pi_{N}T^{-1}_{\bm{W}}=\bm{T}_{\bm{W}}^{-1}=T_{\bm{W}}^{-1}$** 3. 3.

$||T_{\bm{W}}^{-1}||_{L^{2}(D^{k})}\leq\tilde{w}_{\max}.$** 4. 4.

$(T_{W}^{-1}\bm{v},\bm{w})_{L^{2}(D^{k})}$* forms an inner product on $(P^{N})^{m}\times(P^{N})^{m}$, which is equivalent to the $L^{2}$ inner product with equivalence constants $C_{1}=\tilde{w}_{\min},~{}C_{2}=\tilde{w}_{\max}$.*

The proof of Lemma 1 is given in detail in [28].

Appendix C Convergence analysis of WADG scheme for poroelastic wave equations

Using the results in Section 4.2, we can extend the semi-discrete convergence analysis in [48, 41, 49] in to linear poroelastic wave propagation on meshes of affine elements. Let $\bm{U}$ and $\bm{U}_{h}$ represent the exact and discrete WADG solutions, respectively. We assume that $\bm{U}$ and $\displaystyle\frac{\partial\bm{U}}{\partial t}$ are sufficiently regular such that

[TABLE]

with ${\left\|\bm{U}\right\|}^{2}_{W^{N+1,2}{(\Omega_{h})}}=\sum_{k}{\left\|\bm{U}\right\|}^{2}_{W^{N+1,2}{(D^{k})}}$.

The WADG formulation can be written as

[TABLE]

[TABLE]

[TABLE]

where $\bm{A}_{0}(\bm{x})=\text{diag}(\bm{Q}_{s},\bm{Q}_{v})$, $\bm{U}$ and $\bm{V}$ are group variables, defined as $\bm{U}=(\bm{\tau},\bm{v})$ and $\bm{V}=(\bm{h},\bm{g})\in(V_{h}(\Omega_{h}))^{d}\times(V_{h}(\Omega_{h}))^{N_{d}}$.

The DG formulation in (11) is consistent and thus

[TABLE]

$\forall~{}\bm{V}\in(V_{h}(\Omega_{h}))^{d}\times(V_{h}(\Omega_{h}))^{N_{d}}$. We decompose the error $\bm{U}-\bm{U}_{h}$ into a projection error $\bm{\epsilon}$ and a discretization error $\bm{\eta}$.

[TABLE]

We assume that $\bm{U}_{h}(\bm{x},0)$ is the $L^{2}$ projection of the exact initial condition, such that $\bm{\eta}|_{t=0}=0$. We also introduce a consistency error $\delta=\bm{A}_{0}\bm{U}-T^{-1}_{\bm{A}_{0}^{-1}}\bm{U}$ resulting from the approximation of $\bm{A}_{0}\bm{U}$ by a weight-adjusted inner product

[TABLE]

where we have used that $T^{-1}_{\bm{A}_{0}^{-1}}\Pi_{N}$. Subtracting the DG and WADG formulations in (34) and (35) and setting $\bm{V}=\bm{\eta}$ yields

[TABLE]

where we have used that $a(\bm{\eta},\bm{\eta})=0$ from skew-symmetry. We bound $a(\bm{\epsilon},\bm{\eta})+b(\bm{\epsilon},\bm{\eta})$ in (36) by integrating by parts the stress equation and using the component-wise $L^{2}$ orthogonality of $\bm{\epsilon}$ to derivatives of $\bm{\eta}$. This reduces to

[TABLE]

where $C_{\tau}$ is maximum of $(\alpha_{\bm{\tau}},\alpha_{\bm{v}})$ and we have also used the property of negative sime-definite matrix $\bm{D}$. Using the Young’s inequality with $\alpha=C_{\tau}/2$ yields the following bound

[TABLE]

Applying this to (36) and using Cauchy-Schwarz on $\left(\displaystyle\frac{\partial\delta}{\partial t},\bm{\eta}\right)_{L^{2}(\partial D^{k})}$ yields

[TABLE]

Using a $hp$ trace inequality [50]

[TABLE]

We now use that that $\bm{Q}_{s}$ and $\bm{Q}_{v}$ are positive definite such that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We define $A_{\max}=\max\left({s}_{\text{max}}\tilde{s}_{\text{max}},{v}_{\text{max}}\tilde{v}_{\text{max}}\right),~{}A_{\min}=\min\left(\tilde{s}_{\text{min}},\tilde{v}_{\text{min}}\right)$ and substitute the result of (38) in (36) performing time integration from $[0,T]$ and using modified Gronwall’s inequality (Lemma 1.10 in [51]), (38) reduces to

[TABLE]

Using the triangle inequality yields

[TABLE]

where $C_{1}$ and $C_{2}$ depends on $A_{\min}$ and $A_{\max}$. These estimates show that $L^{2}$ errors decrease proportionally to $h^{N+1/2}$. In practice we often observe the rate of convergence to be $O(h^{N+1})$.

Appendix D Analytical solution of diffusive part of the poroelastic system

From (9), the diffusive part of the system is expressed as

[TABLE]

The component wise expression for intermediate velocity vector is

[TABLE]

where $\beta_{v}^{(i)}=\displaystyle\frac{\rho_{f}}{\rho_{f}^{2}-\rho m_{i}}$ and $\beta_{q}^{(i)}=\displaystyle\frac{\rho}{\rho_{f}^{2}-\rho m_{i}}$. The above system can be solved analytically, giving

[TABLE]

where $\lambda_{s}^{(i)}=-(\eta/\kappa_{i})\beta_{q}^{(i)}$. The state of intermediate vector $[v_{1}^{*},~{}v_{2}^{*},~{}v_{3}^{*},~{}q_{1}^{*},~{}q_{2}^{*},~{}q_{3}^{*}]$ is the input for the non-stiff being solved with DG scheme.