On mathematical and numerical modelling of multiphysics wave propagation with polytopal Discontinuous Galerkin methods

In this work we review discontinuous Galerkin finite element methods on polytopal grids (PolydG) for the numerical simulation of multiphysics wave propagation phenomena in heterogeneous media. In particular, we address wave phenomena in elastic, poro-elastic, and poro-elasto-acoustic materials. Wave propagation is modeled by using either the elastodynamics equation in the elastic domain, the acoustics equations in the acoustic domain and the low-frequency Biot's equations in the poro-elastic one. The coupling between different models is realized by means of (physically consistent) transmission conditions, weakly imposed at the interface between the subdomains. For all models configuration, we introduce and analyse the PolydG semi-discrete formulation, which is then coupled with suitable time marching schemes. For the semi-discrete problem, we present the stability analysis and derive a-priori error estimates in a suitable energy norm. A wide set of verification tests with manufactured solutions are presented in order to validate the error analysis. Examples of physical interest are also shown to demonstrate the capability of the proposed methods.


Introduction
Multiphysics wave propagation in heterogeneous media is a very attractive research topic and, in recent decades, it has registered considerable interest in the mathematical, geophysical and engineering communities.Mathematical models for wave propagation phenomena range from the linear transport equation, to the non-linear system of Navier-Stokes equations.They appear in many different scientific disciplines, such as acoustic engineering [78], vibroacoustics [59], aeronautical engineering, [37], biomedical engineering [55], and computational geosciences; see [36] for a comprehensive review.
Thanks to the ongoing development of increasingly advanced high-performance computing facilities, the use of digital twins models for solving wave propagation problems has given a notable impulse towards a deeper understanding of these phenomena.Numerical methods designed for wave simulations must simultaneously account for the following three distinguishing features: accuracy, geometric flexibility and scalability.Accuracy is essential to correctly reproduce the physical phenomenon, and allows to minimize numerical dispersion and dissipation errors that would deteriorate the quality of the solution.Geometric flexibility is required since the computational domain usually features complicated geometrical shapes as well as sharp media contrasts.Scalability is demanded to solve on parallel machines real computational models featuring several hundred of millions or even billions of unknowns.
In this work we consider wave propagation problems arising from geophysics and we discuss and analyze several models, with increasing complexity, employed in this scientific area.We first present models of elastodynamics, then of poro-elasticity, and finally coupled poro-elasto-acoustics models.
Elastodynamics and viscoelastodynamics models are typically used for the study of seismic waves that propagate across the globe and are generated by earthquakes, volcanic activity, or artificial explosions.As far as the elastodynamics equations are concerned, the most used numerical methods are finite differences [38,62,67], finite elements [23], finite volumes [30,44,45,66], and spectral elements in either their conforming [51,56,77] or discontinuous setting [11,42,48].
Poro-elastodynamics models are used to describe the propagation of pressure and elastic waves through a porous medium.Pressure waves propagate through the saturating fluid inside pores, while elastic ones through the porous skeleton.In the pioneering work by Biot [25] general equations of waves propagation in poro-elastic materials were introduced.More recently, in [75] it is proposed a model of seismic waves in saturated soils, distinguishing between in-phase (fast) movements between solid and fluid and out-phase (slow) ones.Poro-elasto-acoustic problem model acoustic/sound waves impacting a porous material and consequently propagating through it.The coupling between acoustic and poro-elastic domains, realized by means of physically consistent transmission conditions at the interface, is discussed in [54] and [40].
The aim of this work is to introduce and analyze a discontinuous Galerkin method on polygonal/polyhedral (polytopal, for short) grids (PolydG) for the numerical discretization of multiphysics waves propagation through heterogeneous materials.The geometric flexibility and arbitrary-order accuracy featured by the proposed scheme are crucial within this context as they ensure the possibility of handling complex geometries and an intrinsic high-level of precision that are necessary to correctly represent the solutions.
For early results in the field of dG methods we refer the reader to [7,12,20,31,33,35,41] for secondorder elliptic problems, to [32] for parabolic differential equations, to [10] for flows in fractured porous media, and to [18] for fluid structure interaction problems.In the framework of dG methods for hyperbolic problems, we mention [53,72] for scalar wave equations on simplicial grids and the more recent PolydG discretizations designed in [15] for elastodynamics problems, in [16] for non-linear sound waves, in [4,5] for coupled elasto-acoustic problems, and in [6] for poro-elasto-acoustic wave propagation.
The results of the present work review and extend the analysis conducted in [15] and [6].In particular, in Section 5 we provide a novel stability estimate for the poro-elastic case requiring minimal regularity on problem data and showing explicitly the dependence on the model coefficients, final simulation time, and initial conditions.Additionally, in Section 6 we generalize the PolydG semi-discrete formulation derived in [6, Section 3] to the heterogeneous case, namely we allow all the model coefficients to be discontinuous over the computational domain.
The remaining part of the paper is structured as follows: in Section 2 we review the differential models for wave propagation in heterogeneous Earth's media while in Section 3 we define the discrete setting used in the paper.The elastodynamic model and its numerical discretization through a PolydG method is recalled in Section 4, while Sections 5 and 6 discuss the numerical analysis of a PolydG method for wave propagation problems in poro-elastic and coupled poro-elastic-acoustic media, respectively.Section 4, 5, and 6 contain at the end suitable verification test cases to validate the theoretical error bounds.Different numerical tests of physical interest are introduced and discussed in Section 7. Finally, in Section 8 we draw some conclusions and discuss some perspective about future work.

Notation
In the following, for an open, bounded domain  ⊂ R  ,  = 2, 3, the notation L 2 () is used in place of [ 2 ()]  , with  ∈ {2, 3}.The scalar product in  2 () is denoted by (•, •)  , with associated norm •  .Similarly, H ℓ () is defined as [ ℓ ()]  , with ℓ ≥ 0, equipped with the norm • ℓ, , assuming conventionally that H 0 () ≡ L 2 ().In addition, we will use H (div, ) to denote the space of L 2 () functions with square integrable divergence.In order to take into account essential boundary conditions, we also introduce the subspaces with Γ  ⊂  having strictly positive Hausdorff measure.Given  ∈ N and a Hilbert space H, the usual notation   ([0, ]; H) is adopted for the space of H-valued functions, -times continuously differentiable in [0, ].Finally, the notation   stands for  ≤  , with  > 0 independent of the discretization parameters, but possibly dependent on the physical coefficients and the final time .

Modelling seismic waves
A seismic event is the result of a sudden release of energy due to the rupture of a more fragile part of the Earth's crust called the fault.The deformation energy, accumulated for tens and sometimes hundreds of years along the fault, is transformed into kinetic energy that radiates, in the form of waves, in all directions through the layers of the Earth.Seismic waves are therefore energy waves that produce an oscillatory movement of the ground during their passage.Seismic waves are subdivided into two main categories: volume waves and surface waves.The former can be decomposed into compression waves (P) and shear waves (S).The (faster) P waves are transmitted both in liquids and in solids, while the (slower) S waves travel only in solid media.P waves induce a ground motion aligned with the wave field direction while S waves induce ground a motion in a plane perpendicular to the wave propagation field.
More and more frequently mathematical models are used for the study and analysis of ground motion.The solution of these models through appropriate numerical methods can provide important information for the evaluation of the seismic hazard of a given region and for the planning of the territory in order to limit the socio-economic losses linked to the seismic event.In the following we consider the differential model that aims at describing the propagation of seismic wave within Earth's interior.
Let Ω be a bounded domain modeling the portion of the Earth where the passage of seismic waves occurs, and let Ω be its boundary that can be decomposed into three disjoint parts Γ  , Γ  , and Γ  .The values of the displacement (Dirichlet conditions), the values of tractions (Neumann conditions), and the values of fictitious tractions (absorbing conditions) are imposed on Γ  , Γ  , and Γ  , respectively.For a temporal interval (0, ], with  > 0, the equation governing the displacement field u(x, ) of a dynamically disturbed elastic medium can be expressed as where  is the mass density, f defines a suitable seismic source and σ is the stress tensor that models the constitutive behaviour of the material.Possible definition for σ and f will be discussed in the sequel.Equation ( 4) is completed by prescribing suitable boundary conditions as well as initial conditions.For the latter, by choosing u(•, 0) =   u(•, 0) = 0, we suppose the domain to be at rest at the initial observation time.

Seismic waves in viscoelastic media
The stress tensor σ in (4) can be defined in different ways to properly model the behavior of the soil.Before presenting the main constitutive laws that can adopted for seismic wave propagation analysis we introduce: (i) the strain tensor , defined as the symmetric gradient (u) = (∇u + ∇  u)/2, and (ii) the fourth-order (symmetric and positive definite) stiffness tensor D, encoding the mechanical properties of the medium.It is expressed in term of the first and the second Lamé coefficients, namely  and , respectively.For an elastic material the generalized Hooke's law defines the most general linear relation among all the components of the stress and strain tensor.In the most general case, i.e. a fully anisotropic material, equation ( 5) contains 21 material parameters.However in our case, i.e., for a perfectly isotropic material, (5) can be reduced as where I is the identity tensor.
Pure elastic constitutive laws are not physically representative in the field of application of interest.A first model for visco-elastic media can be handled by modifying the equation of motion according to [58].In the approach, the inertial term   u in ( 4) is replaced by   u + 2   u +  2 u where  is an attenuation parameter.As a matter of fact, with this substitution, i.e., all frequency components are equally attenuated with distance, resulting in a frequency proportional quality factor  > 0 [65].A second attenuation model is obtained by considering materials "endowed with memory" in the sense that the state of stress at the instant  depends on all the deformations undergone by the material in previous times.This behaviour can be expressed through an integral equation of the form where the stress σ is determined by the entire strain history.Implicit in this law is the dependence on time of the Lamé parameters  and , cf.[57,63].We remark that, by using (8) it is possible to obtain an almost constant quality factor  in a suitable frequency range, cf.[63].
We conclude this section by addressing proper boundary conditions to supplement equation (7).Several conditions can be set to correctly define the interaction between the wave and the domain boundary.Dirichlet conditions are employed to prescribe the behaviour of the displacement field, i.e. u = g  on Γ  , while Neumann conditions σn = g  on Γ  represent the distribution of surface loads.Here n denotes the outward pointing normal unit vector with respect to Ω.
For geophysical applications, since the domain of interest Ω represents a portion of the Earth the following boundary conditions are commonly adopted: (i) free-surface condition, i.e. σn = 0 for the top Earth's surface and (ii) transparent boundary conditions σn = t for the remaining lateral and bottom surfaces.The latter consists in modeling the absorbing boundary layers by introducing a fictitious traction term t = t(u,   u), consisting of a linear combination of displacement space and time derivatives.Examples can be found in [21,68].In Figure 1 we report an illustrative example of domain Ω together with boundary conditions.

Seismic waves in porous media
Modeling wave propagation through fluid-saturated porous rock is crucial for the characterization of the seismic response of geologic formations.In this case, the effects stemming from the interaction between the viscous fluid and the solid matrix have to be taken into account.In the framework of Biot's poroelasticity theory [24,25], the total stress tensor σ additionally depends on the pore pressure  according to the following relation with σ(u) defined as in (6) and 0 <  ≤ 1 denoting the Biot coefficient.Adding to the momentum balance equation ( 4) the inertial term corresponding to the filtration displacement w = (w  − u), where  > 0 is the reference porosity and w  the fluid displacement, leads to Here, the average density  is given by  =   + (1 − )   , where   > 0 is the saturating fluid density and   > 0 is the solid density.To derive Biot's wave equations in Section 5, the rheology of the porous material (9) and the momentum balance (10) are combined with the dynamics of the fluid system described by Darcy's law and the conservation of fluid mass in the pores.Two major differences have been observed when dealing with poro-elastic media instead of elastic ones: (i) the attenuation due to wave-induced fluid flow and (ii) the presence of an additional compressional wave of the second kind (slow P-wave), which becomes a diffusive mode in the low-frequency range, cf.[36].As observed in [43], this slow P-wave is mainly localized near the material heterogeneities or the source.

Modelling the seismic source
Seismic wave can be generated by different natural and artificial sources.Depending on the problem's configuration, one can consider a single point-source, an incidence plane wave or a finite-size rupturing fault.
We can define a point-wise force f acting on a point x 0 ∈ Ω in the   ℎ direction as where e  is the unit vector of the   ℎ Cartesian axis, (•) is the delta distribution, and  (•) is a function of time.The expression of  (•) can be selected among different waveforms.Here, we report one of the most employed one, i.e. the Ricker wavelet [70], defined as where  0 is the wave amplitude,   is the peak frequency of the signal and  0 is a fixed reference time.
To define a vertically incident plane wave one can consider a uniform distribution of body forces along the plane  =  0 of the form f (x, ) =  ()e  ( −  0 ).The latter generates a displacement in the   ℎ direction given by ū where  (•) is the Heaviside function and  (that can be equal either to   = √︁  + 2/ or   = √︁ /) is the wave velocity, see [52].By taking the derivative of (13) with respect to time and evaluating the result at  =  0 we can express  () as  () = 2  ū  .Finally, we introduce one of the most important seismic input for seismic wave propagation that is the so called double-couple source force.A point double-couple or moment-tensor source localized in the computational domain Ω is often adopted to simulate small local or near-regional earthquakes.Its mathematical representation is based on the seismic moment tensor m(x, ), defined in [1] as where n  and s  denote the fault normal and the rake vector along the fault, respectively. 0 (x, ) describes the time history of the moment release at x and  is the force elementary volume.The equivalent body force distribution is finally obtained through the relation f (x, ) = −∇ • m(x, ), see [47].

Discrete setting for PolyDG methods
In this section we define the notation related to the subdivision of the computational domain Ω by means of polytopic meshes.We introduce a polytopic mesh T ℎ made of general polygons (in 2d) or polyhedra (in 3d).We denote such polytopic elements by , define by |  | their measure and by ℎ  their diameter, and set ℎ = max  ∈ T ℎ ℎ  .We let a polynomial degree   ≥ 1 be associated with each element  ∈ T ℎ and we denote by  ℎ : The discrete space is introduced as follows: , where P  ℎ (T ℎ ) is the space of piecewise polynomials in Ω of total degree less than or equal to   in any  ∈ T ℎ .
In order to deal with polygonal and polyhedral elements, we define an interface of T ℎ as the intersection of the ( − 1)-dimensional faces of any two neighboring elements of T ℎ .If  = 2, an interface/face is a line segment and the set of all interfaces/faces is denoted by F ℎ .When  = 3, an interface can be a general polygon that we assume could be further decomposed into a set of planar triangles collected in the set F ℎ .We decompose the faces of T ℎ into the union of internal () and boundary () faces, respectively, i.e.: F ℎ = F  ℎ ∪ F  ℎ .Moreover we further split the boundary faces as , meaning that on F  ℎ (resp.F  ℎ and F  ℎ ) Dirichlet (resp.Neumann and absorbing) boundary conditions are applied.Following [34], we next introduce the main assumption on T ℎ .Definition 1.A mesh T ℎ is said to be polytopic-regular if for any  ∈ T ℎ , there exists a set of nonoverlapping -dimensional simplices contained in , denoted by {   }  ⊂ , such that for any  ⊂ , the following condition holds: Assumption 1.The sequence of meshes {T ℎ } ℎ is assumed to be uniformly polytopic regular in the sense of Definition 1.
We remark that this assumption does not impose any restriction on either the number of faces per element nor their measure relative to the diameter of the element they belong to.Under Assumption 1, the following trace-inverse inequality holds: see [34,Section 3.2] for the detailed proof and a complete discussion on inverse estimates.In order to avoid technicalities, we also make the following ℎ-local bounded variation property assumption.

Assumption 2. For any pair of neighboring elements 𝜅
Next, following [19], for sufficiently piecewise smooth scalar-, vector-and tensor-valued fields , v and τ , respectively, we define the averages and jumps on each interior face  ∈ F  ℎ shared by the elements  ± ∈ T ℎ as follows: where ⊗ is the tensor product in R 3 , • ± denotes the trace on  taken within the interior of  ± , and n ± is the outward unit normal vector to  ± .Accordingly, on boundary faces  ∈ F  ℎ , we set Finally, we introduced some important concepts employed for the convergence analysis of PolydG methods presented in the sequel, namely, the mesh covering T ♯ and the Stein extension operator E. Indeed, the latter are used to extend standard ℎ-interpolation estimates on simplices to polytopal elements.We refer the reader to [8,32,34,35] for all the details.
A covering T ♯ = {K  } related to the polytopic mesh T ℎ is a set of shape regular −dimensional simplices K  such that for each  ∈ T ℎ there exists a K  ∈ T ♯ such that  ⊂ K  .We suppose that there exits a covering T ♯ of T ℎ and a positive constant  Ω , independent of the mesh parameters, such that max and ℎ   ℎ  for each pair  ∈ T  and   ∈ T ♯ with  ⊂ T ℎ .This latter assumption assures that, when the computational mesh T ℎ is refined, the amount of overlap present in the covering T ♯ remains bounded.
For an open bounded domain Σ ⊂ R  and a polytopic mesh T ℎ over Σ satisfying Assumption 1, we can introduce the Stein extension operator E :   () →   (R  ) [76], for any  ∈ T ℎ and  ∈ N 0 , such that Ẽ |  =  and Ẽ ,R   , .The corresponding vector-valued version mapping H  () onto H  (R  ) acts component-wise and is denoted in the same way.In what follows, for any  ∈ T ℎ , we will denote by K  the simplex belonging to T ♯ such that  ⊂ K  .

Time integration
We introduce here the time integration scheme used for the numerical simulations shown in the following sections.First, we anticipate that, fixing a suitable basis for the discrete space, all the semi-discrete PolydG formulations that we will introduce in the following can be written in the general abstract form: where the precise definition of the unknown  ℎ , the right-hand side  ℎ , and the matrices M ℎ , D ℎ , and A ℎ will be given in the forthcoming sections.Assuming that M ℎ is invertible, we have Then, we discretize the interval [0, ] by introducing a timestep Δ > 0, such that ∀  ∈ N,  +1 −   = Δ and define   ℎ =  ℎ (  ) and   ℎ =  ℎ (  ).Finally, to integrate in time (16) we apply the Newmark− scheme defined by introducing a Taylor expansion for  ℎ and  ℎ =  ℎ , respectively: being L  ℎ = L ℎ (  ,   ℎ ,   ℎ ) and the Newmark parameters   and   satisfy the constraints 0 ≤   ≤ 1, 0 ≤ 2  ≤ 1.The typical choices are   = 1/2 and   = 1/4, for which the scheme is unconditionally stable and second order accurate.We also remark that, when L ℎ = L ℎ (  ,   ℎ ),   = 0, and   = 1/2, the Newmark scheme (18) reduces to the leap-frog scheme which is explicit and second order accurate.We next address in detail the PolydG semi-discrete approximation of the problems we are considering.

Elastic wave propagation in heterogeneous media
Hereafter, for the sake of presentation, we will consider the linear visco elastodyamics model, i.e. equations ( 7) and (6).We suppose Ω = Γ  ∪ Γ  and we consider homogeneous Dirichlet and Neumann boundary conditions on Γ  and Γ  , respectively.The system of equations can be recast as The case with non homogenous Neumann conditions is treated in [3], while absorbing conditions are considered in [68].Finally, we refer to [11,71] for a detailed analysis of viscoelastic attenuation models.We suppose the mass density  and the Lamé parameters  and  to be strictly positive bounded functions of the space variable x, i.e. , ,  ∈  ∞ (Ω).We also suppose the forcing term f to be regular enough, i.e., f ∈  2 ((0, ]; L 2 (Ω)) and that the initial conditions (u 0 , u 1 ) ∈ H 1 0 (Ω) × L 2 (Ω).The weak formulation of problem (19) reads as follows: for all  ∈ (0, ] find u = u() ∈ H 1 0 (Ω) such that where for any u, v ∈ H 1 0 (Ω) we have set Problem ( 20) is well-posed and its unique solution u ∈  ((0, ];

Semi-discrete formulation
Using the notation introduced in Section 3, we define the PolyDG semi-discretization of problem (20): for all  ∈ (0, ], find u ℎ = u ℎ () ∈ V ℎ such that for any v ℎ ∈ V ℎ , supplemented with the initial conditions (u ℎ (0),   u ℎ (0)) = (u 0 ℎ , u 1 ℎ ), where u 0 ℎ , u 1 ℎ ∈ V ℎ are suitable approximations of u 0 and u 1 , respectively.Here, we also assume the tensor D and the density  to be element-wise constant over T ℎ .The bilinear form A  ℎ : for all u, v ∈ V ℎ .Here, we adopt the compact notation 23) is defined face-wise as where for any  ∈ T ℎ (here | • | 2 is the operator norm induced by the  2 -norm on R  , where  denotes the dimension of the space of symmetric second-order tensors, i.e.,  = 3 if  = 2,  = 6 if  = 3), and  0 is a (large enough) positive parameter at our disposal.
By fixing a basis for V ℎ and denoting by U ℎ the vector of the expansion coefficients in the chosen basis of the unknown u ℎ , the semi-discrete formulation (22) can be written equivalently as: with M denoting the mass matrix in V ℎ , A  the stiffness matrix corresponding to the bilinear form A  , and with initial conditions U ℎ (0) = U 0 and U ℎ (0) = U 1 .Note that F ℎ is the vector representations of the linear functional (f , v ℎ ) Ω .Formulation ( 25) can be recast in the form ( 16)-( 17) by setting

Stability and convergence results
In this section we recall the stability and convergence results for the semidiscrete PolyDG formulation (22).We refer the reader to [15] and to [9] for all the details.The results are obtained in the following energy norm where with be the approximate solution of (22) obtained with the stability constant  0 defined in (24) chosen sufficiently large.Then, where u ℎ (0 DG, , being u 0,ℎ , u 1,ℎ ∈ V ℎ suitable approximation of the initial conditions u 0 and u 1 , respectively.The proof of the previous stability estimate can be found for instance in [9,15].From (28) it is possible to conclude that the PolyDG approximation is dissipative.Indeed, when f = 0 (no external forces) the energy of the system at rest u 0 ℎ E is not conserved through time.
Concerning the convergence results of the PolyDG scheme we report in the following the main result.We refer the reader to [15] for the details and for the proof of the following theorem.Theorem 3. Let Assumption 1 and Assumption 2 be satisfied and assume that the exact solution u of (20) is sufficiently regular.For any time  ∈ [0, ], let u ℎ ∈ V ℎ be the PolyDG solution of problem (22) obtained with a penalty parameter  0 appearing in (24) sufficiently large.Then, for any time  ∈ (0, ] the following bound holds where with   = min(   + 1,   ) for all  ∈ T ℎ .The hidden constant depends on the material parameters and the shape-regularity of the covering T ♯ , but is independent of ℎ  ,   .

Verification test
We solve the wave propagation problem (19) in Ω = (0, 1) 2 , choosing  =  =  =  = 1 and assuming that the exact solution u is given by Dirichlet boundary conditions and initial conditions are set accordingly.We set the final time  = 1 and chose a time step Δ = 10 −4 of the leap-frog scheme, cf.(18).The penalty parameter  0 appearing in (24) has been set equal to 10.We compute the discretization error by varying the polynomial degree   = , for any  ∈ T ℎ , and the number of polygonal elements   .
In Figure 3 (left), we report the computed energy error u − u ℎ  at final time  as a function of the mesh size ℎ.We retrieve the algebraic convergence proved in (29) for a polynomial degree  = 2, 3, 4. Next, we report the computed L 2 -error  u L 2 (Ω) = u − u ℎ L 2 (Ω) at time  obtained on a shape-regular polygonal grid (cf. Figure 2) versus the polynomial degree , which varies from 1 to 5, in semilogarithmic scale.We fix the number of polygonal elements as   = 160.In this case we observe an exponential converge in , as shown in Figure 3

Poro-elastic media
In this section, we consider a poro-elastic material occupying a polyhedral domain Ω  ⊂ Ω modeled by equations ( 9) and (10).The low-frequency Biot's system [25] can be written as in Ω  × {0}, (w,   w) = (w 0 , w 1 ) in where the density   is given by   =  −1   with tortuosity  > 1,  represents the dynamic viscosity of the fluid,  is the absolute permeability, and  denotes the Biot modulus.As in the previous section, we assume that the model coefficients   ,   ,  −1 ,  ∈  ∞ (Ω  ) are strictly positive scalar fields and that the source term f , g and the initial conditions (w 0 , w 1 ) are regular vector fields, namely f , g ∈  2 ((0, ]; L 2 (Ω  )) and (w 0 , w 1 ) ∈ H 0 (div, Ω  ) × L 2 (Ω  ).The third and fourth equations in (31) correspond to the dynamic Darcy's law and the conservation of fluid mass, respectively.For the sake of simplicity, in (31) we have also assumed that the clamped region Γ  ⊂ Ω  is impermeable and a null pore pressure condition is prescribed on the Neumann boundary We remark that more general boundary conditions can be treated up to minor modifications.
In what follows, we focus on the two-displacement formulation of the low frequency poro-elasticity problem [61], that is obtained by inserting the expression of the total stress σ and the pore pressure  in the other equations in (31).The corresponding weak formulation reads: for all  ∈ (0, ] find with A  : H 1 0 (Ω  ) × H 1 0 (Ω  ) → R defined as the restriction to Ω  of the function in (21) and the bilinear forms M  , A  defined as for all (u, w), (v, z) ∈ H 1 0 (Ω  ) ×H 0 (div, Ω  ).The well-posedness of the low-frequency poro-elasticity problem (32) has been established in [46,Section 5.2] in the framework of semigroup theory.

Semi-discrete formulation
Proceeding as in Section 4.1, we derive the semi-discrete PolydG approximation of problem (32).We introduce a polytopic mesh T  ℎ of Ω  satisfying Assumptions 1 and 2 and denote by F  ℎ the set of faces of T  ℎ .Here, we consider the same polynomial space for both the discrete solid displacement u ℎ and filtration displacement w ℎ , i.e. u ℎ , w ℎ ∈ V  ℎ = (P  ℎ (T  ℎ ))  , and we assume that all the model coefficients are piecewise constant over T  ℎ .The PolydG semi-discrete problem consists in finding, for all  ∈ (0, ], the solution (u ℎ (), 23) and the bilinear form A  ℎ defined such that for all w, z ∈ V  ℎ and the penalization function  ∈  ∞ (F  ℎ ) is given by where  0 is a positive user-dependent parameter.We remark that, owing to the H (div)-regularity of the filtration displacement w solving (32), the penalization term in (35) acts only on the normal component of the jumps.Problem (34) is completed with suitable initial conditions (u ℎ (0), w ℎ (0),   u ℎ (0), We conclude this section by observing that the algebraic representation of the semi-discrete formulation (34) is given by with [ ℎ ,  ℎ ,  ℎ ,  ℎ ] (0) = [ 0 ,  0 ,  1 ,  1 ] and [ ℎ ,  ℎ ] T corresponding to the vector representation of the right-hand side of (34).Recalling the notation introduced in Section 3.1 and setting can be rewritten in the form (16).

Stability and convergence results
The aim of this section is to establish an a priori estimate for the solution of problem (34).First, we define for all u, w ∈  1 ( [0, ]; V  ℎ ) the energy function with   =   (1−)
Remark 5. We observe that, proceeding as in [26, Lemma 7], it is possible to obtain a stability estimate for problem (34) requiring The key step is based on estimating the term ∫  0 (f ,   u ℎ ) Ω  by using partial integration and the discrete Korn's first inequality [29, Lemma 1].
For the sake of conciseness, we decide not to present here the convergence analysis for the PolydG formulation of the poro-elastic problem (34).However, an error estimate can be readily deduced from Theorem 59 below, in the case in which the exact solution on the acoustic part of the domain is null.

Verification test
We consider problem (31) in Ω  = (−1, 0) × (0, 1) and choose as exact solution As before, Dirichlet boundary conditions and initial conditions are set accordingly.The model problem is solved on a sequence of polygonal meshes as the one shown in Figure 5 (left), with physical parameters shown in Figure 5 (right).The final time  has been set equal to 0.25, considering a timestep of Δ = 10 −4 for the Newmark- scheme,   = 1/2 and   = 1/4, cf.(18).The penalty parameters  0 and  0 appearing in definitions ( 24) and ( 36), respectively, have been chosen equal to 10.
In Figure 6 (left) we report the computed energy error (u − u ℎ , w − w ℎ ) E , cf. ( 59), as a function of the mesh size ℎ for a polynomial degree  = 2, 3, 4. In this case we retrieve the rate of convergence O (ℎ  ) as proved in (59).In Figure 6 (right) we plot the computed L 2 -errors for the elastic u and filtration w displacements as a function of the polynomial degree  in a semilog-scale.We fix the number of polygonal elements as   = 100.We observe an exponential rate of convergence since the solution ( 44) is analytic.

Poro-elastic-acoustic media
In this section, we present the PolydG discretization of the poro-elasto-acoustic interface problem.We refer the reader to [6] for the rigorous mathematical analysis of the model problem and the detailed derivation of the proposed method.In what follows, we assume that Ω is decomposed into two disjoint, polygonal/polyhedral subdomains: Ω = Ω  ∪ Ω  , cf. Figure 7.
The two subdomains share part of their boundary, resulting in the interface Γ  = Ω  ∩ Ω  .We set where the surface measures of Γ  , Γ  , and Γ  are assumed to be strictly positive.The outer unit normal vectors to Ω  and Ω  are denoted by n  and n  , respectively, so that n  = −n  on Γ  .
The subdomain Ω  represents a poro-elasto medium whose dynamical behavior is described by Biot's equations (31).In the fluid domain Ω  , we consider an acoustic wave with constant velocity  > 0 and mass density   > 0 such that   ,  −2 ∈  ∞ (Ω  ).For a given source term ℎ ∈  2 ((0, ];  2 (Ω  )), the acoustic potential  satisfies with ( 0 ,  1 ) ∈  1 0 (Ω  ) ×  2 (Ω  ).To close the coupled poro-elasto-acoustic problem, some interface conditions on Γ  are needed.Here, we consider physically consistent transmission conditions (see, e.g., [54] and [39]) expressing the continuity of normal stresses, continuity of pressure, and conservation of mass: The parameter  : Γ  → [0, 1] denotes the hydraulic permeability at the interface and models different pores configurations, cf. Figure 7.In the open pores region  −1 (1) ⊂ Γ  the second equation in (46) reduces to  =   , while in the sealed pores subset  −1 (0) we have w • n  = 0, implying that  −1 (0) is impermeable.Finally, the imperfect pores region  −1 ((0, 1)) models an intermediate state between open and sealed pores.For later use, we split the interface into two disjoint (possibly non-connected) subsets We remark that the first and second conditions in (46) plays the role of a Neumann and a Robin-like conditions for system (31), respectively.Similarly, the third equation in (46) acts as a Neumann condition for problem (45).The existence and uniqueness of a strong solution to the poro-elasto-acoustic problem coupling equations ( 31), (45), and ( 46) is proved in [6, Appendix A].

Semi-discrete formulation
We decompose the polytopic regular mesh T ℎ as T ℎ = T  ℎ ∪ T  ℎ , where T  ℎ and T  ℎ are aligned with Ω  and Ω  , respectively.In a similar way, we decompose F ℎ as , and F  ℎ and F  ℎ denote the faces of T  ℎ and T  ℎ , respectively, not laying on Γ  .The discrete spaces are selected as follows: given element-wise constant polynomial degrees  ℎ : T  ℎ → N * and  ℎ : . Finally, we also assume that the coefficients   and  are piecewise constant over T  ℎ and  is piecewise constant over F  ℎ .Under this assumption, we can decompose the set of mesh faces belonging to Γ  as The semi-discrete PolydG formulation of problem (48) consists in finding, for all  ∈ (0, ], the discrete solution As initial conditions we take the  2 -orthogonal projections onto (V  ℎ × V  ℎ ×   ℎ ) 2 of the initial data (u 0 , w 0 ,  0 , u 1 , w 1 ,  1 ).For all u, v, w, z ∈ V  ℎ and ,  ∈   ℎ , the bilinear forms A ℎ and C ℎ appearing in (50) are given by with A  ℎ : V  ℎ × V  ℎ → R defined as in (23) and Notice that the bilinear form A  ℎ is different from A  ℎ defined in (35).Indeed, the definition of A  ℎ in (53) also takes into account the essential condition z • n  = 0 on Γ   embedded in the definition of the functional space W  .The stabilization function  ∈  ∞ (F  ℎ ) is defined such that with  0 > 0 being a user-dependent parameter.

Stability and convergence results
In this section, we present the main stability and convergence results proved in [6].First, we introduce the energy norm defined such that, for all (u, w, ) with • E defined in (38) and • DG, : ∇ 2 The stability of the semi-discrete PolydG problem ( 50) is a consequence of Proposition 6 below, which also implies that the formulation is dissipative.Indeed, in the case of null external source terms, it follows from estimate ( 58) that (u ℎ , w ℎ ,  ℎ ) () E (u ℎ , w ℎ ,  ℎ ) (0) E for any  > 0. The proof of the following result is based on taking 50), using the skew-symmetry of the coupling terms, and then reasoning as in Proposition 4 (see [6,Theorem 3.4] for the details).Proposition 6.For sufficiently large penalty parameters  0 ,  0 ,  0 and for any  ∈ (0, ], the solution with hidden constant depending on time  and on the material properties, but independent of the interface parameter . In what follows, we report the main result concerning the error analysis of the PolydG discretization (50).To infer the error estimate of Theorem 7 below, an additional assumption on the interface permeability  is required.
ℎ  , with hidden constant independent of .
We remark that the previous assumption is used only for establishing the error estimate below but, according to our observation, it is not needed in practical applications.We refer the reader to [6,Theorem 4.3] for the detailed proof of the following result.large penalization parameters  0 ,  0 and  0 .Then, for any time  ∈ (0, ], the discretization error where , with   = min(   + 1,   ) and   = min(  + 1,   ) for all  ∈ T ℎ .The hidden constant depends on time , the material properties, and the shape-regularity of the covering T ♯ , but is independent of the discretization parameters and of .

Examples of physical interest 7.1 Two layered media
In this section we consider a wave propagation problem in heterogeneous media taken from [64].The aim of this test is to show how different assumptions on the model can determine and change the behavior of the wave propagation.
The domain of interest is Ω = (0, 4.8) 2 km 2 and consists of two layers as depicted in Figure 10.In the first case (a) the layers are perfectly elastic, cf.Table 1, while in the second case (b) the layers are assumed to be poro-elastic, cf.Table 2.A point-wise source f , cf. (11), acting in the − direction is located in the upper part of the domain at point x = (2.4,2.7) km.The time evolution of the latter is given by a Ricker-wavelet (12) with amplitude  0 = 1 m, time-shift  0 = 0.3 s and peak-frequency   = 5 Hz.For both models (a) and (b) we use a polygonal mesh with characteristic size ℎ = 10 −2 and a polynomial degree  = 3.We set homogeneous Dirichlet conditions on the boundary and use null initial conditions.To integrate in time model (a) we chose the leap-frog scheme while for model (b) the Newmark- scheme with parameters   and   as in the previous section.We fix the final time  = 1 s and chose Δ = 10 In Figure 11 we report selected snapshots of the computed magnitude of the velocity field |   u ℎ () | for models (a) and (b).As expected, the propagation of the wave in the elastic domain is regular and refraction phenomena are not very evident (due to a low contrast between the wave speeds).On the contrary, when porous media are accounted for, the refraction effects are more pronounced.This is in agreement with the findings in [64].

Wave propagation in layered poro-elastic-acoustic media
As a final test cases we consider the domain reproduced in Figure 12 where an acoustic layer is in contact with a heterogeneous poro-elastic body.
For the acoustic domain we set   = 1500 [kg/m 3 ] and  = 1000 [m/s].Physical parameters for the poro-elastic domain are chosen as in Table 2 where, for this case, the property of the former "Lower Layer" are assigned to the first poro-elastic subdomain, while those of the former "Upper Layer" to the second poro-elastic subdomain, cf.

Conclusions
In this work we have presented a review of the development of PolyDG methods for multiphysics wave propagation phenomena in elastic, poro-elastic and poro-elasto-acoustic media.
After having recalled the theoretical background of the analysis of PolyDG methods we analysed the well-posedness and stability of different numerical formulations and proved ℎ-version a priori error estimates for the semi-discrete scheme.Time integration of the latter is obtained based on employing the leap-frog or the a Newmark methods.Numerical experiments have been designed not only to verify the theoretical error bounds but also to demonstrate the flexibility in the process of mesh design offered by polytopic elements.In this respect, numerical tests of physical interest have been also discussed.
To conclude, PolyDG methods allow a robust and flexible numerical discretization that can be successfully applied to wave propagation problems.Future developments in this direction include the study of multi-physics problems such as fluid-structure (with poro-elastic or thermo-elastic structure) interaction problems (we refer, e.g., to [18,79] for preliminary results) as well as the exploitation of algorithms to design agglomeration-based multigrid methods and preconditioners for the efficient iterative solution of the (linear) system of equations stemming from PolyDG discretizations (see [13,14,17,27,28] for seminal results).

Figure 1 :
Figure 1: Example of domain Ω with boundary Ω divided into a Dirichlet Γ  , a Neumann Γ  and an absorbing Γ  part.

Figure 4 :
Figure 4: Example of a porous domain Ω  together with mixed boundary conditions on Γ  and Γ   .

Figure 6 :
Figure 6: Test case of Section 5.3.Computed energy error as a function of the mesh size ℎ for polynomial degree  = 2, 3, 4. The rate of convergence is also reported in the last row, cf.(59) (left).Computed L 2 -errors  u L 2 (Ω  ) = u − u ℎ L 2 (Ω  ) and  w L 2 (Ω  ) = w − w ℎ L 2 (Ω  ) as a function of the polynomial degree  in a semilogarithmic scale for   = 100 polygonal elements (right).

Figure 10 :
Figure 10: Test case of Section 7.1.Computational domain: the location of the point-source force is superimposed in red.

Figure 12 .
In this numerical example we chose the dynamic viscosity  equal to 0.001.Boundary and initial conditions have been set equal to zero both for the poroelastic and the acoustic domain.Forcing terms are null in Ω  , while in Ω  we consider a force of the form ℎ =  (, )(), where  is a Ricker wavelet of the form(12) with  0 = 1 [Hz m 3 ],   = 39.4784[Hz 2 ] and  0 = 0.75 s.The function  (, ) is defined as  (, ) = 1, if (, ) ∈ 4=1 (x  , ), while  (, ) = 0, otherwise, where (x  , ) is the circle centered in x  and with radius .Here, we set x 1 = (13097, 8868) m, x 2 = (16673, 8868) m, x 3 = (27079, 8868) m, x 4 = (29324, 8868) m and  = 100 m.Notice that, the support of the function  (, ) has been reported in Figure12, superimposed with a sample of the computational mesh employed.Simulations have been carried out by considering: a mesh consisting in  = 6356 triangles, subdivided into   = 2380 and   = 3976 triangles for the acoustic and poroelastic domain, respectively; a Newmark scheme with time step Δ = 10 −2 s and   = 1/2 and   = 1/4 in a time interval [0, 4] s; a polynomial degree   =   =  = 4.In Figure13, we show the computed pressure  ℎ considering the interface permeability  = 1.The latter value models an open pores condition at the interface, cf.(46).Remark that  ℎ =    ℎ in the acoustic domain while  ℎ = −(∇ • u ℎ + ∇ • w ℎ ) in the poro-elastic one.As one can see, the pressure wave correctly propagates from the acoustic domain to the poro-elastic one: the continuity at the interface boundary can be appreciated.Finally, we note how the second porous layer (sound absorbing material) produces a damping of the pressure field.

Figure 12 :
Figure 12: Test case of Section 7.2.Computational domain.Location of the acoustic sources are also superimposed.

Funding.
This work has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no.896616 (project PDGeoFF: Polyhedral Discretisation Methods for Geomechanical Simulation of Faults and Fractures in Poroelastic Media).The authors are members of the INdAM Research Group GNCS and this work is partially funded by INdAM-GNCS.Paola F. Antonietti has been partially funded by the research project PRIN n. 201744KLJL funded by MIUR.Conflict of interest.The authors have no conflicts of interest to declare that are relevant to the content of this article.Availability of data and material.The datasets generated during the current study are available from Ilario Mazzieri upon reasonable request.Authors' contributions.All authors contributed to the study conception and design.The first draft of the manuscript was written by Ilario Mazzieri and Michele Botti.All authors commented on previous versions of the manuscript and approved the final one.Ethics approval.Not applicable.

Table 1 :
Test case of Section 7.1.Physical parameters for the elastic medium.

Table 2 :
Test case of Section 7.1.Physical parameters for the poro-elastic medium.