Abstract

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.

keywords:

Msc:

35L50, 65M06

Pacs:

43.20.-Gp, 46.40.-f

††journal: Journal of Computational Physics

1 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 fc, 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 BIOT56-A (). An adequate modeling of this diffusive mode remains a major challenge in real applications. At frequencies greater than fc, 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 BIOT56-B (); MASSON06 (). Experimental observations of the slow wave in the low-frequency range PLONA80 () and in the high-frequency range CHANDLER81 () 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 DAI95 (); ZENG01 (); WENZLAU09 (), finite elements ZHAO05 (), discontinuous Galerkin methods PUENTE08 (), boundary elements SCHANZ01 (), pseudospectral methods CARCIONE07 () and spectral element methods MORENCY08 (). A recent review of computational poroelasticity can be found in CARCIONE10 (). 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 SCHWARTZKOPFF04 () to integrate Biot’s equations. A flux-conserving space-time mesh refinement BERGER98 () is implemented around the interfaces to capture the slow compressional wave. Lastly, an immersed interface method LOMBARD04 (); LOMBARD06 () 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 CHIAVASSA10 (), 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.

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 us and uf, the elastic stress tensor σ, and the acoustic pressure p. In one hand, the constitutive laws are:

(1)

where I is the identity, ξ is the rate of fluid change, and ε is the strain tensor

ξ=−∇.(ϕ(uf−us)),ε=12(∇us+∇uTs).

(2)

The symmetry of σ in (1) implies compatibility conditions between spatial derivatives of ε, leading to the Beltrami-Michell equation RICE76 (); COUSSY95 ()

where vs=∂us∂t=(vs1,vs2)T is the elastic velocity, and w=ϕ∂∂t(uf−us)=(w1,w2)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

fc=ηϕ2πaκρf.

(5)

If f≥fc, more sophisticated models are required BIOT56-B (); LU05 (). In practice, the viscosity of the fluid is always non-zero; nevertheless, considering η=0 can be relevant for two reasons:

if f≫fc, the viscous forces are smaller than the inertial forces FENG83a (); MORENCY08 () 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.

2.2 Evolution equations

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

The energy of poroelastic waves can be deduced from (6). Without any source terms (F=0) and setting Cε=σ+βpI, it is proven in THESE_EZZIANI () that

E(t)=12∫R2(ρvs2+Cε:ε)dS+12∫R2(ρww2+1mp2)dS+∫R2ρfvs.wdS

(9)

is an energy that satisfies

dEdt=−∫R2ηκw2dS.

(10)

Consequently, E is conserved when the viscous effects are neglected (η=0) and is a decreasing function otherwise.

2.3 Heterogeneous media

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

.

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′=dxdτ and y′=dydτ 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 GUREVICH99 ():

[vs]=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.

The eigenvalues of A and B in (8) are real: ±¯¯cpf, ±¯¯cps, ±¯¯cs, and 0 (multiplicity 2), where ¯¯cpf>max(¯¯cs,¯¯cps)>0. If η≠0, the spectral radius R(S)=ηκρχ can be very large and then the system (8) is stiff.

A plane wave dei(ωt−k.r) is injected in (6), where k=ke 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:

where the roots are ±ks, Re{ks}>0. Based on (13) and (14), the phase velocities c=ω/Re{k} and the attenuations α=Im{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 cpf(f), cps(f) and cs(f) are monotonically increasing functions, tending asymptotically towards the eigenvalues ¯¯cpf, ¯¯cps and ¯¯cs. 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≪fc, then cps(f)≪¯¯cps, and the slow compressional wave tends towards a static diffusive mode CHANDLER81 (). At greater frequencies, cps 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 BOURBIE87 (). 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.

3 Numerical modeling

3.1 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¯¯cpf,2R(S)),

(15)

where Θ is obtained by a Von-Neumann analysis of stability when S=0. In (15), 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 ¯¯cpfΔt/Δx≈10−12≪Θ, which is intractable for computations.

We follow here a more efficient strategy based on second-order Strang’s splitting LEVEQUE07 (), solving alternatively the hyperbolic system

∂∂tU+A∂∂xU+B∂∂yU=0,

(16)

and then the diffusive system with a source term

∂∂tU=−SU+F.

(17)

The linear system (16) is solved by applying any scheme for hyperbolic systems, giving Un+1/2i,j. In the numerical experiments performed in section 4, a fourth-order ADER scheme is used SCHWARTZKOPFF04 (), which involves a centered stencil of 25 nodes. On Cartesian grids, this scheme amounts to a fourth-order Lax-Wendroff scheme LORCHER05 (). It is dispersive of order 4 and dissipative of order 6, and its stability limit is Θ=1STRIKWERDA99 (); HDR-LOMBARD (). 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)

vn+1k=vn+1/2k+ρfρ(1−e−ηκρχT)wn+1/2k,wn+1k=e−ηκρχTwn+1/2k,

(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: ¯¯cpfΔ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 PUENTE08 (), 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 CHIAVASSA10 () 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.

3.2 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 ¯¯cpf≫¯¯cps 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 BERGER84 (); BERGER98 (), which is more naturally coupled to the flux-conserving scheme developed to solve (16). 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 RODRIGUEZ05 (); THESE_RODRIGUEZ (). 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 HDR-LOMBARD (). In the case of the Lax-Wendroff scheme applied to the scalar advection equation, the stability of the coupling is proven in BERGER85 () 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(f0)=cpf(f0)/f0 and λps(f0)=cps(f0)/f0 are deduced from the dispersion analysis, where f0 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(f0)Δx/q≈λpf(f0)Δx⇒q≈cpf(f0)cps(f0).

(19)

3.3 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 (12) are not enforced numerically if no special treatment is applied. Lastly, the smoothness requirements to solve (16) 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 PIRAUX01 (); LOMBARD04 (); LOMBARD06 (); HDR-LOMBARD (). 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 (12). The
parameter r is discussed at the end of this section.

Let us consider a point M(xI,yJ)∈Ω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 Urk with nv=4(r+1)(r+2) components:

Urk=limM→P,M∈Ωk(UT,...,∂l∂xl−m∂ymUT,...,∂r∂yrUT)T,

(20)

where l=0,...,r and m=0,...,l. Following this formalism, the zero-th order jump conditions (12) are written

C01U01=C00U00,

(21)

where the matrices of the jump conditions C0k depend on the local geometry of Γ:

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 (16). For example, we obtain

∂∂t(C00U00)=−C00A0∂∂xU00−C00B0∂∂yU00,

(23)

where A0 and B0 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.

ddτ(C00U00)=(ddτC00)U00+C00(x′∂∂xU00+y′∂∂yU00).

(24)

From (21), (23) and (24), we build matrices C1k such that C11U11=C10U10, which provides first-order jump conditions. By iterating this process r times, r-th order interface conditions are obtained

Cr1Ur1=Cr0Ur0,

(25)

where Crk are nc×nv matrices (k=0,1), and nc=3(r+1)(r+2). The computation of matrices Crk 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:

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 Urk. For this purpose, we define the vectors Vrk such that

Urk=GrkVrk,

(27)

where Grk are nv×(nv−nb) matrices, and nb=r(r−1)/2 if r≥2, nb=0 otherwise. Based on (20) and (26), an algorithm to compute the non-zero components of Grk is proposed in B.

where Srk=CrkGrk are nc×(nv−nb) matrices. Since the system (28) is underdetermined, the solution is not unique, and hence it can be written

Vr1=((Sr1)−1Sr0|KSr1)(Vr0Λr),

(29)

where (Sr1)−1 is the least-squares pseudo-inverse of Sr1, KSr1 is the matrix filled with the kernel of Sr1, and Λr is a set of nv−nc−nb Lagrange multipliers that represent the coordinates of Vr1 onto the kernel. A singular value decomposition of Sr1 is used to build (Sr1)−1 and the kernel KSr1NRC ().

Figure 3: Irregular point M(xI,yJ)∈Ω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 +.

where I8 is the 8×8 identity matrix, l=0,...,r and m=0,...,l. The modified value at (xI,yJ) is a smooth extension of the solution on the other side of Γ (figure 3), and it writes:

U∗I,J=ΠrI,JUr0=ΠrI,JGr0Vr0.

(31)

The vector Vr0 in (31) remains to be estimated in terms of the boundary conditions and of the numerical values at surrounding grid points. For this purpose, we consider the disc D centered at P with a radius d, that contains Nd grid points. At the grid points of D∩Ω0, r-th order Taylor expansion of the solution at P gives

U(xi,yj,tn)=Πri,jUr0+O(Δxr+1),=Πri,jGr0(1|0)(Vr0Λr)+O(Δxr+1).

(32)

At the grid points of D∩Ω1, r-th order Taylor expansion of the solution at P and the boundary conditions (29) give

where M is a convenient 8Nd×(2nv−2nb−nc) matrix. To ensure that the system (34) is overdetermined, the radius d of the disc is chosen to satisfy

ε(d,r)=8Nd2nv−2nb−nc≥1.

(35)

Exact values in (34) 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 (31). As a consequence, they can be removed and the (nv−nb)×8Nd restriction ¯¯¯¯¯¯¯¯¯¯¯M−1 of M−1 is defined by

Table 1: Quantities involved in the computation of the modified values (section 3.3).

Comments and practical details.

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.

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 (37), and hence in the scheme.

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 (35), 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 Nd≈20 and ε≈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 s-th order overall accuracy GUSTAFSSON75 (). As a consequence, a fourth-order ADER scheme (s=4) requires a third-order immersed interface method (r=3) to maintain fourth-order convergence.

A GKS analysis of stability has been performed in 1D in the case of an inviscid saturating fluid HDR-LOMBARD (). 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.

3.4 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:

4 Numerical experiments

4.1 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 DAI95 (), 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.

Parameters

Ω0

Ω1

ρs (kg/m3)

2650

2211

μ (Pa)

2.926109

3.539109

ρf (kg/m3)

1040

1040

η (Pa.s)

1.510−3

10−3

ϕ

0.335

0.05

a

2

2

κ (m2)

10−11

5.10−12

λf (Pa)

6.1425109

4.689109

β

0.9558

0.0527

m (Pa)

6.491109

9.852109

¯¯cpf (m/s)

2384.1

2350.4

¯¯cps (m/s)

758.9

486.4

¯¯cs (m/s)

1229.0

1290.0

fc (Hz)

3844.9

765.1

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

Once the spatial mesh sizes Δx and Δy are chosen on the coarse grid, the time step follows from the CFL number in Ω0: ¯¯cpf0Δ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.

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 C. In the case of an incident plane wave, the exact expression given in 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 MARTIN08 (); ZENG01 (). 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.

4.2 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] m2 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 (41), with parameters: γ=10−3, f0=40 Hz, t=3.310−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.

(a)

(b)

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.

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 l2 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.