Guillaume Chiavassa
email: guillaume.chiavassa@centrale-marseille.fr

Bruno Lombard
email: lombard@lma.cnrs-mrs.fr

Time domain numerical modeling of wave propagation in 2D heterogeneous porous media

Keywords: porous media, elastic waves, Biot's model, time splitting, finite difference methods, Cartesian grid, immersed interface method 2000 MSC: 35L50, 65M06 PACS: 43.20.-Gp, 46.40.-f

This paper deals with the numerical modeling of wave propagation in porous media described by Biot's theory. The viscous efforts between the fluid and the elastic skeleton are assumed to be a linear function of the relative velocity, which is valid in the low-frequency range. The coexistence of propagating fast compressional wave and shear wave, and of a diffusive slow compressional wave, makes numerical modeling tricky. To avoid restrictions on the time step, the Biot's system is splitted into two parts: the propagative part is discretized by a fourth-order ADER scheme, while the diffusive part is solved analytically. Near the material interfaces, a space-time mesh refinement is implemented to capture the small spatial scales related to the slow compressional wave. The jump conditions along the interfaces are discretized by an immersed interface method. Numerical experiments and comparisons with exact solutions confirm the accuracy of the numerical modeling. The efficiency of the approach is illustrated by simulations of multiple scattering.

Introduction

The propagation of waves in porous media has crucial implications in many areas, such as the characterization of industrial foams, spongious bones and petroleum rocks. The most widely used model describing the propagation of mechanical waves in a saturated porous medium was proposed by Biot in 1956. A major achievement in Biot's theory was the prediction of a second (slow) compressional wave, besides the (fast) compressional wave and the shear wave classically propagated in elastic media.

Two regimes are distinguished, depending on the frequency of the waves. At frequencies smaller than a critical frequency f c , the fluid flow inside the pores is of Poiseuille type, and the viscous efforts between the fluid and the solid depend linearly on the relative velocity. In this case, the slow compressional wave is almost static and highly attenuated [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I: Low-frequency range[END_REF]. An adequate modeling of this diffusive mode remains a major challenge in real applications. At frequencies greater than f c , inertial effects begin to dominate the shear forces, resulting in an ideal flow profile except in the viscous boundary layer, and the slow wave propagates [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. II: High-frequency range[END_REF][START_REF] Masson | Finite-difference modeling of Biot poroelastic equations at seismic frequencies[END_REF]. Experimental observations of the slow wave in the low-frequency range [START_REF] Plona | Observation of a second bulk compressional wave in a porous medium at ultrasonic frequencies[END_REF] and in the high-frequency range [START_REF] Chandler | The equivalence of quasi-static flow in fluid-saturated porous media and Biot's slow wave in the limit of zero frequency[END_REF] have confirmed Biot's theory. In the current study, we focus on the low-frequency range.

Until the 1990's, Biot's equations were mainly studied in the harmonic regime. Various time-domain methods have been proposed since, based on finite differences [START_REF] Dai | Wave propagation in heterogeneous porous media: a velocity-stress, finite-difference method[END_REF][START_REF] Zeng | The application of the perfectly matched layer in numerical modeling of wave propagation in poroelastic media[END_REF][START_REF] Wenzlau | Finite-difference modeling of wave propagation and diffusion in poroelastic media[END_REF], finite elements [START_REF] Zhao | An explicit finite element method for Biot dynamic formulation in fluid-saturated porous media and its application to a rigid foundation[END_REF], discontinuous Galerkin methods [START_REF] De La Puente | Discontinuous Galerkin methods for wave propagation in poroelastic media[END_REF], boundary elements [START_REF] Schanz | Application of 3-D time domain boundary element formulation to wave propagation in poroelastic solids[END_REF], pseudospectral methods [START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF] and spectral element methods [START_REF] Morency | Spectral-element simulations of wave propagation in porous media[END_REF]. A recent review of computational poroelasticity can be found in [START_REF] Carcione | Computational poroelasticity -a review[END_REF]. Nevertheless, none of the methods proposed in the low-frequency range give a complete answer to the following difficulties:

• the viscous effects greatly influence numerical stability, imposing a restrictive time step. In some physically relevant cases, computations cannot be carried out in a reasonable time;

• the wavelength of the slow compressional wave is much smaller than that of the other waves. Consequently, one faces the following alternative: either a coarse grid well-suited to the fast wave is chosen, and the slow wave is badly discretized; either a fine mesh is used, and the computational cost increases dramatically;

• maximum computational efficiency is obtained on a Cartesian grid; in counterpart, the interfaces are coarsely discretized, which yields spurious solutions. Alternatively, unstructured meshes adapted to the interfaces provide accurate description of geometries and jump conditions; however, the computational effort greatly increases, due to the cost of the mesh generation and to the CFL condition of stability.

The aim of the present study is to develop an efficient numerical strategy to remove these drawbacks. A time-splitting is used along with a fourth-order ADER scheme [START_REF] Schwartzkopff | Fast high order ADER schemes for linear hyperbolic equations[END_REF] to integrate Biot's equations. A flux-conserving spacetime mesh refinement [START_REF] Berger | Adaptive mesh refinement using wave-propagation algorithms for hyperbolic systems[END_REF] is implemented around the interfaces to capture the slow compressional wave. Lastly, an immersed interface method [START_REF] Lombard | Numerical treatment of two-dimensional interfaces for acoustic and elastic waves[END_REF][START_REF] Lombard | Numerical modeling of elastic waves across imperfect contacts[END_REF] is developed to provide a subcell resolution of the interfaces and to accurately enforce the jump conditions between the different porous media. As illustrated by the simulations, the combination of these numerical methods highlights the importance of an accurate modeling of the slow wave.

This article, which generalizes a previous one-dimensional work [START_REF] Chiavassa | Numerical modeling of 1-D transient poroelastic waves in the low-frequency range[END_REF], is organized as follows. Biot's model is briefly recalled in section 2. The numerical methods are described in section 3. Section 4 presents numerical experiments and comparisons with exact solutions. In section 5, conclusions are drawn and future perspectives are suggested.

Physical modeling

Biot's model

Biot's model describes the propagation of mechanical waves in a porous medium consisting of a solid matrix saturated with fluid circulating freely through the pores [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I: Low-frequency range[END_REF][START_REF] Bourbié | Acoustics of Porous Media[END_REF][START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF][START_REF] Carcione | Computational poroelasticity -a review[END_REF]. It is assumed that

• the wavelengths are large compared with the diameter of the pores;

• the amplitudes of perturbations are small;

• the elastic and isotropic matrix is fully saturated by a single fluid phase. This model relies on 10 physical parameters: the density ρ f and the dynamic viscosity η of the fluid; the density ρ s and the shear modulus µ of the elastic skeleton; the porosity 0 < φ < 1, the tortuosity a ≥ 1, the absolute permeability κ, the Lamé coefficient λ f and the two Biot's coefficients β and m of the saturated matrix. The unknowns are the elastic and acoustic displacements u s and u f , the elastic stress tensor σ, and the acoustic pressure p. In one hand, the constitutive laws are:

σ = (λ f tr ε -β m ξ) I + 2 µ ε, p = m (-β tr ε + ξ) , (1)
where I is the identity, ξ is the rate of fluid change, and ε is the strain tensor

ξ = -∇. (φ (u f -u s )) , ε = 1 2 ∇ u s + ∇ u T s . (2)
The symmetry of σ in (1) implies compatibility conditions between spatial derivatives of ε, leading to the Beltrami-Michell equation [START_REF] Rice | Some basic stress diffusion solutions for fluid-saturated elastic porous media with compressible constituents[END_REF][START_REF] Coussy | Mechanics of Porous Continua[END_REF]]

∂ 2 σ 12 ∂ x ∂ y = θ 0 ∂ 2 σ 11 ∂ x 2 + θ 1 ∂ 2 σ 22 ∂ x 2 + θ 2 ∂ 2 p ∂ x 2 + θ 1 ∂ 2 σ 11 ∂ y 2 + θ 0 ∂ 2 σ 22 ∂ y 2 + θ 2 ∂ 2 p ∂ y 2 , θ 0 = - λ 0 4 (λ 0 + µ) , θ 1 = λ 0 + 2 µ 4 (λ 0 + µ) , θ 2 = µ β 2 (λ 0 + µ) , (3)
where λ 0 = λ f -β 2 m is the Lamé coefficient of the dry matrix.

On the other hand, the conservation of momentum yields

     ρ ∂ v s ∂ t + ρ f ∂ w ∂ t = ∇σ, ρ f ∂ v s ∂ t + ρ w ∂ w ∂ t + η κ w = -∇p, (4)
where v s = ∂us ∂t = (v s1 , v s2 ) T is the elastic velocity, and w = φ ∂ ∂ t (u fu s ) = (w 1 , w 2 ) T is the filtration velocity. To be valid, the second equation of (4) requires that the spectrum of the waves lies mainly in the low-frequency range, involving frequencies lower than

f c = η φ 2 π a κ ρ f . (5)
If f ≥ f c , more sophisticated models are required [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. II: High-frequency range[END_REF][START_REF] Lu | Wave field simulation for heterogeneous porous media with singular memory drag force[END_REF]. In practice, the viscosity of the fluid is always non-zero; nevertheless, considering η = 0 can be relevant for two reasons:

• if f ≫ f c , the viscous forces are smaller than the inertial forces [START_REF] Feng | High-frequency acoustic properties of a fluid/porous solid interface. I. New surface mode[END_REF][START_REF] Morency | Spectral-element simulations of wave propagation in porous media[END_REF] and can be neglected to a first approximation;

• the exact solutions of poro-elastodynamic equations are computed more accurately if the saturating fluid is inviscid, which is attractive to validate the numerical methods.

Evolution equations

A velocity-stress formulation is followed: from ( 1) and ( 4), we obtain the system of PDEs

                                                               ∂ v s1 ∂ t - ρ w χ ∂ σ 11 ∂ x + ∂ σ 12 ∂ y - ρ f χ ∂ p ∂ x = ρ f χ η κ w 1 , ∂ v s2 ∂ t - ρ w χ ∂ σ 12 ∂ x + ∂ σ 22 ∂ y - ρ f χ ∂ p ∂ y = ρ f χ η κ w 2 , ∂ w 1 ∂ t + ρ f χ ∂ σ 11 ∂ x + ∂ σ 12 ∂ y + ρ χ ∂ p ∂ x = - ρ χ η κ w 1 , ∂ w 2 ∂ t + ρ f χ ∂ σ 12 ∂ x + ∂ σ 22 ∂ y + ρ χ ∂ p ∂ y = - ρ χ η κ w 2 , ∂ σ 11 ∂ t -(λ f + 2 µ) ∂ v s1 ∂ x -β m ∂ w 1 ∂ x -λ f ∂ v s2 ∂ y -β m ∂ w 2 ∂ y = f σ 11 , ∂ σ 12 ∂ t -µ ∂ v s2 ∂ x + ∂ v s1 ∂ y = f σ 12 , ∂ σ 22 ∂ t -λ f ∂ v s1 ∂ x -β m ∂ w 1 ∂ x -(λ f + 2 µ) ∂ v s2 ∂ y -β m ∂ w 2 ∂ y = f σ 22 , ∂ p ∂ t + m β ∂ v s1 ∂ x + ∂ w 1 ∂ x + β ∂ v s2 ∂ y + ∂ w 2 ∂ y = f p , (6)
where f σ 11 , f σ 12 , f σ 22 and f p are force densities,

ρ w = a φ ρ f , ρ = φ ρ f +(1-φ) ρ s , and χ = ρ ρ w -ρ 2 f > 0. Setting U = (v s1 , v s2 , w 1 , w 2 , σ 11 , σ 12 , σ 22 , p) T , F = (0, 0, 0, 0, f σ 11 , f σ 12 , f σ 22 , f p ) T , (7)
equations ( 6) are written as a first-order non-homogeneous linear system

∂ ∂ t U + A ∂ ∂ x U + B ∂ ∂ y U = -S U + F, (8)
where A, B and S are 8 × 8 real matrices detailled in Appendix A.

The energy of poroelastic waves can be deduced from [START_REF] Bourbié | Acoustics of Porous Media[END_REF]. Without any source terms (F = 0) and setting C ε = σ + β p I, it is proven in [START_REF] Ezziani | Modélisation de la propagation d'ondes dans les milieux viscoélastiques et poroélastiques[END_REF] that

E(t) = 1 2 R 2 (ρv s 2 +C ε : ε) dS+ 1 2 R 2 (ρ w w 2 + 1 m p 2 ) dS+ R 2 ρ f v s .w dS (9)
is an energy that satisfies

d E d t = - R 2 η κ w 2 dS. (10)
Consequently, E is conserved when the viscous effects are neglected (η = 0) and is a decreasing function otherwise. The physical parameters defined in section 2.1 are piecewise constant and can be discontinuous across interfaces. In what follows, we will focus on two domains Ω 0 and Ω 1 , which are separated by a stationary interface Γ described by a parametric equation (x(τ ), y(τ )) (figure 1). At any point P on Γ, the unit tangential vector t and the unit normal vector n are

t = 1 x ′ 2 + y ′ 2 x ′ , y ′ T , n = 1 x ′ 2 + y ′ 2 y ′ , -x ′ T . (11)
The derivatives x ′ = d x d τ and y ′ = d y d τ are assumed to be continuous everywhere along Γ, and to be differentiable as many times as required.

The evolution equations ( 6) must be completed by a set of jump conditions. The simple case of perfect bonding and perfect hydraulic contact along Γ is considered here, modeled by the jump conditions [START_REF] Gurevich | Interface conditions for Biot's equations of poroelasticity[END_REF]:

[v s ] = 0, [w.n] = 0, [σ.n] = 0, [p] = 0. ( 12
)
Enforcing these conditions is one of the main objective of the immersed interface method presented in section 3.3.

Dispersion analysis

(a) (b) The eigenvalues of A and B in (8) are real: ±c pf , ±c ps , ±c s , and 0 (multiplicity 2), where c pf > max(c s , c ps ) > 0. If η = 0, the spectral radius R(S) = η κ ρ χ can be very large and then the system (8) is stiff. A plane wave d e i(ω t-k. r) is injected in [START_REF] Bourbié | Acoustics of Porous Media[END_REF], where k = k e and d are the wavevector and the polarization, respectively; r is the position, ω = 2 π f is the angular frequency and f is the frequency. If d is collinear with k, the dispersion relation of compressional waves is obtained:

0
A k 4 + B(ω) k 2 + C(ω) = 0, A = κ m λ f + 2 µ -β 2 m , C(ω) = χ κ ω 4 -i η ρ ω 3 , B(ω) = -κ ((λ f + 2 µ) ρ w + m (ρ -2 β ρ f )) ω 2 + i η (λ f + 2 µ) ω, (13)
where the roots ±k pf and ±k ps satisfy 0 < ℜe {k pf } < ℜe {k ps }. If d is orthogonal with k, the dispersion relation of the shear wave is obtained:

k = 1 √ µ A C -B 2 C 1/2 , A = ω 2 (ρ + φ ρ f (a -2)) -i ω φ 2 η κ , B = -ω 2 φ ρ f (a -1) + i ω φ 2 η κ , C = ω 2 φ ρ f a -i ω φ 2 η κ , (14)
where the roots are ±k s , ℜe {k s } > 0. Based on ( 13) and ( 14), the phase velocities c = ω/ℜe {k} and the attenuations α = ℑm {k} of each wave are defined. In the remainder of this article, the subscripts pf , ps and s denote fast compressional, slow compressional and shear waves, respectively. The phase velocities c pf (f ), c ps (f ) and c s (f ) are monotonically increasing functions, tending asymptotically towards the eigenvalues c pf , c ps and c s . If η = 0, the three waves are non-dispersive and non-dissipative, and the energy of poroelastic waves (9) is conserved. If η = 0, the fast compressional wave and the shear wave are weakly dispersive and dissipative. The slow compressional wave, however, is highly modified by the viscosity of the saturating fluid. If f ≪ f c , then c ps (f ) ≪ c ps , and the slow compressional wave tends towards a static diffusive mode [START_REF] Chandler | The equivalence of quasi-static flow in fluid-saturated porous media and Biot's slow wave in the limit of zero frequency[END_REF]. At greater frequencies, c ps is larger but the attenuation increases. These properties are summarized in figure 2.

In the low-frequency range, the direct contribution of the slow compressional wave to the overall wave propagation processes is therefore negligible when considering an homogeneous medium. However, the influence of the slow wave becomes crucial in heterogeneous media [START_REF] Bourbié | Acoustics of Porous Media[END_REF]. The slow compressional wave, generated during the interaction between the propagative waves and the scatterers, remains localized around the interfaces. Consequently, this slow wave has a major influence on the balance equations at the interfaces, modifying crucially the behavior of fast compressional and shear diffracted waves. An accurate computation of the slow wave is therefore necessary, as shown in the numerical tests of section 4.

Numerical modeling

Numerical scheme

To integrate the system (8), a uniform grid is introduced; the spatial mesh sizes are ∆ x, ∆ y, and the time step is ∆ t. A straightforward discretization of (8) by an explicit time scheme leads to the following stability condition:

∆t ≤ min Θ ∆x c pf , 2 R(S) , ( 15
)
where Θ is obtained by a Von-Neumann analysis of stability when S = 0. In [START_REF] Dupuy | Wave propagation in heterogeneous porous media formulated in the frequencyspace domain[END_REF], the bound induced by the spectral radius of S can be very restrictive. In sandstone saturated with bitumen, for example, the maximal CFL number is roughly c pf ∆t/∆x ≈ 10 -12 ≪ Θ, which is intractable for computations.

We follow here a more efficient strategy based on second-order Strang's splitting [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF], solving alternatively the hyperbolic system

∂ ∂ t U + A ∂ ∂ x U + B ∂ ∂ y U = 0, (16)
and then the diffusive system with a source term

∂ ∂ t U = -S U + F. (17)
The linear system ( 16) is solved by applying any scheme for hyperbolic systems, giving U n+1/2 i,j

. In the numerical experiments performed in section 4, a fourth-order ADER scheme is used [START_REF] Schwartzkopff | Fast high order ADER schemes for linear hyperbolic equations[END_REF], which involves a centered stencil of 25 nodes. On Cartesian grids, this scheme amounts to a fourth-order Lax-Wendroff scheme [START_REF] Lörcher | Lax-Wendroff-type schemes of arbitrary order in several space dimensions[END_REF]. It is dispersive of order 4 and dissipative of order 6, and its stability limit is Θ = 1 [START_REF] Strikwerda | Finite Difference Schemes and Partial Differential Equations[END_REF][START_REF] Lombard | Modélisation Numérique de la Propagation et de la Diffraction d'Ondes Mécaniques[END_REF]. Other single-grid schemes can be used without any restrictions.

Since the physical parameters do not vary with time, the diffusive system ( 17) is solved exactly. For simplicity, null force density is taken: F = 0. In this case, p and σ are unchanged, whereas the velocities become (k = 1, 2)

v n+1 k = v n+1/2 k + ρ f ρ 1 -e -η κ ρ χ T w n+1/2 k , w n+1 k = e -η κ ρ χ T w n+1/2 k , (18)
where T depends on the time step (see section 3.4). The splitting ( 16)-( 17) along with exact integration (18) recovers the optimal condition of stability:

c pf ∆ t/∆ x ≤ Θ.
Since the matrices A and B do not commute with S, the theoretical order of convergence falls from 4 to 2 when the viscosity is non-negligible.

Using a fourth-order accurate scheme such as ADER 4 is nevertheless advantageous, compared with a second-order scheme such as Lax-Wendroff: the stability limit is improved, and numerical artifacts (dispersion, attenuation, anisotropy) are greatly reduced.

In [START_REF] De La Puente | Discontinuous Galerkin methods for wave propagation in poroelastic media[END_REF], the authors notice that the first-order splitting does not lead to a correct representation of the slow mode at low frequencies. Nevertheless, the numerous one-dimensional examples provided in [START_REF] Chiavassa | Numerical modeling of 1-D transient poroelastic waves in the low-frequency range[END_REF] demonstrate that the second-order splitting accurately represents the static mode when a sufficient number of discretization points per wavelength is used. This can be obtained by using a local space-time refinement presented in the following section.

Mesh refinement

The slow wave has much smaller spatial scales of evolution than the wavelength of the other waves. A very fine grid is therefore required to account for its evolution. Since the use of a fine uniform grid on the whole computational domain is out of reach, grid refinement provides a good alternative. In addition, the slow wave remains localized near the interfaces (section 2.4), and hence grid refinement is necessary only around these places. Lastly, even if the slow wave propagates (η = 0) the property c pf ≫ c ps is usually satisfied: consequently, a fine mesh near the interface is still useful to perform accurate extrapolations, as required by the immersed interface method (section 3.3).

We adopt here a space-time mesh refinement approach based on flux conservation [START_REF] Berger | Adaptative mesh refinement for hyperbolic partial differential equations[END_REF][START_REF] Berger | Adaptive mesh refinement using wave-propagation algorithms for hyperbolic systems[END_REF], which is more naturally coupled to the flux-conserving scheme developed to solve [START_REF] Ezziani | Modélisation de la propagation d'ondes dans les milieux viscoélastiques et poroélastiques[END_REF]. The refined zones are rectangular Cartesian patches with mesh sizes ∆ x / q, ∆ y / q, where the integer q is the refinement factor. To reduce the cpu time and to limit the numerical dispersion on the coarse grid, a local time step ∆ t / q is used [START_REF] Joly | An error analysis of conservative space-time mesh refinement methods for the one-dimensional wave equation[END_REF][START_REF] Rodriguez | Raffinement de maillage spatio-temporel pour les équations de l'élastodynamique[END_REF]. When one time step is done on the coarse grid, q time substeps are done on the refined zone. The extrapolated values required to couple coarse and fine grids are obtained by linear interpolation in space and time on the numerical values at the surrounding nodes [START_REF] Lombard | Modélisation Numérique de la Propagation et de la Diffraction d'Ondes Mécaniques[END_REF]. In the case of the Lax-Wendroff scheme applied to the scalar advection equation, the stability of the coupling is proven in [START_REF] Berger | Stability of interfaces with mesh refinement[END_REF] whatever q.

The additional cost induced by mesh refinement can become prohibitive, both concerning the memory requirements and the computational time, because of the q substeps inside one time step. The value of q must be estimated carefully in terms of the physical parameters. For this purpose, the wavelengths λ pf (f 0 ) = c pf (f 0 )/f 0 and λ ps (f 0 ) = c ps (f 0 )/f 0 are deduced from the dispersion analysis, where f 0 is the central frequency of the source. The number of fine grid nodes per wavelength of the slow compressional wave and the number of coarse grid nodes per wavelength of the fast compressional wave must then be roughly equal:

λ ps (f 0 ) ∆x / q ≈ λ pf (f 0 ) ∆x ⇒ q ≈ c pf (f 0 ) c ps (f 0 ) . ( 19
)

Immersed interface method

The discretization of the interfaces requires special care. A straightforward stair-step discretization of the interfaces introduces a first-order geometrical error and yields spurious numerical diffractions. In addition, the jump conditions [START_REF] Chiavassa | Numerical modeling of 1-D transient poroelastic waves in the low-frequency range[END_REF] are not enforced numerically if no special treatment is applied. Lastly, the smoothness requirements to solve [START_REF] Ezziani | Modélisation de la propagation d'ondes dans les milieux viscoélastiques et poroélastiques[END_REF] are not satisfied, decreasing the convergence rate of the ADER scheme.

To remove these drawbacks while maintaining the efficiency of Cartesian grid methods, we adapt an immersed interface method previously developed in acoustics and elastodynamics [START_REF] Piraux | A new interface method for hyperbolic problems with discontinuous coefficients: one-dimensional acoustic example[END_REF][START_REF] Lombard | Numerical treatment of two-dimensional interfaces for acoustic and elastic waves[END_REF][START_REF] Lombard | Numerical modeling of elastic waves across imperfect contacts[END_REF][START_REF] Lombard | Modélisation Numérique de la Propagation et de la Diffraction d'Ondes Mécaniques[END_REF]. At the irregular points where the ADER's stencil crosses the interface Γ, the scheme will use modified values of the solution, instead of the usual numerical values. The modified values are extrapolations, based on the local geometry of Γ and on r successive derivatives of the jump conditions [START_REF] Chiavassa | Numerical modeling of 1-D transient poroelastic waves in the low-frequency range[END_REF]. The parameter r is discussed at the end of this section.

Let us consider a point M(x I , y J ) ∈ Ω 1 and its orthogonal projection P onto Γ (figure 3). The algorithm to build the modified value at M is divided into four steps.

Step 1: high-order interface conditions. On the side Ω k (k = 0, 1), the boundary values of the spatial derivatives of U up to the r-th order are put in a vector U r k with n v = 4 (r + 1) (r + 2) components:

U r k = lim M →P, M ∈Ω k U T , ..., ∂ l ∂ x l-m ∂ y m U T , ..., ∂ r ∂ y r U T T , (20)
where l = 0, ..., r and m = 0, ..., l. Following this formalism, the zero-th order jump conditions (12) are written

C 0 1 U 0 1 = C 0 0 U 0 0 , (21)
where the matrices of the jump conditions C 0 k depend on the local geometry of Γ:

C 0 k (τ ) =         1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 y ′ -x ′ 0 0 0 0 0 0 0 0 y ′ -x ′ 0 0 0 0 0 0 0 y ′ -x ′ 0 0 0 0 0 0 0 0 1         . ( 22
)
The jump condition ( 21) is differentiated with respect to time t, and then the time derivatives are replaced by spatial derivatives thanks to the conservation law [START_REF] Ezziani | Modélisation de la propagation d'ondes dans les milieux viscoélastiques et poroélastiques[END_REF]. For example, we obtain

∂ ∂ t (C 0 0 U 0 0 ) = -C 0 0 A 0 ∂ ∂ x U 0 0 -C 0 0 B 0 ∂ ∂ y U 0 0 , (23)
where A 0 and B 0 are the matrices in Ω 0 . The jump condition ( 21) is also differentiated in terms of τ . Taking advantage of the chain-rule, we obtain e.g.

d d τ (C 0 0 U 0 0 ) = d d τ C 0 0 U 0 0 + C 0 0 x ′ ∂ ∂ x U 0 0 + y ′ ∂ ∂ y U 0 0 . (24)
From ( 21), ( 23) and ( 24), we build matrices

C 1 k such that C 1 1 U 1 1 = C 1 0 U 1 0
, which provides first-order jump conditions. By iterating this process r times, r-th order interface conditions are obtained

C r 1 U r 1 = C r 0 U r 0 , (25)
where C r k are n c × n v matrices (k = 0, 1), and n c = 3 (r + 1) (r + 2). The computation of matrices C r k is a tedious task when r ≥ 2, that can be greatly simplified using computer algebra tools.

Step 2: high-order Beltrami-Michell equations. The equation ( 3) is satisfied anywhere in a poroelastic medium. Under sufficient smoothness requirements, it can be differentiated with respect to x and y, as many times as required:

∂ j σ 12 ∂ x j-i-1 ∂ y i+1 = θ 0 ∂ j σ 11 ∂ x j-i ∂ y i + θ 1 ∂ j σ 22 ∂ x j-i ∂ y i + θ 2 ∂ j p ∂ x j-i ∂ y i +θ 1 ∂ j σ 11 ∂ x j-i-2 ∂ y i+2 + θ 0 ∂ j σ 22 ∂ x j-i-2 ∂ y i+2 + θ 2 ∂ j p ∂ x j-i-2 ∂ y i+2 , (26)
where j ≥ 2 and i = 0, • • • , j -2. The equations ( 26) are also satisfied along Γ. They can be used therefore to reduce the number of independent components in U r k . For this purpose, we define the vectors V r k such that

U r k = G r k V r k , (27)
where 20) and ( 26), an algorithm to compute the non-zero components of G r k is proposed in Appendix B.

G r k are n v × (n v -n b ) matrices, and n b = r (r -1)/2 if r ≥ 2, n b = 0 otherwise. Based on (
Step 3: high-order boundary values.

Based on ( 25) and ( 27), the vectors of independent boundary values satisfy

S r 1 V r 1 = S r 0 V r 0 , (28)
where

S r k = C r k G r k are n c × (n v -n b ) matrices.
Since the system ( 28) is underdetermined, the solution is not unique, and hence it can be written

V r 1 = (S r 1 ) -1 S r 0 | K S r 1 V r 0 Λ r , (29)
where (S r 1 ) -1 is the least-squares pseudo-inverse of S r 1 , K S r 1 is the matrix filled with the kernel of S r 1 , and Λ r is a set of n v -n c -n b Lagrange multipliers that represent the coordinates of V r 1 onto the kernel. A singular value decomposition of S r 1 is used to build (S r 1 ) -1 and the kernel K S r 1 [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF].

Step 4: construction of modified values.

Let Π r ij be the matrix of r-th order 2D Taylor expansions Π r i,j = 1, ...,

1 l ! (l -m) ! (x i -x P ) l-m (y j -y P ) m , ..., (y j -y P ) r r ! I 8 , (30)
where I 8 is the 8 × 8 identity matrix, l = 0, ..., r and m = 0, ..., l. The modified value at (x I , y J ) is a smooth extension of the solution on the other side of Γ (figure 3), and it writes:

U * I,J = Π r I,J U r 0 = Π r I,J G r 0 V r 0 . (31)
The vector V r 0 in [START_REF] Martin | An unsplit convolutional perfectly matched layer improved at grazing incidence for seismic wave propagation in poroelastic media[END_REF] remains to be estimated in terms of the boundary conditions and of the numerical values at surrounding grid points. purpose, we consider the disc D centered at P with a radius d, that contains N d grid points. At the grid points of D ∩ Ω 0 , r-th order Taylor expansion of the solution at P gives

U(x i , y j , t n ) = Π r i,j U r 0 + O(∆ x r+1 ), = Π r i,j G r 0 (1 | 0) V r 0 Λ r + O(∆ x r+1 ). ( 32
)
At the grid points of D ∩ Ω 1 , r-th order Taylor expansion of the solution at P and the boundary conditions [START_REF] Lu | Wave field simulation for heterogeneous porous media with singular memory drag force[END_REF] give

U(x i , y j , t n ) = Π r i,j U r 1 + O(∆ x r+1 ), = Π r i,j G r 1 (S r 1 ) -1 S r 0 | K S r 1 V r 0 Λ r + O(∆ x r+1 ). (33)
Equations ( 32) and [START_REF] Masson | Finite-difference modeling of Biot poroelastic equations across all frequencies[END_REF] are written in the matrix form

(U(., t n )) D = M V r 0 Λ r +    O(∆ x r+1 ) . . . O(∆ x r+1 )    , ( 34
)
where M is a convenient 8

N d × (2 n v -2 n b -n c ) matrix.
To ensure that the system (34) is overdetermined, the radius d of the disc is chosen to satisfy

ε(d, r) = 8 N d 2 n v -2 n b -n c ≥ 1. ( 35
)
Exact values in [START_REF] Morency | Spectral-element simulations of wave propagation in porous media[END_REF] are replaced by numerical ones, and the Taylor rests are removed. The least-squares inverse of M is denoted by M -1 . The Lagrange multipliers Λ k are accounted in the construction of M, but are not involved in the definition of the modified value [START_REF] Martin | An unsplit convolutional perfectly matched layer improved at grazing incidence for seismic wave propagation in poroelastic media[END_REF]. As a consequence, they can be removed and the

(n v -n b ) × 8 N d restriction M -1 of M -1 is defined by V r 0 = M -1 (U n ) D . (36)
Lastly, the modified value follows from ( 31) and ( 36):

U * I,J = Π r I,J G r 0 M -1 (U n ) D . ( 37
)
Comments and practical details.

1. A similar algorithm is applied at each irregular point along Γ. The sizes of the matrices involved are summarized in table 1. Since the jump conditions do not vary with time, the evaluation of the matrices in ( 37) is done during a preprocessing step. Only small matrix-vector products are therefore required at each time step. After optimization of the computer codes, this additional cost is made negligible, lower than 1% of the time-marching. 2. The matrix M in (34) depends on the subcell position of P inside the mesh and on the jump conditions at P , involving the local geometry and the curvature of Γ at P . Consequently, all these insights are incorporated in the modified value [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF], and hence in the scheme.

Quantity Size

n v 4 (r + 1) (r + 2) n c 3 (r + 1) (r + 2) n b r (r -1) / 2 if r ≥ 2, 0 else C r k n c × n v G r k n v × (n v -n b ) S r k n c × (n v -n b ) Π r k 8 × n v 1 (n v -n b ) × (n v -n b ) 0 (n v -n b ) × (n v -n c -n b ) M 8 N d × (2 n v -2 n b -n c ) M -1 (n v -n b ) × 8 N d
Table 1: Quantities involved in the computation of the modified values (section 3.3).

3.

The simulations indicate that overestimation of ε in ( 35) has a crucial influence on the stability of the immersed interface method. Various strategies can be used to ensure [START_REF] Piraux | A new interface method for hyperbolic problems with discontinuous coefficients: one-dimensional acoustic example[END_REF], for instance an adaptive choice of d depending on the local geometry of Γ at P . We adopt here a simpler strategy, based on a constant radius d. Taking r = 2, numerical experiments have shown that d = 3.2 ∆ x is a good candidate, while d = 4.5 ∆ x is used when r = 3. In this case, we obtain typically N d ≈ 20 and ε ≈ 4. 4. The order r plays an important role on the accuracy of the coupling between the immersed interface method and a s-th order scheme. If r ≥ s, then a s-th order local truncation error is obtained at the irregular points. This condition can be slightly relaxed: r = s-1 still ensures a sth order overall accuracy [START_REF] Gustafsson | The convergence rate for difference approximations to mixed initial boundary value problems[END_REF]. As a consequence, a fourth-order ADER scheme (s = 4) requires a third-order immersed interface method (r = 3) to maintain fourth-order convergence. 5. A GKS analysis of stability has been performed in 1D in the case of an inviscid saturating fluid [START_REF] Lombard | Modélisation Numérique de la Propagation et de la Diffraction d'Ondes Mécaniques[END_REF]. Extending this approach to 2D problems with viscous saturating fluids is out of reach. Various numerical experiments, however, indicate the stability of the method under the usual CFL condition (section 3.1), if two requirements are satisfied: (i) the number of grid nodes used for extrapolations is sufficiently large, as stated in point 3; (ii) the Beltrami-Michell equations ( 27) are used.

Summary of the algorithm

The numerical strategy presented in this section couples three numerical methods: a finite difference numerical scheme with splitting (section 3.1), a space-time mesh refinement (section 3.2), and an immersed interface method (section 3.3). To clarify the interactions between these methods, the global algorithm is summarized as follows: ⊲ Preprocessing -Detection of irregular grid points -Computation of extrapolation matrices in [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF] -Initialization of the solution at t = 0 -Diffusive step [START_REF] Gurevich | Interface conditions for Biot's equations of poroelasticity[END_REF] where T = ∆ t / 2 on the coarse grid -Diffusive step [START_REF] Gurevich | Interface conditions for Biot's equations of poroelasticity[END_REF] where T = ∆ t / (2 q) on the refined grids ⊲ Time iterations

• Coarse grid:

-Computation of modified values (37) if present -Solving the propagative step ( 16) -Diffusive step [START_REF] Gurevich | Interface conditions for Biot's equations of poroelasticity[END_REF] where T = ∆ t

• On each refined grid, q subtime iterations:

-Space-time interpolations at the grid boundaries -Computation of modified values [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF] -Solving the propagative step ( 16) -Diffusive step [START_REF] Gurevich | Interface conditions for Biot's equations of poroelasticity[END_REF] where T = ∆ t / q ⊲ End of time iterations -Diffusive step [START_REF] Gurevich | Interface conditions for Biot's equations of poroelasticity[END_REF] where T = ∆ t / 2 on the coarse grid -Diffusive step [START_REF] Gurevich | Interface conditions for Biot's equations of poroelasticity[END_REF] where T = ∆ t / (2 q) on the refined grids

Numerical experiments

Configurations

Five tests are proposed along this section. In Test 1, the convergence order of the ADER scheme coupled with the immersed interface method is measured. Test 2 illustrates the different kind of waves in homogeneous media, and also the influence of the local space-time refinement. Test 3 investigates the numerical stability of the global algorithm. Diffraction of a plane wave by one (Test 4) and four (Test 5) cylindrical scatterers illustrates the accuracy and the physical relevance of the proposed numerical methods.

The physical parameters given in table 2 correspond to Cold Lake sandstone and shale saturated with water [START_REF] Dai | Wave propagation in heterogeneous porous media: a velocity-stress, finite-difference method[END_REF], respectively. In some experiments, an inviscid saturating fluid is artificially considered: η = 0 Pa.s, the other parameters being unchanged. As recalled in section 2.1, this limit-case has physical significance only in the high-frequency range. It is mainly addressed here for a numerical purpose. Once the spatial mesh sizes ∆ x and ∆ y are chosen on the coarse grid, the time step follows from the CFL number in Ω 0 : c pf 0 ∆ t/ max (∆ x, ∆ y) = 0.95 < 1. If η = 0, the maximum CFL number induced by ( 15) is equal to 0.5: consequently, the simulations done here with the splitting ( 16)-( 17) are twice faster than with unsplitted methods.

Parameters Ω 0 Ω 1 ρ s (kg/m 3 )
The grids are excited by two means: either a plane fast compressional wave, either a point source that generates cylindrical waves. Details of the excitation method are given in Appendix C. In the case of an incident plane wave, the exact expression given in Appendix C is also enforced numerically on the edges of the computational domain. No special attention is paid to simulate outgoing waves, for instance with Perfectly-Matched Layers [START_REF] Martin | An unsplit convolutional perfectly matched layer improved at grazing incidence for seismic wave propagation in poroelastic media[END_REF][START_REF] Zeng | The application of the perfectly matched layer in numerical modeling of wave propagation in poroelastic media[END_REF]. In all the presented tests, the size of the domain and the number of iterations are chosen to avoid the spurious reflections of diffracted waves with the outer frontiers.

Test 1: convergence measurements

In Test 1, we focus on the coupling between the ADER scheme (section 3.1) and the immersed interface method (section 3.3). For this purpose, we consider a domain [0, 400] m 2 cut by a plane interface with slope 60 degrees. The saturating fluids are inviscid: exact solutions can be computed very accurately without Fourier synthesis, and splitting errors of the scheme are avoided. The source is plane (C.2), with parameters: γ = 10 -3 , f 0 = 40 Hz, t = 3.3 10 -2 s, and θ = -30 degrees (figure 4-a). Consequently, the incident wave propagates normally to the interface, leading to a 1-D configuration; from a numerical point of view, however, the problem is fully bidimensional.

The computations are done on a uniform grid of N × N points, during N/4 time steps. Comparisons with the exact values of the pressure p are done on the subdomain [50, 350] m × [150, 250] m, in order to avoid spurious effects induced by the edges of the computational domain (figure 4-b). The measures involve reflected and transmitted fast and slow compressional waves generated by the interface (no shear wave is generated in 1D); these waves are highly sensitive to the discretization of the jump conditions.

Errors in l 2 norm and convergence rates are reported in table 3 and drawn on figure 5. Various values of r are investigated; r = 0 means that no immersed interface method is applied; in this case, first-order accuracy is obtained. As stated in section 3.3, fourth-order accuracy is maintained if r = 4 -1 = 3, i.e if third-order extrapolations are used in the immersed interface method. In the present test case, r = 2 is sufficient to obtain the same level of accuracy on a large range of grid size. Nevertheless, this could be untrue in other contexts, and hence we will always use r = 3 in the following simulations.

Test 2: mesh refinement

In the second test, we focus on the coupling between the ADER scheme (section 3.1) and the mesh refinement (section 3.2). For this purpose, a homogeneous medium Ω 0 on a domain [-250, 250] Figure 6 shows snapshots of the fields after 280 iterations. At this time t ≃ 4.4 / f 0 , the fast compressional wave has also crossed the grid G 2 . In the left column, one takes η = 0, and consequently all the waves generated by the source propagate radially. In the right column, the viscosity is non zero, and the slow compressional wave remains localized around the source point. This wave is clearly observed on the pressure field, but is also present in the stress fields although not visible on the corresponding panels of figure 6.

N r = 0 order r = 1 order r = 2 order r = 3 order 400 4.894 10 0 -6.527 10 0 -6.107 10 0 -5.067 10 0 -800 1.247 10 0 1.973 1.667 10 0 1.961 1.065 10 0 2.520 8.642 10 -1 2.552 1200 6.520 10 -1 1.599 6.758 10 -1 2.242 2.273 10 -1 3.809 1.936 10 -1 3.690 1600 4.770 10 -1 1.086 3.617 10 -1 2.173 7.995 10 -2 3.632 6.835 10 -2 3.619 2000 3.818 10 -1 0.998 2.254 10 -1 2.119 3.344 10 -2 3.906 2.888 10 -2 3.861 2400 3.128 10 -1 1.093 1.509 10 -1 2.201 1.564 10 -2 4.168 1.414 10 -2 3.917 2800 2.696 10 -1 0.964 1.097 10 -1 2.069 9.157 10 -3 3.473 7.697 10 -3 3.945 3200 2.390 10 -1 0.902 8.431 10 -2 1.971 5.598 10 -3 3.685 4.504 10 -3 4.013 3600 2.115 10 -1 1.038 6.573 10 -2 2.114 3.613 10 -3 3.718 2.840 10 -3 3.915 4000 1.928 10 -1 0.879 5.362 10 -2 1.933 2.485 10 -3 3.552 1.866 10 -3 3.986 Table 3: Test 1. Convergence rate in l 2 norm. No immersed interface method (r = 0), linear (r = 1), quadratic (r = 2) or cubic (r = 3) immersed interface method.

In order to evaluate the influence of the refined areas, a reference solution is computed using a fine mesh ∆x = ∆y = 1/5 m on the whole computational domain. As observed in figure 7-(a),(b), the local space-time refinement applied in G 1 and G 2 does not lead to significant spurious reflections during the propagation of the waves. On the contrary and as expected, the refinement around the source point leads to a much better resolution of the slow wave in the viscous case: see figure 7-(c).

In conclusion, mesh refinement coupled with the ADER scheme (with or without splitting) accurately represents the behavior of the different poroelastic waves in an homogeneous medium, even though they present a spatially complex structure, as observed in figure 7-(d).

Test 3: stability of the complete algorithm

Some theoretical results based on GKS theory exist concerning the stability of immersed interface methods or local mesh refinement [START_REF] Gustafsson | The convergence rate for difference approximations to mixed initial boundary value problems[END_REF][START_REF] Berger | Stability of interfaces with mesh refinement[END_REF][START_REF] Trefethen | Stability of finite-difference models containing two boundaries or interfaces[END_REF][START_REF] Lin | Some grid refinement schemes for hyperbolic equations with piecewise constant coefficients[END_REF][START_REF] Berger | Adaptive mesh refinement using wave-propagation algorithms for hyperbolic systems[END_REF][START_REF] Lombard | Modélisation Numérique de la Propagation et de la Diffraction d'Ondes Mécaniques[END_REF]. However, these analyses have been done mainly in the case of one-dimensional model problems and basic numerical schemes. The present algorithm combines more sophisticated numerical methods in 2D, and hence GKS analysis is out of reach. The only reasonable way to confirm the stability of the full method is to perform numerical perturbation tests.

For this purpose, we consider a heterogeneous porous media made up by a matrix Ω 0 with a cylindrical scatterer Ω 1 of radius r = 10 m, centered in a domain [-50, 50] m 2 . In both Ω 0 and Ω 1 , the viscosity of the saturating fluids are taken into account (η = 0). All the components of U (7) are initialized randomly at each grid points, and the source F is set to zero. The computational domain [-50, 50] m 2 is discretized by a coarse grid. The scatterer is included in a grid G 1 of size 15 × 15 m 2 with a large refinement factor q = 9.

In figure 8, the energy of poroelastic waves ( 9) is displayed during 3 10 4 iterations. Theoretically, this energy should slowly decrease, depending on w in [START_REF] Chandler | The equivalence of quasi-static flow in fluid-saturated porous media and Biot's slow wave in the limit of zero frequency[END_REF]. At the beginning of the simulations, a large decrease of E is observed. It is logically induced by the random non-smooth initial field, which generates a large numerical dissipation. After roughly 1000 time steps, a smooth field is reached and the mechanical energy slowly decreases. This confirms the stability of the full numerical method, involving the ADER scheme with splitting, mesh refinement, and immersed interface method.

Test 4: diffraction of a plane wave by a cylinder

A cylindrical scatterer Ω 1 of radius 40 m is centered at point (0, 0) in the computational domain [-200, 200] m 2 . The source is a plane wave (C.2) initially in medium Ω 0 , with parameters: γ = 1, f 0 = 40 Hz, t = -2.09 10 -2 s, and θ = 0 degree. The initial conditions are illustrated in figure 9 (a-b).

In figure 9, the viscosity has been canceled in both media Ω 0 and Ω 1 . In such a configuration, the analytical solution of (7) can be computed using Fourier and Bessel expansions and is used to validate the simulations. The diffracted waves propagate with velocities c pf , c ps and c s given in table 2. To ensure the same number of points per wavelength for all the diffracted waves (see section 3.2), the computational domain is locally refined by inserting the cylinder in the grid G 1 = [-45, 45] 2 with a refinement factor q = 5, deduced from relation [START_REF] Gustafsson | The convergence rate for difference approximations to mixed initial boundary value problems[END_REF] and table 2. The influence of this refinement combined with the immersed interface method is clearly visible in figure 9.

Without refinement nor interface method (q = 1, r = 0, left column), the waves created during the interaction with the scatterer are polluted by spurious numerical artifacts. With q = 5 and r = 3 (right column), these non-physical perturbations disappear and the reflected-transmitted waves are correctly computed. Comparisons with the exact solution presented on figure 10 confirm the accuracy of the simulation. Without mesh refinement nor immersed interface method, inaccurate results are obtained, especially for the reflected fast compressional wave where a shift of about 5 m is observed (figure 10-(b)).

In figure 10 -(c) and (d), the time evolution of the pressure registered at the point (-60, 60) is presented from t = 0.028 s, in order to avoid the incident initial wave. We observe the fast reflected wave, followed by the slow wave and a combination of diffracted fast and slow waves. Once again, the full strategy captures accurately all the temporal variations of the solution, while the basic algorithm gives very poor results, especially concerning the waves diffracted by the scatterer (see figure 10-(d)).

The same configuration is now considered by taking into account the viscosity of the saturating fluids (see table 2). Based on the dispersion analysis performed in section 2.4, the values of the phase velocity at f 0 = 40 Hz give a refinement factor q ≈ 22 in the grid G 1 to satisfy our refinement criterion [START_REF] Gustafsson | The convergence rate for difference approximations to mixed initial boundary value problems[END_REF]. Snapshots of p and σ 11 after 180 iterations are given on figure 11, showing different structures than in the inviscid case presented in figure 9. The diffracted fast compressional and shear waves propagate, while the slow compressional wave remains localized around the scatterer.

The spatial structure of this diffusive wave is clearly visible in figure 11-(c),(d). Figure 11-(e) focuses on the pressure component of the slow wave around the interface. It is observed that the numerical method applied without mesh refinement nor interface treatment gives inaccurate results about this static wave. This is also the case of the propagative reflected compressional wave shown in figure 11-(f). The computed fields converge when the refinement factor in G 1 increases. With refinement factors q larger than 15, the results are mostly indistinguishable, and hence they are not represented here. With q = 15, the refined grid G 1 contains 1365 2 discretization points and involves 19184 irregular points in the immersed interface algorithm.

The time evolution of the mechanical energy ( 9) is illustrated in figure 12. The times t i = 0.008 s and t f ≃ 0.066 s correspond respectively to the instants when the plane wave begins to interact with the scatterer and when it has crossed the scatterer entirely. In the case of inviscid saturating fluids (blue circles on figure 12), the energy is almost conserved, as expected, and only a decrease of 0.15 % is observed. In the viscous case (blue diamond on figure 12), the energy decreases when the plane wave interacts with the scatterer, and remains constant otherwise. This behavior is logical, since the rate of decrease of energy in [START_REF] Chandler | The equivalence of quasi-static flow in fluid-saturated porous media and Biot's slow wave in the limit of zero frequency[END_REF] is governed by w. In homogeneous medium, w is extremely small compared with other fields, and hence d E d t ≈ 0. On the other hand, during the interaction between the wave and the scatterer, the filtration velocity of the slow wave generated at the interface has an important amplitude. The theoretical decrease of energy, obtained by a numerical integration of equation [START_REF] Chandler | The equivalence of quasi-static flow in fluid-saturated porous media and Biot's slow wave in the limit of zero frequency[END_REF], is shown using a dotted line in figure 12 and is very close to the observed decrease of E. This confirms that our numerical strategy accurately models the dissipation of mechanical energy. 9) for inviscid and viscous case. Numerical integration of relation ( 10) is also represented by the dashed curve.

Test 5: multiple scatterers

In the last test, the ability of the proposed numerical strategy to handle complex geometries is illustrated. Four cylindrical scatterers of medium Ω 1 are inserted in the matrix Ω 0 . Each cylinder is surrounded by a refined grid. The refinement factor is q = 5 in the inviscid case, and q = 10 in the case of viscous saturating fluids. The pressure field is represented in figure 13. The behavior of the fast compressional waves is qualitatively the same in both cases, unlike the slow waves:

• in the inviscid case (figure 13, left column), the slow waves propagate and interact with the other scatterers, which generates new sets of reflected-transmitted waves. Accurate computation of these successive interactions is obtained thanks to the combination of the local refinement procedure and the immersed interface method;

• in the viscous case (figure 13, right column), the slow waves remain localized around the interfaces, and hence they do not participate directly to the scattering process. At time t 2 , they have mostly disappeared due to their physically diffusive behavior. As indicated by test 4, a fine modeling of these waves is nevertheless necessary in order to compute accurately the propagative diffracted waves.

t 1 t 1 t 2 t 2

Conclusion

Numerical modeling of 2-D poroelastic waves was addressed numerically in the time-domain. The evolution equations were issued from Biot's theory, assuming viscous efforts of Poiseuille type, which is valid essentially in the low-frequency range. Three numerical tools were combined: a fourth-order scheme with time-splitting, a space-time mesh refinement, and an immersed interface method. The resulting method provides highly accurate simulations of wave propagation in realistic configurations. Numerical experiments have indicated that even if the slow waves remain localized around sources or interfaces, the manner in which they are computed influences the propagative waves, that are measured in practical applications. Academic cylindrical geometries have been considered here; however, more complex smooth geometries can be handled without any restrictions, for example cubic splines. Various extensions of the present work are suggested:

• The numerical methods presented here make it possible to simulate relevant physical experiments. We have especially in mind the modeling of multiple scattering in random media. Based on simulated data, the properties of the effective medium equivalent to the disordered medium under study can be deduced [START_REF] Chekroun | Comparisons between multiple scattering methods and direct numerical simulations for elastic wave propagation in concrete[END_REF]. This numerical approach compares advantageously with the methods usually followed by physicists: real experiments are expensive, and analytical methods are restricted to very small concentrations of scatterers [START_REF] Luppé | Coherent waves in a multiply scattering poro-elastic medium obeying Biot's theory[END_REF]. The example given in section 4.6 is obviously preliminary: interaction of a plane wave with hundreds of scatterers needs to be addressed, which requires the parallelization of the algorithms.

• Poroelastic media in perfectly bonded contact are considered here. More realistic conditions can be studied and properly enforced by the immersed interface method, for instance sliding and imperfect bonding [START_REF] Lombard | Numerical modeling of elastic waves across imperfect contacts[END_REF], or imperfect hydraulic contact [START_REF] Bourbié | Acoustics of Porous Media[END_REF]. Comparison of numerical simulation with experimental results could help to validate or improve the models of contact, which constitutes a current issue in poroelasticity. Generalizing our approach to the interface between a poroelastic medium and a fluid is another direction of work with large applications, for instance in biomechanics.

• Incorporation of attenuation in the elastic skeleton is required to model the real processes of dissipation [START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF]. For this purpose, memory variables need to be introduced in the evolution equations [START_REF] Bourbié | Acoustics of Porous Media[END_REF].

• Lastly, the numerical modeling of the transient Biot equations in the full range of validity of poroelasticity constitutes a natural extension of the present work. At frequencies greater than f c in (5), a correction of the viscosity proportional with the square root of frequency needs to be introduced, as described e.g. by the JKD model [START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media[END_REF][START_REF] Dupuy | Wave propagation in heterogeneous porous media formulated in the frequencyspace domain[END_REF]. In the time domain, the new evolution equations involve fractional derivatives of order 1/2, whose efficient numerical evaluation is a major challenge [START_REF] Lu | Wave field simulation for heterogeneous porous media with singular memory drag force[END_REF][START_REF] Haddar | Efficient solution of a wave equation with fractional order dissipative terms[END_REF][START_REF] Masson | Finite-difference modeling of Biot poroelastic equations across all frequencies[END_REF].

                     
e i (ω t-k pf .r) ĥ(ω), (C.2) where γ is an amplitude factor, k pf is the wavenumber [START_REF] Coussy | Mechanics of Porous Continua[END_REF] and

Y pf (ω) = ((1 -φ) ρ s + ρ f β (a -1)) ω 2 -(λ f + 2 µ -m β 2 ) k 2 pf + i ω φ β η κ ρ f (a β -φ) ω 2 -i ω φ β η κ .
(C.3) When η = 0, k pf depends linearly on ω, and consequently Y pf and the vector column in (C.2) no more involve ω: a straightforward time-domain expression of the incident plane wave can be obtained.

As a second source, we also implement force densities acting on σ 12 in (6). The only non-null component in [START_REF] Carcione | Some aspects of the physics and numerical modeling of Biot compressional waves[END_REF] is

f σ 12 = g(x, y) h(t) (C.4)
where g is a truncated gaussian centered at point (x s , y s ): and h is a Ricker signal:

h(t) =      2 π 2 f 2 0 t - 1 f 0 2 -1 exp -2 π 2 f 2 0 t - 1 f 0 2 if 0 < t

1 ΓFigure 1 :

 11 Figure 1: Interface Γ between two poroelastic media Ω 0 and Ω 1 .

Figure 3 :

 3 Figure 3: Irregular point M (x I , y J ) ∈ Ω 1 and its orthogonal projection P onto Γ. The grid nodes used to compute U * I,J are inside the circle with radius d and centered on P ; they are denoted by +.

Figure 4 :

 4 Figure 4: Test 1. Snapshots of p at the initial instant (a) and at the instant of measure (b). The white rectangle denotes the zone where convergence errors are measured.

  m 2 is excited by the force density (C.4). The parameters of the source are: x s = y s = 0, f 0 = 40 Hz, ζ = c pf 0 15 f 0 , r 0 = 2 ζ, and γ = 10 3 . The computational domain is discretized on a coarse mesh of 500 2 points. Two locally refined areas are added: one around the source point G 1 = [-25, 25] 2 , and one at G 2 = [80, 130] × [80, 120]. Both grids are refined by a factor q = 5, leading to 255 2 points in G 1 and 205 × 255 points in G 2 .

Figure 5 :

 5 Figure 5: Test 1. Error measured in l 2 norm versus the number of points N x , with various order r of the immersed interface method. Dotted line corresponds to 4-th order slope.

12 Figure 6 :

 126 Figure 6: Test 2. Snapshots of p, σ 11 and σ 12 of the fast (p f ), slow (p s ), and shear (s) waves emitted by a source point. Left column: viscosity η = 0. Right column: viscous case. Dashed areas indicate the location of the refined grids G 1 and G 2 .

Figure 7 :

 7 Figure 7: Test 2 in the inviscid case, corresponding to the left column of figure 6. Stress σ 11 along lines y = 30 m (a) and y = 130 m (b). Comparison between the fine grid reference solution and the solution obtained with a coarse grid and two refined areas of factor q = 5 (c). 3D representation of σ 11 (d).

Figure 8 :

 8 Figure 8: Test 3. Time evolution of the computed total energy (9).

Figure 9 :

 9 Figure 9: Test 4, inviscid fluids. (a): p at initial time. (b): p along the line y = 0; vertical lines denote the frontiers of the scatterer. (c-d): p after 180 iterations. (c): no refinement nor interface method (q = 1, r = 0). (d): refinement factor q = 5 and third order immersed interface method (r = 3). Dotted square represents the frontiers of the refined domain. The cross denotes the location of the receiver.

Figure 10 :

 10 Figure 10: Test 4, inviscid fluids. (a): p on the line y = 0 after 180 iterations. (b): zoom around the reflected fast compressional wave. (c): time evolution of p recorded at receiver (x = -60, y = 60) from t = 0.028 s to the final instant. (d) zoom on the time evolution of scattered waves.

Figure 11 :Figure 12 :

 1112 Figure 11: Test 4, viscous fluids. Snapshot of p (a) and σ 11 (b) after 180 iterations, with a refinement factor q = 15 and a third order interface method (r = 3). Snapshot of p around the scatterer (c), (d). Zoom of p on the line y = 0 around the reflected fast compressional wave (e) and around the diffusive slow wave around the right interface (f).

Figure 13 :

 13 Figure 13: Test 5. Snapshots of p at two different times t 1 ant t 2 , corresponding to 140 and 280 iterations. Left column: inviscid fluids; right column: viscous saturating fluids.

( 1 -) k pf ω 2 µ sin θ cos θ k pf ω λ f + 2 µ

 12 β m = 2 m-1 , ω 0 = 2π f 0 ; the coefficients a m are: a 1 = 1, a 2 = -21/32, a 3 = 63/768, a 4 = -1/512, ensuring C 6 smoothness. The support of the incident plane wave lies initially in Ω 0 . If η = 0, this wave is slightly dispersive and its time-domain expression follows from a Fourier synthesis:Û(r, ω) = γ Y pf ) cos θ φ (1 -Y pf ) sin θ k pf ω λ f + 2 µ cos 2 θ + β m φ (Y pf -1sin 2 θ + β m φ (Y pf -1)k pf ω m (β + φ (Y pf -1))

2 if

 2 r = (x -x s ) 2 + (y -y s ) 2 ≤ r 0 ,

Table 2 :

 2 Physical parameters of the matrix (Ω 0 ) and of the scatterer (Ω 1 ), corresponding to sandstone and shale saturated with water, respectively.

		2650	2211
	µ (Pa)	2.926 10 9	3.539 10 9
	ρ f (kg/m 3 ) 1040	1040
	η (Pa.s)	1.5 10 -3	10 -3
	φ	0.335	0.05
	a	2	2
	κ (m 2 )	10 -11	5. 10 -12
	λ f (Pa)	6.1425 10 9 4.689 10 9
	β	0.9558	0.0527
	m (Pa)	6.491 10 9	9.852 10 9
	c pf (m/s)	2384.1	2350.4
	c ps (m/s)	758.9	486.4
	c s (m/s)	1229.0	1290.0
	f c (Hz)	3844.9	765.1

  This source generates cylindrical waves of all types: fast and slow compressional waves, shear waves.

	<	2 f 0	,
	0 otherwise.
	(C.6)

Appendix A. Matrices involved in system [START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF] Based on ( 6) and [START_REF] Carcione | Some aspects of the physics and numerical modeling of Biot compressional waves[END_REF], the matrices in [START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF] write

Appendix B. Algorithm for the Beltrami-Michell conditions

The following algorithm is proposed to compute the non-zero components of matrices G r k (k = 0, 1) involved in the Step 2 of the immersed interface method (see section 3.3):

for γ = 0, ..., r, for δ = 0, ..., γ if δ = 0 then for ε = 1, ..., 8

Appendix C. Implementation of sources

Two sources are considered. The first one involves a plane right-going fast compressional wave, whose wavevector k makes an angle θ with the