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

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

Guillaume Chiavassa guillaume.chiavassa@centrale-marseille.fr Bruno Lombard lombard@lma.cnrs-mrs.fr M2P2, UMR 6181 - CNRS - Ecole Centrale Marseille, Technopôle de Chateau-Gombert, 38 rue Frédéric Joliot-Curie, 13451 Marseille, France Laboratoire de Mécanique et d’Acoustique, UPR 7051 CNRS, 31 chemin Joseph Aiguier, 13402 Marseille, France

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.

porous media, elastic waves, Biot’s model, time splitting, finite difference methods, Cartesian grid, immersed interface method
35L50, 65M06
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 , 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 , 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.

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 BIOT56-A (); 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:


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


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


where is the Lamé coefficient of the dry matrix.

On the other hand, the conservation of momentum yields


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 low-frequency range, involving frequencies lower than


If , more sophisticated models are required BIOT56-B (); LU05 (). In practice, the viscosity of the fluid is always non-zero; 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 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


where , , and are force densities, , , and . Setting


equations (6) are written as a first-order non-homogeneous linear system


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


is an energy that satisfies


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

2.3 Heterogeneous media

Figure 1: Interface between two poroelastic media and


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


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 ():


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)
Figure 2: Phase velocities (a) and attenuations (b) of the solutions to Biot’s model corresponding to the porous medium in table 2. : fast compressional wave; : slow compressional wave; : shear wave. In (a), the horizontal dotted lines refer to the eigenvalues , and of and .

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:


where the roots and satisfy . If is orthogonal with , the dispersion relation of the shear wave is obtained:


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 non-dispersive and non-dissipative, 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 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 , and the time step is . A straightforward discretization of (8) by an explicit time scheme leads to the following stability condition:


where is obtained by a Von-Neumann 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 second-order Strang’s splitting LEVEQUE07 (), solving alternatively the hyperbolic system


and then the diffusive system with a source term


The linear system (16) is solved by applying any scheme for hyperbolic systems, giving . 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 STRIKWERDA99 (); 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: . In this case, and are unchanged, whereas the velocities become ()


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

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:


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 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: high-order 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:


where and . Following this formalism, the zero-th order jump conditions (12) are written


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


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


where and are the matrices in . The jump condition (21) is also differentiated in terms of . Taking advantage of the chain-rule, we obtain e.g.


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


where are matrices (), and . The computation of matrices is a tedious task when , 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 and , as many times as required:


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


where are matrices, and if , otherwise. Based on (20) and (26), an algorithm to compute the non-zero components of is proposed in B.

Step 3: high-order boundary values.
Based on (25) and (27), the vectors of independent boundary values satisfy


where are matrices. Since the system (28) is underdetermined, the solution is not unique, and hence it can be written


where is the least-squares pseudo-inverse 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 ().

Figure 3: Irregular point and its orthogonal projection onto . The grid nodes used to compute are inside the circle with radius and centered on ; they are denoted by .

Step 4: construction of modified values.
Let be the matrix of -th order 2D Taylor expansions


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:


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


At the grid points of , -th order Taylor expansion of the solution at and the boundary conditions (29) give


Equations (32) and (33) are written in the matrix form


where is a convenient matrix. To ensure that the system (34) is overdetermined, the radius of the disc is chosen to satisfy


Exact values in (34) are replaced by numerical ones, and the Taylor rests are removed. The least-squares 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


Lastly, the modified value follows from (31) and (36):

Quantity Size
if , 0 else
Table 1: Quantities involved in the computation of the modified values (section 3.3).

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 in (34) depends on the subcell position of inside the mesh and on the jump conditions at , involving the local geometry and the curvature of at . Consequently, all these insights are incorporated in the modified value (37), and hence in the scheme.

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

  4. 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 fourth-order ADER scheme () requires a third-order immersed interface method () to maintain fourth-order convergence.

  5. 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:


  • Detection of irregular grid points

  • Computation of extrapolation matrices in (37)

  • Initialization of the solution at

  • Diffusive step (18) where on the coarse grid

  • Diffusive step (18) where on the refined grids

Time iterations

  • Coarse grid:

    • Computation of modified values (37) if present

    • Solving the propagative step (16)

    • Diffusive step (18) where

  • On each refined grid, subtime iterations:

    • Space-time interpolations at the grid boundaries

    • Computation of modified values (37)

    • Solving the propagative step (16)

    • Diffusive step (18) where

End of time iterations

  • Diffusive step (18) where on the coarse grid

  • Diffusive step (18) where on the refined grids

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

(kg/m) 2650 2211
(kg/m) 1040 1040
0.335 0.05
2 2
(m/s) 2384.1 2350.4
(m/s) 758.9 486.4
(m/s) 1229.0 1290.0
(Hz) 3844.9 765.1
Table 2: Physical parameters of the matrix () and of the scatterer (), corresponding to sandstone and shale saturated with water, respectively.

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

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
Table 3: Test 1. Convergence rate in norm. No immersed interface method , linear , quadratic or cubic immersed interface method.
Figure 5: Test 1. Error measured in norm versus the number of points , with various order of the immersed interface method. Dotted line corresponds to 4-th order slope.

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, first-order accuracy is obtained. As stated in section 3.3, fourth-order accuracy is maintained if , i.e if third-order 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.

4.3 Test 2: mesh refinement

Figure 6: Test 2. Snapshots of , and of the fast (), slow (), and shear () waves emitted by a source point. Left column: viscosity . Right column: viscous case. Dashed areas indicate the location of the refined grids and