High-Resolution Finite Volume Modeling of Wave Propagation in Orthotropic Poroelastic Media

High-Resolution Finite Volume Modeling of Wave Propagation in Orthotropic Poroelastic Media

Grady I. Lemoine111Department of Applied Mathematics, University of Washington, Guggenheim Hall Box 352420, Seattle, WA 98195. Supported in part by NIH grant 5R01AR53652-2 and NSF grant DMS-0914942 333Corresponding author, email: gl@uw.edu    M. Yvonne Ou222Department of Mathematical Sciences, University of Delaware, 501 Ewing Hall, Newark, DE 19716. Supported in part by NSF-DMS Mathematical Biology Grant 0920852 and NIH CBER PILOT Grant 322 159    Randall J. LeVeque111Department of Applied Mathematics, University of Washington, Guggenheim Hall Box 352420, Seattle, WA 98195. Supported in part by NIH grant 5R01AR53652-2 and NSF grant DMS-0914942

Poroelasticity theory models the dynamics of porous, fluid-saturated media. It was pioneered by Maurice Biot in the 1930s through 1960s, and has applications in several fields, including geophysics and modeling of in vivo bone. A wide variety of methods have been used to model poroelasticity, including finite difference, finite element, pseudospectral, and discontinuous Galerkin methods. In this work we use a Cartesian-grid high-resolution finite volume method to numerically solve Biot’s equations in the time domain for orthotropic materials, with the stiff relaxation source term in the equations incorporated using operator splitting. This class of finite volume method has several useful properties, including the ability to use wave limiters to reduce numerical artifacts in the solution, ease of incorporating material inhomogeneities, low memory overhead, and an explicit time-stepping approach. To the authors’ knowledge, this is the first use of high-resolution finite volume methods to model poroelasticity. The solution code uses the clawpack finite volume method software, which also includes block-structured adaptive mesh refinement in its amrclaw variant. We present convergence results for known analytic plane wave solutions, achieving second-order convergence rates outside of the stiff regime of the system. Our convergence rates are degraded in the stiff regime, but we still achieve similar levels of error on the finest grids examined. We also demonstrate good agreement against other numerical results from the literature. To aid in reproducibility, we provide all of the code used to generate the results of this paper, at https://bitbucket.org/grady_lemoine/poro-2d-cartesian-archive.


oroelastic, wave propagation, finite-volume, high-resolution, running head, stiff relaxation, operator splitting

1 Introduction

Poroelasticity theory is a homogenized model for solid porous media containing fluids that can flow through the pore structure. This field was pioneered by Maurice A. Biot, who developed his theory of poroelasticity from the 1930s through the 1960s; a summary of much of Biot’s work can be found in his 1956 and 1962 papers [4, 5, 6]. Biot theory uses linear elasticity to describe the solid portion of the medium (often termed the skeleton or matrix), linearized compressible fluid dynamics to describe the fluid portion, and Darcy’s law to model the aggregate motion of the fluid through the matrix. While it was originally developed to model fluid-saturated rock and soil, Biot theory has also been used in underwater acoustics [10, 30, 31], and to describe wave propagation in in vivo bone [19, 20, 29].

Biot theory predicts rich and complex wave phenomena within poroelastic meterials. Three different types of waves appear: fast P waves analogous to standard elastic P waves, in which the fluid and matrix show little relative motion, and typically compress or expand in phase with each other; shear waves analogous to elastic S waves; and slow P waves, where the fluid expands while the solid contracts, or vice versa. The slow P waves exhibit substantial relative motion between the solid and fluid compared to waves of the other two types. The viscosity of the fluid dissipates poroelastic waves as they propagate through the medium, with the fast P and S waves being lightly damped and the slow P wave strongly damped. The viscous dissipation also causes slight dispersion in the fast P and S waves, and strong dispersion in the slow P wave.

A variety of different numerical approaches have been used to model poroelasticity. Carcione, Morency, and Santos provide a thorough review of the previous literature [14]. The earliest numerical work in poroelasticity seems to be that of Garg [27], using a finite difference method in 1D. Finite difference and pseudospectral methods have continued to be popular since then, with further work by Mikhailenko [39], Hassanzadeh [33], Dai et al. [22], and more recently Chiavassa and Lombard [17], among others. Finite element approaches began being used in the 1980s, with Santos and Oreña’s work [43] being one of the first. Boundary element methods have also been used, such as in the work of Attenborough, Berry, and Chen [2]. Spectral element methods have also been used in both the frequency domain [24] and the time domain [40]. With the recent rise of discontinuous Galerkin methods, DG has been applied to poroelasticity in several works, such as that of de la Puente et al. [23]. There have also been semi-analytical approaches to solving the poroelasticity equations, such as that of Detournay and Cheng [26], who analytically obtain a solution in the Laplace transform domain, but are forced to use an approximate inversion procedure to return to the time domain. Finally, there has been significant work on inverse problems in poroelasticity, for which various forward solvers have been used; of particular note is the paper of Buchanan, Gilbert, and Khashanah [9], who used the finite element method (specifically the FEMLAB software package) to obtain time-harmonic solutions for cancellous bone as part of an inversion scheme to estimate poroelastic material parameters, and the later papers of Buchanan and Gilbert [7, 8], where the authors instead used numerical contour integration of the Green’s function. Numerical work in the 1970s and 1980s focused on isotropic poroelasticity, with the earliest work on anisotropic poroelasticity being by Carcione in 1996 [12].

A major theme in time-domain numerical modeling of poroelasticity has been the difficulty of handling the viscous dissipation term, which has its own intrinsic time scale and causes the poroelasticity system to be stiff, at least if low-frequency waves are being considered. (The time scales associated with dissipation are independent of frequency, so at higher frequencies there is less separation between them and the time scale associated with wave motion — in other words, the system is not stiff if sufficiently high-frequency waves are being considered.) This viscous dissipation is particularly problematic for the slow P wave. While it can still be addressed if Biot’s equations are solved as a unified system (e.g. with an implicit time-integration method), since the viscous dissipation term is easy to solve analytically in isolation, operator splitting approaches have also been a popular way to address this issue; we adopt this approach as well. Carcione and Quiroga-Goode used operator splitting in conjunction with a pseudospectral method [15], and Chiavassa and Lombard use it with a finite difference approach [17]. De la Puente et al. also investigate an operator splitting approach, and encounter difficulties in obtaining a fast rate of convergence due to the stiffness of the dissipation term, which pushed the problems they investigated toward the diffusive limit of the Biot system [23].

Our work in this paper solves a velocity-stress formulation of Biot linear orthotropic poroelasticity theory using Cartesian-grid high-resolution finite volume methods. These methods are memory-efficient explicit techniques designed to model hyperbolic systems, and can include wave limiters designed to reduce the effect of numerical artifacts in the solution. Since they are based on solving Riemann problems at grid interfaces, it is also straightforward to include material inhomogeneities between cells using these methods. To our knowledge, this is the first use of finite volume methods to model poroelasticity. We employ the clawpack finite volume method package [45], which offers the use of operator splitting to model viscous dissipation, and includes optional Berger-Colella-Oliger block-structured adaptive mesh refinement [3] to improve solution efficiency for larger problems. We verify our code both against known analytical solutions and against other numerical solutions from the poroelasticity literature.

2 Poroelasticity theory for transversely isotropic materials

Biot’s equations of poroelasticity are a complicated system of PDEs that exhibit rich and varied behaviors. The reader is encouraged to refer to a detailed treatment of the subject, such as chapter 7 of Carcione’s book [13], but we also provide an overview here. After a review of the basic relations modeling the behavior of a poroelastic medium in sections 2.1 through 2.4, section 2.5 presents the linear first-order system of PDEs that forms the basis of our numerical model. Section 2.6 then defines an energy norm for our state vector that solves some of the scaling issues associated with modeling poroelasticity in SI units. The matrix created for this energy norm allows us to explore some useful properties of the model: section 2.7 provides a concise proof that our first-order system is hyperbolic, and section 2.8 exhibits a strictly convex entropy function for the system, which has implications for the correctness of our numerical solutions that we explore later in section 3.3.

2.1 Stress-strain relation

We assume the constituent material of the solid matrix is isotropic and the anisotropy of the solid matrix results purely from the microstructure. We also assume that the anisotropy has a specific form — that the medium is orthotropic, possessing three orthogonal planes of symmetry, and is isotropic with respect to some axis of symmetry. This type of anisotropy is common in engineering composites [28], and in biological materials [21], as well as being present in certain types of stone. Let the -axis be this axis of symmetry. The elastic stiffness tensor of such an orthotropic, transversely isotropic medium contains five independent components. In its principal axes, and using shorthand notation, can be arranged as


with the stress tensor and engineering strains arranged into vectors and of the form


Note the factor of 2 applied to the shear strains to convert from the tensor strain to the engineering strain. With this shorthand notation, the stress-strain relation is


2.2 Energy densities and the dissipation potential

One useful property of the poroelasticity system is that it admits an energy density, from which we can define an energy norm that we will use extensively.

This section is mainly based on Biot’s 1956 papers [4, 5] and chapter 7 of Carcione’s book [13]. All formulations are in terms of the following variables:

  1. , the displacement vector of the solid matrix

  2. , the relative motion of the fluid scaled by the porosity, where is the displacement vector of the pore fluid and is the porosity of the medium

  3. , the variation in fluid content

2.2.1 Strain energy

In terms of the undrained elasticity tensor and the strain components of the solid matrix and , , the strain energy density of the Biot model for 3D transversely isotropic materials, with -axis being the axis of symmetry, is given by


where for a transversely isotropic material and the undrained elastic constants are related to those of the dry matrix, , via


Here and are the bulk moduli of the constituent material of the solid matrix and of the pore fluid, respectively. The total stress (solid matrix plus pore pressure) acting on a volume element of the medium are given by the derivatives of strain energy with respect to the associated strains,


or, in short notation,


Similarly, the pore pressure is


Most of the work in this paper is for 2D plane-strain conditions in the - plane. For these conditions, we will later need a function that maps from the pressure and the in-plane stress components , , at a point to the strain energy density at that point. To obtain this, we will first derive the strains as a function of stress for plane-strain conditions by setting in equations (11) and (12), then solving for the remaining strains and the variation in fluid content:


where the matrices , , and are




Substituting (13) through (15) into (4) yields the strain energy as a function of stress for plane strain conditions, which can be expressed as the quadratic form


Note the use of the drained elastic coeffcients (, etc.) rather than the undrained coefficients (). We expect this quadratic form to be positive-definite on physical grounds — if it were not, then it would be possible to deform the medium or change its fluid content without doing work.

2.2.2 Kinetic energy and the dissipation potential

For anisotropic poroelastic media, the kinetic energy density has the form


where is the induced mass matrix. Assume all the three matrices can be diagonalized in the same coordinate system so that , and . By looking at the special case of no relative motion between fluid and solid, it can be shown that the mass coefficients satisfy the two equations (given as (7.169) in Carcione [13])


Defining the velocity variables


the kinetic energy density can be expressed as


where and are the constituent solid density and pore fluid density, respectively and is the bulk density of the medium. The induced mass parameters are related to the tortuosity by


Assuming the pore fluid flow is of the Poiseuille type, the dissipation potential in an anisotropic medium in terms of the dynamic viscosity of pore fluid and the permeability tensor is


Assume that has the same principal directions as , and with eigenvalues . Then we have


2.3 Equations of motion

The equations of motion for the solid part are given in terms of the energy densities and dissipation potential as


or equivalently


Similarly, the equation of motion for the fluid part is


or equivalently


which reduces to


The equations of motion for the fluid-solid composite are obtained by adding (25) with (27), giving


2.4 Governing equations for plane-strain case

Since the material is assumed isotropic in the - plane, we consider the plane strain problem in the - plane. The governing equations are obtained by suppressing the -component (subscript 2) of , , , and those terms of with or in , , , , and . (While nonzero out-of-plane stresses do arise in a plane-strain problem, they do not produce in-plane motion, and can be ignored for purposes of studying the in-plane dynamics of the medium.) We obtain two types of governing equation:

  • Stress-strain relations, obtained by differentiating (12) and (11) for with respect to time,


    where , , , and are the solid and fluid external sources.

  • Equations of motion


    where .

2.5 Governing equations as a linear first-order system

Solving (34) and (36) for and , we obtain


where . Similarly, (35) and (37) lead to


where . Combining the stress-strain relations (30)-(33) with equations (38)-(41), we obtain the linear first-order system






It is this system (42) that forms the basis for our numerical work. Note that while the coefficient matrices , , and are defined in the material principal axes, we can extend this system to model media where the principal axes are different from the global - axes through an appropriate transformation of the state variables in that come from vector and tensor quantities, and application of the chain rule in the partial derivatives with respect to the spatial variables. In such cases, we refer to the principal directions as the and axes to distinguish them from the computational and axes.

It is worth noting that (42) is not just a generic “black box” equation, but one of a very specific type: a first-order hyperbolic system with a stiff relaxation source term. We will prove hyperbolicity in section 2.7, and the source term shows itself to be of relaxation type by having only zero and negative eigenvalues — in the absence of the spatial derivative terms, it would cause the solution to decay exponentially toward the null space , and even with the other terms present we can expect it to keep the solution close to . (Whether the relaxation term really is stiff depends on the other time scales of the particular problem being solved, but it is stiff for some of the problems considered here.) On the subject of time scales, because is extremely sparse, we can immediately read off the eigenvalues associated with dissipation in the and axes — respectively, and — and so define the characteristic time for decay in each axis as the negative inverse of these eigenvalues,


2.6 Energy norm

Let be the total mechanical energy per unit volume in a representative element, where is the kinetic energy from (20), and is the strain energy from (16). For the subsequent analysis, we will use its Hessian with respect to the state variables in , which is the symmetric matrix


We know that is a positive-definite matrix because it is the Hessian of the positive-definite quadratic form . In fact, because has no linear terms in the state variables, we can write it compactly in terms of its Hessian as


For many poroelastic materials, the components of are very badly scaled relative to each other when expressed in common units – for example, waves in the geological materials of Table 1 typically have stress components about seven orders of magnitude larger than their velocity components when expressed in SI base units. This makes using the usual vector norms on problematic, but we can fix this issue by using to define an energy norm,


We define the norm using the Hermitian conjugate-transpose (superscript ) rather than the simple vector transpose in case we later want to take the energy norm of a complex vector. This norm has the advantage of scaling the elements of in a physically relevant fashion, and producing a result that is physically meaningful and has consistent units. Incidentally, this also lets us write the energy density in the even more compact form


Note that because involves the elastic moduli, density, etc. of the medium, it is different for different materials in a heterogeneous domain. Because the energy norm is derived from the energy density at a point, however, energy norms computed in different materials can still be meaningfully compared.

2.7 Hyperbolicity

Equipped with the matrix , we can now prove that the left-hand side of (42) forms a hyperbolic system.

First, consider the system formed by setting the left-hand side of (42) equal to zero,


This system is hyperbolic if the linear combination is diagonalizable with real eigenvalues for all real and . Suppose is an eigenvector of with eigenvalue . Then


or, multiplying by ,


Since is nonsingular, any pair that satisfies the generalized eigenproblem (55) also satisfies the original eigenproblem (54).

It is not obvious how this transformation is helpful, but if we examine the component matrices and of , after substantial algebra we can discover they are symmetric:


Thus is a symmetric matrix, and (55) is a real symmetric generalized eigenproblem with a positive-definite matrix on its right-hand side. As such, it has purely real eigenvalues and a full set of linearly independent eigenvectors. is therefore diagonalizable and has pure real eigenvalues, which means that (53) is a hyperbolic system.

2.8 Entropy function

We can also show that the energy density is a strictly convex entropy function of the system (42), in a sense similar to that of Chen, Levermore, and Liu [16]. Adapting the definition of [16] to the notation used here, a function is a strictly convex entropy function for the system (42) if it satisfies the following conditions:

  1. is symmetric for all scalars and

  2. for all

  3. For , the following are equivalent:

  4. is positive-definite

Here the primes indicate gradients with respect to , so is the Hessian of with respect to . The definition of Chen, Levermore, and Liu includes an additional clause in item 3 related to an operator we call ( in their notation) that maps from to the conserved quantities of the relaxation part of the system, — namely, that conditions 3(a) and 3(b) should also be equivalent to for some appropriately-sized vector . Rather than take this as part of the definition of a strictly convex entropy function, it is more convenient here to take it as a requirement on .

Suppose . From the preceding sections we already know that conditions 1 and 4 are satisfied, so it only remains to prove conditions 2 and 3. Since , condition 2 reduces to for all . Using equations (46) and (49) for and , and recalling , we find that


By inspection, is a symmetric negative-semidefinite matrix, so for all and condition 2 is satisfied.

For condition 3, note that 3(a) implies 3(b) since if , necessarily . To see that implies , note that from (57), we have if and only if . Since only has nonzero entries in the columns corresponding to and , if and only if . Therefore condition 3 holds, and is a strictly convex entropy function for the poroelastic system (42).

3 Finite volume solution method

3.1 Wave propagation

We solve the equations of poroelasticity using a Cartesian grid finite volume approach. This section describes the basics of the finite volume method used here, as well as specifics of how we apply this method to poroelasticity. For a comprehensive discussion of this class of finite volume methods, see LeVeque’s book [36].

The class of finite volume method we use here updates cell averages at every step by solving a Riemann problem between each pair of adjacent grid cells. Thinking of one cell as the “left” cell of the problem, and the other as the “right” cell, the Riemann solution process produces the left-going and right-going fluctuations and — the changes in the cell variables caused by the left-going and right-going waves — along with a set of waves with speeds that are used to implement higher-order correction terms. With these correction terms included, the methods used here are second-order accurate. Where solutions are not smooth, wave limiters can be used on the higher-order terms to prevent spurious oscillations. While limiters can reduce the asymptotic order of accuracy of the solution, they often decrease the actual value of the error, depending on the norm being used to measure it and on the grid resolution. They can also improve the qualitative behavior of the solution by suppressing dispersive errors, leading to improved estimates of quantities such as wave arrival times, and keeping total variation from increasing.

For a homogeneous first-order hyperbolic system such as (53), the left-going and right-going fluctuations are related to the waves and wave speeds by


For a linear problem, such as linear poroelasticity, the waves are simply eigenvectors of the flux Jacobian matrix (for instance, the matrix of (44) for waves propagating in the direction in a material with its principal axes aligned with the coodinate axes) associated with the material through which the wave propagates — that is, they have the form , where is the eigenvector and is a scalar that gives the strength of the wave. Each wave speed is the corresponding eigenvalue of the flux Jacobian.

A quantity of critical importance in these solution methods is the CFL number . Informally, the CFL number is the ratio of the distance a wave travels in one timestep to the width of a grid cell; more formally, for a Cartesian grid the global CFL number is


Here and are the speeds of waves generated from the Riemann problems in the and directions, and are the grid spacings, and is the time step size. The methods used here are stable for ; since comes from a maximum over all waves, this means that our stability is limited by the fast P wave.

Because the poroelasticity equations we use here are a linear system, solution of the Riemann problem is straightforward. There is one complication, however. We wish to consider domains composed of multiple materials — in fact, our code has the capability for each grid cell to be made of a different material — so the coefficient matrices are only piecewise constant. We always choose the grid boundaries to coincide with the material boundaries, but we must still solve Riemann problems between domains with different coefficient matrices and . Suppose we are solving a Riemann problem in the -direction, with in the left cell and in the right cell. We know that in the left cell, the Riemann solution consists of waves with strength in the directions of the eigenvectors of , corresponding to the negative eigenvalues of ; similarly, in the right cell we will have waves with strength in the directions of eigenvectors or , corresponding to the positive eigenvalues of . There will also be a stationary discontinuity at the cell interface, which will lie in the null space of and . (The fact that has the same null space for any poroelastic material greatly simplfies matters here, and the corresponding eigenvectors and will not carry a subscript identifying them with the left or right material.) The total jump in across all the waves and the stationary discontinuity must add up to the difference in between the left and right states, , so we require


We can thus compute the wave strengths as . In practice, since the strength of the stationary discontinuity is never used directly, we never compute and . The same analysis holds for a Riemann problem in the -direction. This approach corresponds to an open-pore condition between the two poroelastic media, as described by Deresiewicz and Skalak [25] and validated by Gurevich and Schoenberg [32]. Other interface conditions are possible in poroelasticity, such as the closed or partially open pore conditions of Deresiewicz and Skalak, or the loose contact condition of Sharma [44]; we do not model these in this work, but they would be straightforward to incorporate into the Riemann solution process.

Because the eigenstructure of poroelasticity is somewhat complex, we do not compute the eigensystems of and analytically. Instead, we use LAPACK [1] to compute the eigenvalues and eigenvectors of and for every poroelastic material present in the model, and the matrices for each Riemann solve direction and every pair of left and right materials that could occur. For efficiency, we pre-compute these quantities for all materials used (or all possible pairs of materials in the case of ) before starting the solution proper, and look them up using a material number stored with each cell during the Riemann solves.

Boundary conditions were implemented using the usual ghost cell approach [36]. We set ghost cell values using either zero-order extrapolation, for boundaries where waves should flow outward and not return, or by setting the ghost cell values equal to the exact solution at the centers of those cells, when we verified our code against known analytic solutions.

3.2 Operator splitting

We include the dissipative part of the poroelasticity equations using operator splitting. Since the matrix is constant, we can use the exact solution operator to advance the solution by a time increment ; not only is this the most accurate solution available for this part of the system, it is also unconditionally stable and allows the time step to be chosen based solely on stability for the wave propagation part of the system.

The software framework we use offers either Godunov or Strang splitting as a run-time option. Godunov splitting is formally first-order accurate in time and uses a single full-length step of the source term operator per time step, while Strang splitting is second-order and uses two half-steps of the source term. For many practical problems Godunov splitting is a good choice because it displays similar similar error to Strang — the coefficient of the first-order error term is often small — while being less computationally intensive. However, we primarily use Strang splitting here because it displays substantially greater accuracy for the particular poroelasticity problems we solve, and because our source term is computationally cheap compared to the wave propagation part of the system. For comparison, we also show results for Godunov splitting.

3.3 Stiff regime and subcharacteristic condition

For some cases we consider, the time step is much larger than the characteristic time scales associated with the solution of . These cases fall outside the regime where asymptotic leading-order error estimates are relevant, and for them the source term is stiff, a known source of difficulty in operator splitting approaches for hyperbolic equations [18, 37]. Based on a conjecture of Pember [41], we expect to avoid spurious solutions for this stiff relaxation system if the poroelasticity equations (42) satisfy a subcharacteristic condition, where the wave speeds for the reduced equations obtained by restricting the full system to the equilibrium manifold of the dissipation term (i.e. zero fluid velocity relative to the solid matrix) interlace with the wave speeds for the full system. The appropriate generalization of Pember’s conjecture to systems of more than two equations is not obvious, but based on the principle that information should propagate more slowly (certainly no more quickly!) in the reduced system than in the full system, we expect to avoid spurious solutions if for all possible wave propagation directions the speeds and of the reduced system are strictly less than the speed of the fastest wave of the full system,


Here is the speed of a fast P wave. We ignore the negative eigenvalues for both the full and reduced systems, because they are simply the negatives of the wave speeds and will automatically satisfy a similar inequality. We also ignore the zero eigenvalues of the full system, since they correspond to eigencomponents of the solution that are left unchanged in the wave propagation part of the solution process, and are evolved according to the exact solution operator in the dissipation part.

To construct the appropriate reduced system, we follow the derivation of Chen, Levermore, and Liu [16]. First we examine the dissipation part of the system in isolation,


System (62) has six conserved quantities , which are related to the state variables by the matrix


The fact that is a vector of conserved quantities of (62) follows immediately from the fact that , so that . Given the conserved quantities , the unique equilibrium of (62) that satisfies and