HighResolution Finite Volume Modeling of Wave Propagation in Orthotropic Poroelastic Media
Abstract
Poroelasticity theory models the dynamics of porous, fluidsaturated 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 Cartesiangrid highresolution 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 timestepping approach. To the authors’ knowledge, this is the first use of highresolution finite volume methods to model poroelasticity. The solution code uses the clawpack finite volume method software, which also includes blockstructured adaptive mesh refinement in its amrclaw variant. We present convergence results for known analytic plane wave solutions, achieving secondorder 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/poro2dcartesianarchive.
oroelastic, wave propagation, finitevolume, highresolution, 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 fluidsaturated 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 semianalytical 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 timeharmonic 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 timedomain 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 lowfrequency 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 highfrequency 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 timeintegration 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 QuirogaGoode 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 velocitystress formulation of Biot linear orthotropic poroelasticity theory using Cartesiangrid highresolution finite volume methods. These methods are memoryefficient 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 BergerColellaOliger blockstructured 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 firstorder 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 firstorder 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 Stressstrain 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
(1) 
with the stress tensor and engineering strains arranged into vectors and of the form
(2) 
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 stressstrain relation is
(3) 
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:

, the displacement vector of the solid matrix

, 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

, 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
(4) 
where for a transversely isotropic material and the undrained elastic constants are related to those of the dry matrix, , via
(5)  
(6)  
(7)  
(8)  
(9) 
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,
(10) 
or, in short notation,
(11) 
Similarly, the pore pressure is
(12) 
Most of the work in this paper is for 2D planestrain conditions in the  plane. For these conditions, we will later need a function that maps from the pressure and the inplane 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 planestrain conditions by setting in equations (11) and (12), then solving for the remaining strains and the variation in fluid content:
(13) 
where the matrices , , and are
(14) 
and
(15) 
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
(16) 
Note the use of the drained elastic coeffcients (, etc.) rather than the undrained coefficients (). We expect this quadratic form to be positivedefinite 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
(17) 
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])
(18) 
Defining the velocity variables
(19) 
the kinetic energy density can be expressed as
(20)  
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
(21) 
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
(22) 
Assume that has the same principal directions as , and with eigenvalues . Then we have
(23) 
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
(24) 
or equivalently
(25) 
Similarly, the equation of motion for the fluid part is
(26) 
or equivalently
(27) 
which reduces to
(28) 
2.4 Governing equations for planestrain 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 outofplane stresses do arise in a planestrain problem, they do not produce inplane motion, and can be ignored for purposes of studying the inplane dynamics of the medium.) We obtain two types of governing equation:
2.5 Governing equations as a linear firstorder system
Solving (34) and (36) for and , we obtain
(38)  
(39) 
where . Similarly, (35) and (37) lead to
(40)  
(41) 
where . Combining the stressstrain relations (30)(33) with equations (38)(41), we obtain the linear firstorder system
(42) 
where
(43)  
(44)  
(45)  
(46) 
and
(47) 
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 firstorder 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,
(48) 
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
(49) 
We know that is a positivedefinite matrix because it is the Hessian of the positivedefinite quadratic form . In fact, because has no linear terms in the state variables, we can write it compactly in terms of its Hessian as
(50) 
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,
(51) 
We define the norm using the Hermitian conjugatetranspose (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
(52) 
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 lefthand side of (42) forms a hyperbolic system.
First, consider the system formed by setting the lefthand side of (42) equal to zero,
(53) 
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
(54) 
or, multiplying by ,
(55) 
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:
(56) 
Thus is a symmetric matrix, and (55) is a real symmetric generalized eigenproblem with a positivedefinite matrix on its righthand 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:

is symmetric for all scalars and

for all

For , the following are equivalent:


is positivedefinite
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 appropriatelysized 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
(57) 
By inspection, is a symmetric negativesemidefinite 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 leftgoing and rightgoing fluctuations and — the changes in the cell variables caused by the leftgoing and rightgoing waves — along with a set of waves with speeds that are used to implement higherorder correction terms. With these correction terms included, the methods used here are secondorder accurate. Where solutions are not smooth, wave limiters can be used on the higherorder 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 firstorder hyperbolic system such as (53), the leftgoing and rightgoing fluctuations are related to the waves and wave speeds by
(58) 
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
(59) 
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
(60) 
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 openpore 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 precompute 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 zeroorder 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 runtime option. Godunov splitting is formally firstorder accurate in time and uses a single fulllength step of the source term operator per time step, while Strang splitting is secondorder and uses two halfsteps 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 firstorder 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 leadingorder 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,
(61) 
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,
(62) 
System (62) has six conserved quantities , which are related to the state variables by the matrix
(63) 
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