Poromechanics: from Linear to Nonlinear Poroelasticity and Poroviscoelasticity

Résumé — Poromécanique : de la poroélasticité linéaire à la poroélasticité non linéaire et la poroviscoélasticité — Compte tenu des répercussions sur les calculs de productivité et de réserves, une modélisation fiable du comportement des roches est essentielle en ingénierie de réservoir. Cet article étudie plusieurs aspects du comportement poroélastique des roches dans le cadre de la mécanique des milieux poreux saturés définie par Biot. Les lois de comportement de la poroélasticité linéaire et non linéaire sont tout d’abord établies à partir d’une décomposition fondamentale de l’état de contraintes, qui permet de relier de façon claire les modèles linéaire et non linéaire. Les notions de contrainte effective et de compressibilité sont abordées. Une loi de comportement incrémentale linéaire est définie à partir de la loi de comportement non linéaire proposée en introduisant des propriétés élastiques tangentes. Ces coefficients sont par nature exprimés en fonction des déformations et de la pression de pore, mais l’on établit Abstract — Poromechanics: From Linear to Nonlinear Poroelasticity and Poroviscoelasticity — Due to the impact on productivity and oil in place estimates, reliable modeling of rock behavior is essential in reservoir engineering. This paper examines several aspects of rock poroelastic behavior within the framework of Biot’s mechanics of fluid saturated porous solids. Constitutive laws of linear and nonlinear poroelasticity are first determined from a fundamental stress decomposition, which allows to clearly connect linear and nonlinear models. Concept of effective stress and rock compressibility are considered. Linear incremental stress-strain relations are derived from the proposed nonlinear constitutive law by defining tangent elastic properties. These characteristics are naturally functions of strains and pore pressure, but explicit expressions as functions of stresses and pore pressure are established herein. Experiments performed on a reservoir sandstone illustrate these points. A constitutive law of poroviscoelasticity is finally presented and applied to experimental data obtained on clay.


INTRODUCTION
Predicted production and restitution rates rely on the accurate representation of reservoir poromechanical behavior.Most engineering reservoir models use rock compressibility as a constant mechanical key parameter.Yet, various meanings are still attached to this concept.Moreover rock properties can be stress-dependent and/or time-dependent.This paper studies rock poroelastic behavior within the framework of Biot's mechanics of fluid saturated porous solids.
Constitutive laws of linear and nonlinear poroelasticity are first determined from a fundamental stress decomposition, which allows to clearly connect linear and nonlinear models.The physical meaning of Biot's and Terzaghi's effective stresses is clarified.The definition of formation compressibility and expressions of Zimmerman's compressibilities in terms of linear poroelastic properties are recalled.In order to derive linear incremental stress-strain relations from the proposed nonlinear constitutive law, tangent elastic properties depending on strains and pore pressure are defined and expressed as explicit functions of stresses and pore pressure.Experiments performed on reservoir sandstone illustrate these points.A constitutive law of poroviscoelasticity is then presented and applied to experimental data obtained on clay.
The models proposed in this paper are derived for an isotropic homogeneous porous material, saturated with a viscous compressible fluid, under the hypothesis of small perturbations (Coussy, 1994;Charlez, 1991).Only isothermal transformations are considered.

FUNDAMENTAL STRESS DECOMPOSITION
Consider a representative elementary volume of a saturated bulk porous material submitted to stresses σ = and pore pressure p p .The stress decomposition defined in Figure1 allows to derive the constitutive law from the behavior of the dry material (i.e.without any fluid in the porous space) and the behavior of the matrix (solid and unconnected porosity) (Biot, 1973).The total loading state [T] is simply described by the superposition of two elementary loading states: [B] (for bulk) stresses σ = and zero pore pressure; [S] (for solid matrix) stresses -p p 1 = and pore pressure p p (1) where σ = amounts to the effective stress of Terzaghi.Note that the stress convention of the mechanics of continuous media is adopted.Hence -p p 1 = denotes a compressive stress.Component [B] involves no fluid pressure.Hence, from a mechanical point of view, it corresponds to the loading of the dry material by a stress field , whereas component [S] corresponds to the hydrostatic loading of the solid matrix with pressure p p .
Van der Knapp brought out that elastic porous media submitted to isotropic loading exhibit a partial linear behavior with regard to the solid matrix compressibility (Van der Knapp, 1959).Biot extended this work to the general case by introducing the concept of semilinearity (Biot, 1973), which stipulates that the solid matrix strains depend linearly on stresses and pore pressure, whereas the strains due to the effective stresses involve nonlinear modifications of local geometries, such as changes in contact areas, crack closure... Hence, both in linear and nonlinear poroelasticity, the solid matrix is supposed isotropic and linearly elastic with bulk modulus K s and shear modulus G s .Under this assumption, elementary loading state [S] leads to a deformation state: (2) Indeed, loading [S] involves no change of the (eulerian) porosity.Hence , where denotes the strain of the solid matrix.
Thus, to the stress decomposition defined in Figure 1 corresponds the following strain decomposition: (3) where denotes the strains produced by σ = -.

Mechanical Part of the Linear Constitutive Law and Biot's Effective Stress
In addition to the assumption of isotropic and linear elastic behavior of the solid matrix, the law of linear poroelasticity relies on the following assumption: the equivalent dry material is isotropic and linearly elastic with bulk modulus K o and shear modulus G. Deformation state is then given by: (4) Let with: ( Fundamental stress decomposition (Note that, according to the stress convention of the mechanics of continuous media, compressive stresses are negative while pressures are positive).
The mechanical part of the linear constitutive law is obtained by introducing Equations (1), (3) and (5), into Equation ( 4): b is the well-known Biot coefficient and is the effective stress of linear elasticity.
Note the difference between Equations ( 4) and ( 7): Biot's effective stress is the stress that produces total strains ε = , while stress that amounts to Terzaghi's effective stress only produces strains .For an incompressible solid matrix (i.e.1/K s = 0), we simply have .Boutéca and Guéguen (1999) clearly show the different roles played by Terzaghi's and Biot's effective stresses from experiments performed on Tavel limestone for which b = 0.65.

Hydraulic Part of the Linear Constitutive Law
The variation in fluid mass content m per unit of initial volume (taken positive for an incoming fluid) can be expressed as (Coussy, 1994): (8) where φ is the (eulerian) porosity, ε = tr ε = and subscript o denotes the reference state.
After linearization, we get: (9) Loading [S] leads to no porosity change.Taking Equation (2) into account, Equation (9) leads to: (10) where K fl is the fluid bulk modulus.
As loading [B] keeps pore pressure constant, m [B] only results from the change in pore volume.Taking the relation into account, we obtain: (11) The assumptions of linear behavior of the dry material and the solid matrix yield: and , where is the mean stress in the solid matrix.As pore pressure is zero, we have .Taking Equation (5) into account, Equation (11) becomes: (12) Equations ( 10) and ( 12) lead to: (13) The hydraulic part of the linear constitutive law is obtained by introducing Equation (3) into Equation ( 13): (14) where M is Biot modulus: (15)

Rock Compressibility
Formation compressibility is defined as the in situ pore volume strains that follows changes in reservoir pore pressure: Thus, formation compressibility depends on the production induced changes in effective stress state.These changes are usually described by the reservoir stress path, defined as the change in effective minimum horizontal stress divided by the change in effective vertical stress from initial reservoir conditions.
Zimmerman's compressibilities describe the volumetric response of a sandstone to applied isotropic stresses and pore pressure p p (Zimmerman, 1991): The so-defined porous rock compressibilities can be expressed as a function of linear poroelastic properties K o and b (Boutéca, 1992).
Equation (6) yields: From which we get: Expressions of C pc and C pp are derived from the change in pore volume: Taking Equations ( 4), ( 1) and ( 6) into account, we finally obtain: From which we get: Note that for an incompressible solid matrix (i.e.b = 1), we simply have: The actual stress path followed by the reservoir can be significantly different from the one predicted by the isotropic loading or uniaxial strain boundary models (Bévillon, 2000).Reservoir simulations conducted for weakly cemented sandstone show a significant increase in production rate when changing the formation compressibility from a high stress path value to a low value (Ruisten et al., 1996).Moreover, for highly compacting reservoir, the bulk compressibility depends on stresses and pore pressure (see Paragraph 3).The concept of rock compressibility is thus to be handled cautiously.

Mechanical Part of the Nonlinear Constitutive Law
From a mechanical point of view, loading [B] corresponds to the loading of the dry material.Stresses and strains are thence related through the free energy F o of the theory of elasticity: where F o is a function of the three strain invariants: . Biot adopted an expression derived from the second-order elasticity theory (Biot, 1963).The free energy F o is then the sum of the classical potential of linear elasticity W and a nonlinear potential H: where K o is the drained bulk modulus, G is the shear modulus and D, F, N, are nonlinear parameters.
Introducing Equations ( 17) and ( 18) into Equation ( 16) leads to the nonlinear stress-strain relation associated with loading [B], which is to be compared to Equation ( 4): (20) The mechanical part of the nonlinear constitutive law is then given by: ( 21) and is to be compared to Equation (6).

Hydraulic Part of the Nonlinear Constitutive Law
Due to the assumption of linear behavior of the solid matrix (see Paragraph 1), elementary variation in fluid mass content m [S] is still given by Equation ( 10): The expression of elementary variation in fluid mass content m [B] includes an additional nonlinear term.Indeed, the introduction of Equation ( 20) into Equation (11) yields: (23) Equations ( 22) and ( 23) lead to Equation ( 24), which is to be compared to Equation ( 13): (24) The hydraulic part of the nonlinear constitutive law is obtained by introducing Equation (3) into Equation ( 24): (25) and is to be compared to Equation ( 14).

Additional Hypothesis
We make herein the reasonable assumption that there is no effect of pure distortion on the variation in fluid mass content.Introducing Equation (19) into Equation (24) yields: where the invariant accounts for the distortion since it becomes zero for .Hence, the above hypothesis implies: N = 2D.
The nonlinear constitutive law is then written: where When considering the principal directions, Equation (26a) gives: (28a) Other principal stresses are simply derived by cyclic permutation.
Figure 2 Definition of a tangent drained bulk modulus.

Definition of Tangent Properties
Nonlinear elastic behavior is often handled by deriving linear incremental stress-strain relations and introducing tangent elastic properties depending on the actual state (see Eq. ( 29) and Fig. 2).

Tangent Shear Modulus for Axisymmetric Loading
Deriving an explicit expression of the tangent shear modulus is a more complex issue, which involves the deviatoric stress q and the deviatoric strain ε d : (34a) (34b) Note that we have q = qand ε d = ε d .However, a simple expression of G t can be obtained under the assumption of axisymmetric loading, which remains a rather general case.Taking for example axis 1 as axis of symmetry (i.e.ε 22 = ε 33 ) and assuming that |σ 11 | > |σ 22 |, Equations (34a) and (34b) read: Equation (28c) gives then: Differentiating Equation ( 35) finally yields: (36) which allows to derive an explicit expression of the tangent shear modulus as a function of the deviatoric strain ε d : (37)

Linear Incremental Constitutive Law for Axisymmetric Loading
A simple calculation shows that under the assumption of axisymmetric loading, the nonlinear constitutive law can be written in the following linear incremental form: where tangent properties K t o , G t and b t , are given by Equations ( 31), ( 32) and (37).(41) is an "effective stress" in that it relates the drained bulk modulus and Biot coefficient to stresses and pore pressure (Boutéca and Guéguen, 1999).

Expression of G t
In the case of an axisymmetric loading, Equation (37) defines the tangent shear modulus as a function of ε d .As ε d = 0 when q = 0, solving Equation ( 35

Application to Experimental Results
Cyclic loading was performed on a 19.8% porosity reservoir sandstone, which contains mainly quartz (70%), the remaining part including feldspar, mica and clay minerals.The experimental setup consists of a pressure cell and two volumetric pumps (Laurent et al., 1993).The sample is coated with impermeable materials, which allows application of pore pressure and isotropic confining pressure.
Nine experiments were conducted on the same sample.Figures 3a and 3b show the applied stress paths (Boutéca et al., 1998).Experiment 1 is a loading cycle.The loading procedure consists of an increase of the confining pressure at constant pore pressure, followed by a loading sequence repeated several times: confining pressure is first maintained constant while pore pressure is increased, then confining pressure is increased at constant pore pressure.The unloading procedure is performed along the same stress path.Experiments 2 to 5 are loading cycles conducted at constant pore pressure, loading and unloading following the same stress path.Experiment 6 is a loading cycle similar to cycle 1.In experiment 7, confining and pore pressures are increased following the same load path.In experiment 8, pore pressure is decreased while confining pressure is maintained constant.And, in experiment 9, confining pressure is decreased at constant pore pressure.
Stress-strain curve corresponding to experiment 1 is plotted in Figure 4. Despite a slight hysteresis, the reservoir sandstone shows a nonlinear elastic behavior.

Checking of the Semilinearity Assumption
According to the definition of semilinearity, the solid matrix strains depend linearly on stresses and pore pressure.Experiment 7 allows to check this hypothesis.The stress path followed in this experiment involves no effective stress σ -.
Experiment 1. Stress-strain curve (despite a slight hysteresis the reservoir sandstone shows a nonlinear elastic behavior).
Figure 5 Simulation of experiment 7. Stress-strain curve (checking of the semilinearity assumption: linear behavior of the matrix).

Figure 3b
Stress paths of experiments 6 to 9.
According to Equation (2), ∆ε should then be a linear function of ∆p p : (44) where ∆ε = εε o and ∆p p = p p -p o p (index o denoting the initial state).
Figure 5 shows that it is actually the case.One can then consider that Biot's semilinear theory should apply.Equation (44) gives: K s = 33 505 MPa.

Simulation of Loading Conducted at Constant Pore Pressure
Loading conducted at constant pore pressure p p = p o p allows a simple identification of the semilinear model described by Equations ( 1), ( 3), ( 27) and ( 28).The volumetric strain variations are then directly related to the mean effective stress σ -= -(p c -p o p ).So that loading and unloading stages may be compared, the strain variations are calculated with respect to the strain corresponding to the lowest confining pressure (this reference state being denoted by index o).
Writing Equation (28a) for actual and reference states and subtracting the two obtained equations yield: Taking Equations ( 31) and ( 32) into account, we finally get: where and are given by Equations ( 40) and ( 41).
Equations ( 47a) and (47b) depend on three material characteristics: K o , (D + 2F)/3 and K s .Parameters K o and (D + 2F)/3 are derived from Equation (47a) and experiment 9, whereas K s has been derived from experiment 7 (see Table 1).The data of experiment 9 are plotted in Figures 6 and 7 together with the theoretical results computed from Equations ( 47a) and (47b).These curves show that, in experiment 9, the rock is actually behaving as a semilinear porous material according to Biot's theory.Simulation of experiment 9. Fluid mass variation as a function of strain (checking of the semilinear behavior of the rock using the second constitutive equation).

Effect of Cycling
In laboratory experiments, an evolution of the measured rock properties is often observed when the sample undergoes cyclic solicitations.Cycling finally ends up in a non variable rock sample, from which accurate and reproducible measurements can be obtained.
Figure 8 shows the data of experiments 1, 6 and 9, plotted together with the semilinear model identified from experiment 9.The experimental curves appear to gradually move towards the theoretical trend derived from experiment 9. We believe that the behavior exhibited in experiment 9 is the intrinsic behavior of the rock and that this intrinsic behavior was obtained resorting to cycling (Boutéca et al., 1998).
Figure 10 Variation law of the tangent drained bulk modulus.

Tangent Bulk Modulus
Experimental data confirm that the tangent bulk modulus varies according to the pressure difference p c -p p as indicated by Equation ( 40) (Boutéca et al., 1994).This is illustrated in Figure 9, where estimations of K t o computed from experiment 6 are plotted as a function of the confining pressure for different pore pressure values and as a function of p c -p p .
Figure 10 shows the variation law of the tangent drained bulk modulus derived from Equation (40) and which represents the rock intrinsic behavior (see Paragraph 3.6.3).The estimations of K t o computed from experiment 9 have been plotted for comparison.Tangent drained bulk modulus as a function of p c -p p (experiment 6).

POROVISCOELASTICITY
Further details on poroviscoelastic modeling can be found in Coussy (1994).

Constitutive Law
Actual strain and variation in fluid mass content are defined with respect to a reference relaxed state (i.e. a state of thermodynamic equilibrium where no evolution of the state variables occur).They can be decomposed in elastic parts (superscript e), immediately recovered through unloading, and viscous parts (superscript v), not instantaneously recovered through unloading: (48 where φ v = m v /ρ fl o appears as the viscous porosity.Note that φ v is a lagrangian porosity (i.e.defined with respect to the initial volume).
Adopting a second-order expression for the free energy leads to Equation ( 49) and constitutive law (50), where K 0 , K 0 o , G 0 , b 0 and M 0 , are the material instantaneous characteristics (with respect to viscous phenomena and not hydraulic diffusion). (49 When adopting a second-order expression for the viscous dissipative potential as well, we get the following complementary evolution law: (51) where η is the shear viscosity and ζ the bulk viscosity.

Instantaneous and Relaxed Characteristics
Consider an experiment in which the following stress and fluid pressure history is imposed: where Y(t) is the Heaviside step function that reads Y(t) = 0 if t < 0 and Y(t) = 1 if t ≥ 0. Let and m 0 be the instantaneous response.When considering finite steps and ∆p, Equation ( 51) shows that the viscous strain rate can not be infinite and thus that the viscous strain is continuous.Hence the instantaneous viscous strain is zero and the instantaneous response is purely elastic.Equations ( 50) yield then: (52) Characteristics K 0 , K 0 o , G 0 , b 0 and M 0 , can therefore actually be interpreted as instantaneous material properties.
Let and m ∞ be the asymptotic response (i.e.obtained after an infinite time).This response corresponds to a new relaxed thermodynamic state, for which strain rates are zero.Equations ( 50) and (51) lead then to: and M ∞ are the material relaxed properties.Equations ( 54) show that parameters β, κ and γ , connect instantaneous and relaxed characteristics.

Application to Experimental Results
Uniaxial strain tests including drained and undrained loading cycles were performed on a shale.Undrained unloading allows a simple application of Zener model.Consider an initial equilibrium state defined by: Undrained mechanical unloading corresponds to an instantaneous decrease of the axial load, with no fluid flow through the sample upstream and downstream sides.The associated boundary conditions are then given by: The undrained condition gives the additional relation: Since we consider uniaxial strain tests, all the quantities depend only on variables z and t.
The momentum equation ( ) and boundary condition (55a) show that the axial stress is constant: The conservation of fluid mass and Darcy's law, lead to the equation of fluid diffusion: (57) Taking Equation (55d) into account, Equations ( 57) and (55b) show that the variation in pore pressure inside the sample is uniform.
Equations ( 50) and (51) yield: From which we get the following differential system: (60) with Taking Equations (55d) and (56) into account, Equation (60) leads to the differential equations that govern the evolution of pore pressure as a function of time: is the instantaneous pressure variation, is the asymptotic pressure variation.
The differential equations that govern the evolution of pore pressure as a function of time is then: (63) Experimental data show that the pore pressure equilibrium time is small when compared to the strain equilibrium time (see Fig. 11).Thus, the variation in pore pressure can be considered instantaneous with respect to the strain evolution, which amounts to disregard the expressions including the exponential term in Equations ( 62) and ( 63).We finally get: (64) where is the instantaneous axial strain, is the asymptotic axial strain.
Figure 12 Simulation of the axial strain evolution during an undrained unloading (checking of Zener proviscoelastic model).
Experimental data showing the axial strain evolution during an undrained unloading are plotted in Figure 12 together with the theoretical results computed from Equation (64).Zener poroviscoelastic model appears to represent satisfactorily shale time-dependent behavior.

CONCLUSION
Reliable modeling of rock behavior is essential in reservoir engineering due to the impact on productivity and oil in place estimates.This paper has attempted to review the main poroelastic models.
Constitutive laws for linear and nonlinear poroelasticity, and poroviscoelasticity, have been presented and illustrated by experimental applications.The concepts of effective stress and compressibility, and their various meanings, have been considered.Explicit expressions of tangent properties as functions of stresses and pore pressure have been derived for rocks presenting a nonlinear behavior.We hope that this document provides the fundamental theoretical laws needed to represent the poroelastic behavior of most reservoir and cap rocks.
Figure 3aStress paths of experiments 1 to 5.
Figure 6Simulation of experiment 9. Stress-strain relation (checking of the semilinear behavior of the rock using the first constitutive equation).

Figure 7
Figure7 Figure9 p o = p (t = 0) -p o Figure 11Illustration of the almost instantaneous nature of the pore pressure variation in undrained condition.