Time domain numerical modeling of wave propagation in 2D heterogeneous porous media
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 lowfrequency 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 fourthorder ADER scheme, while the diffusive part is solved analytically. Near the material interfaces, a spacetime 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:
porous media, elastic waves, Biot’s model, time splitting, finite difference methods, Cartesian grid, immersed interface methodMsc:
35L50, 65M06Pacs:
43.20.Gp, 46.40.f1 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 , 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 BIOT56A (). An adequate modeling of this diffusive mode remains a major challenge in real applications. At frequencies greater than , 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 BIOT56B (); MASSON06 (). Experimental observations of the slow wave in the lowfrequency range PLONA80 () and in the highfrequency range CHANDLER81 () have confirmed Biot’s theory. In the current study, we focus on the lowfrequency range.
Until the 1990’s, Biot’s equations were mainly studied in the harmonic regime. Various timedomain 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 lowfrequency 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 wellsuited 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 timesplitting is used along with a fourthorder ADER scheme SCHWARTZKOPFF04 () to integrate Biot’s equations. A fluxconserving spacetime 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 onedimensional 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.
2 Physical modeling
2.1 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 BIOT56A (); BOURBIE87 (); CARCIONE07 (); CARCIONE10 (). 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 and the dynamic viscosity of the fluid; the density and the shear modulus of the elastic skeleton; the porosity , the tortuosity , the absolute permeability , the Lamé coefficient and the two Biot’s coefficients and of the saturated matrix. The unknowns are the elastic and acoustic displacements and , the elastic stress tensor , and the acoustic pressure . In one hand, the constitutive laws are:
(1) 
where is the identity, is the rate of fluid change, and is the strain tensor
(2) 
The symmetry of in (1) implies compatibility conditions between spatial derivatives of , leading to the BeltramiMichell equation RICE76 (); COUSSY95 ()
(3) 
where is the Lamé coefficient of the dry matrix.
On the other hand, the conservation of momentum yields
(4) 
where is the elastic velocity, and is the filtration velocity. To be valid, the second equation of (4) requires that the spectrum of the waves lies mainly in the lowfrequency range, involving frequencies lower than
(5) 
If , more sophisticated models are required BIOT56B (); LU05 (). In practice, the viscosity of the fluid is always nonzero; nevertheless, considering can be relevant for two reasons:

if , the viscous forces are smaller than the inertial forces FENG83a (); MORENCY08 () and can be neglected to a first approximation;

the exact solutions of poroelastodynamic equations are computed more accurately if the saturating fluid is inviscid, which is attractive to validate the numerical methods.
2.2 Evolution equations
A velocitystress formulation is followed: from (1) and (4), we obtain the system of PDEs
(6) 
where , , and are force densities, , , and . Setting
(7) 
equations (6) are written as a firstorder nonhomogeneous linear system
(8) 
where , and are real matrices detailled in A.
The energy of poroelastic waves can be deduced from (6). Without any source terms () and setting , it is proven in THESE_EZZIANI () that
(9) 
is an energy that satisfies
(10) 
Consequently, is conserved when the viscous effects are neglected () and is a decreasing function otherwise.
2.3 Heterogeneous media
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 and , which are separated by a stationary interface described by a parametric equation (figure 1). At any point on , the unit tangential vector and the unit normal vector are
(11) 
The derivatives and 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 ():
(12) 
Enforcing these conditions is one of the main objective of the immersed interface method presented in section 3.3.
2.4 Dispersion analysis
(a)  (b) 

The eigenvalues of and in (8) are real: , , , and 0 (multiplicity 2), where . If , the spectral radius can be very large and then the system (8) is stiff.
A plane wave is injected in (6), where and are the wavevector and the polarization, respectively; is the position, is the angular frequency and is the frequency. If is collinear with , the dispersion relation of compressional waves is obtained:
(13) 
where the roots and satisfy . If is orthogonal with , the dispersion relation of the shear wave is obtained:
(14) 
where the roots are , . Based on (13) and (14), the phase velocities and the attenuations of each wave are defined. In the remainder of this article, the subscripts , and denote fast compressional, slow compressional and shear waves, respectively.
The phase velocities , and are monotonically increasing functions, tending asymptotically towards the eigenvalues , and . If , the three waves are nondispersive and nondissipative, and the energy of poroelastic waves (9) is conserved. If , 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 , then , and the slow compressional wave tends towards a static diffusive mode CHANDLER81 (). At greater frequencies, is larger but the attenuation increases. These properties are summarized in figure 2.
In the lowfrequency 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 , and the time step is . A straightforward discretization of (8) by an explicit time scheme leads to the following stability condition:
(15) 
where is obtained by a VonNeumann analysis of stability when . In (15), the bound induced by the spectral radius of can be very restrictive. In sandstone saturated with bitumen, for example, the maximal CFL number is roughly , which is intractable for computations.
We follow here a more efficient strategy based on secondorder Strang’s splitting LEVEQUE07 (), solving alternatively the hyperbolic system
(16) 
and then the diffusive system with a source term
(17) 
The linear system (16) is solved by applying any scheme for hyperbolic systems, giving . In the numerical experiments performed in section 4, a fourthorder ADER scheme is used SCHWARTZKOPFF04 (), which involves a centered stencil of 25 nodes. On Cartesian grids, this scheme amounts to a fourthorder LaxWendroff scheme LORCHER05 (). It is dispersive of order 4 and dissipative of order 6, and its stability limit is STRIKWERDA99 (); HDRLOMBARD (). Other singlegrid 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: . In this case, and are unchanged, whereas the velocities become ()
(18) 
where depends on the time step (see section 3.4). The splitting (16)(17) along with exact integration (18) recovers the optimal condition of stability: .
Since the matrices and do not commute with , the theoretical order of convergence falls from 4 to 2 when the viscosity is nonnegligible. Using a fourthorder accurate scheme such as ADER 4 is nevertheless advantageous, compared with a secondorder scheme such as LaxWendroff: the stability limit is improved, and numerical artifacts (dispersion, attenuation, anisotropy) are greatly reduced.
In PUENTE08 (), the authors notice that the firstorder splitting does not lead to a correct representation of the slow mode at low frequencies. Nevertheless, the numerous onedimensional examples provided in CHIAVASSA10 () demonstrate that the secondorder 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 spacetime 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 () the property 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 spacetime mesh refinement approach based on flux conservation BERGER84 (); BERGER98 (), which is more naturally coupled to the fluxconserving scheme developed to solve (16). The refined zones are rectangular Cartesian patches with mesh sizes , where the integer is the refinement factor. To reduce the cpu time and to limit the numerical dispersion on the coarse grid, a local time step is used RODRIGUEZ05 (); THESE_RODRIGUEZ (). When one time step is done on the coarse grid, 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 HDRLOMBARD (). In the case of the LaxWendroff scheme applied to the scalar advection equation, the stability of the coupling is proven in BERGER85 () whatever .
The additional cost induced by mesh refinement can become prohibitive, both concerning the memory requirements and the computational time, because of the substeps inside one time step. The value of must be estimated carefully in terms of the physical parameters. For this purpose, the wavelengths and are deduced from the dispersion analysis, where 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:
(19) 
3.3 Immersed interface method
The discretization of the interfaces requires special care. A straightforward stairstep discretization of the interfaces introduces a firstorder 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 (); HDRLOMBARD (). 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 successive derivatives of the jump conditions (12). The parameter is discussed at the end of this section.
Let us consider a point and its orthogonal projection onto (figure 3). The algorithm to build the modified value at is divided into four steps.
Step 1: highorder interface conditions.
On the side (), the boundary values of the spatial derivatives of up to the th order are put in a vector with components:
(20) 
where and . Following this formalism, the zeroth order jump conditions (12) are written
(21) 
where the matrices of the jump conditions depend on the local geometry of :
(22) 
The jump condition (21) is differentiated with respect to time , and then the time derivatives are replaced by spatial derivatives thanks to the conservation law (16). For example, we obtain
(23) 
where and are the matrices in . The jump condition (21) is also differentiated in terms of . Taking advantage of the chainrule, we obtain e.g.
(24) 
From (21), (23) and (24), we build matrices such that , which provides firstorder jump conditions. By iterating this process times, th order interface conditions are obtained
(25) 
where are matrices (), and . The computation of matrices is a tedious task when , that can be greatly simplified using computer algebra tools.
Step 2: highorder BeltramiMichell equations.
The equation (3) is satisfied anywhere in a poroelastic medium. Under sufficient smoothness requirements, it can be differentiated with respect to and , as many times as required:
(26) 
where and . The equations (26) are also satisfied along . They can be used therefore to reduce the number of independent components in . For this purpose, we define the vectors such that
(27) 
where are matrices, and if , otherwise. Based on (20) and (26), an algorithm to compute the nonzero components of is proposed in B.
Step 3: highorder boundary values.
Based on (25) and (27), the vectors of independent boundary values satisfy
(28) 
where are matrices. Since the system (28) is underdetermined, the solution is not unique, and hence it can be written
(29) 
where is the leastsquares pseudoinverse of , is the matrix filled with the kernel of , and is a set of Lagrange multipliers that represent the coordinates of onto the kernel. A singular value decomposition of is used to build and the kernel NRC ().
Step 4: construction of modified values.
Let be the matrix of th order 2D Taylor expansions
(30) 
where is the identity matrix, and . The modified value at is a smooth extension of the solution on the other side of (figure 3), and it writes:
(31) 
The vector 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 centered at with a radius , that contains grid points. At the grid points of , th order Taylor expansion of the solution at gives
(32) 
At the grid points of , th order Taylor expansion of the solution at and the boundary conditions (29) give
(33) 
Equations (32) and (33) are written in the matrix form
(34) 
where is a convenient matrix. To ensure that the system (34) is overdetermined, the radius of the disc is chosen to satisfy
(35) 
Exact values in (34) are replaced by numerical ones, and the Taylor rests are removed. The leastsquares inverse of is denoted by . The Lagrange multipliers are accounted in the construction of , but are not involved in the definition of the modified value (31). As a consequence, they can be removed and the restriction of is defined by
(36) 
Lastly, the modified value follows from (31) and (36):
(37) 
Quantity  Size 

if , 0 else  
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 matrixvector 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 timemarching.

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 depending on the local geometry of at . We adopt here a simpler strategy, based on a constant radius . Taking , numerical experiments have shown that is a good candidate, while is used when . In this case, we obtain typically and .

The order plays an important role on the accuracy of the coupling between the immersed interface method and a th order scheme. If , then a th order local truncation error is obtained at the irregular points. This condition can be slightly relaxed: still ensures a th order overall accuracy GUSTAFSSON75 (). As a consequence, a fourthorder ADER scheme () requires a thirdorder immersed interface method () to maintain fourthorder convergence.

A GKS analysis of stability has been performed in 1D in the case of an inviscid saturating fluid HDRLOMBARD (). 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 BeltramiMichell 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 spacetime 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 spacetime 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: Pa.s, the other parameters being unchanged. As recalled in section 2.1, this limitcase has physical significance only in the highfrequency range. It is mainly addressed here for a numerical purpose.
Parameters  

(kg/m)  2650  2211 
(Pa)  
(kg/m)  1040  1040 
(Pa.s)  
0.335  0.05  
2  2  
(m)  
(Pa)  
(Pa)  
(m/s)  2384.1  2350.4 
(m/s)  758.9  486.4 
(m/s)  1229.0  1290.0 
(Hz)  3844.9  765.1 
Once the spatial mesh sizes and are chosen on the coarse grid, the time step follows from the CFL number in : . If , 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 PerfectlyMatched 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 m 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: , Hz, s, and degrees (figure 4a). Consequently, the incident wave propagates normally to the interface, leading to a 1D configuration; from a numerical point of view, however, the problem is fully bidimensional.
(a)  (b) 

The computations are done on a uniform grid of points, during time steps. Comparisons with the exact values of the pressure are done on the subdomain m, in order to avoid spurious effects induced by the edges of the computational domain (figure 4b). 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.
order  order  order  order  

400          
800  1.973  1.961  2.520  2.552  
1200  1.599  2.242  3.809  3.690  
1600  1.086  2.173  3.632  3.619  
2000  0.998  2.119  3.906  3.861  
2400  1.093  2.201  4.168  3.917  
2800  0.964  2.069  3.473  3.945  
3200  0.902  1.971  3.685  4.013  
3600  1.038  2.114  3.718  3.915  
4000  0.879  1.933  3.552  3.986 
Errors in norm and convergence rates are reported in table 3 and drawn on figure 5. Various values of are investigated; means that no immersed interface method is applied; in this case, firstorder accuracy is obtained. As stated in section 3.3, fourthorder accuracy is maintained if , i.e if thirdorder extrapolations are used in the immersed interface method. In the present test case, 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 in the following simulations.