Variations of the in‐situ stress within nonhomogeneous layers of the subsurface rock induced by changes in pore fluid pressure

This paper presents a three‐dimensional analytical solution for deformation and stress change inside and outside a porous layer of rock subjected to prescribed axisymmetric disturbance in the pore fluid pressure. The porous layer is confined between two impermeable layers of rock. Different mechanical moduli are considered for the porous and confining strata. The solution is derived from Navier's equations of static equilibrium while considering the continuity of stress and deformation across the rock layers which are assumed to follow linearly poroelastic constitutive behavior. Results for the case of uniform pore fluid pressure change within the porous layer are demonstrated in terms of rock deformation, as well as changes in the magnitude and orientation of principal stresses, inside and outside the porous layer. Variations of these results with key problem parameters including the stiffness ratio of rock layers and the aspect ratio of disturbed pore pressure volume are presented. Findings indicate the large effect of layers heterogeneity on redistribution and reorientation of the in‐situ stress, yet, less significantly on deformation of the porous layer. Results further reveal the substantial error rates that uniaxial strain models may produce in estimating the stress path of depleted or injected reservoirs.

effective stress ′ ℎ to the corresponding change in the vertical effective stress ′ . 2,13 The effective stress is, ′ = + , where denotes the corresponding total stress component (tensile positive), is the pore fluid pressure, is Biot's coefficient, and is the Kronecker delta function. The uniaxial strain model of poroelastic rock deformation predicts the following linear relation between the horizontal stress change and the pore fluid pressure change inside the depleted reservoir.
Equation (1) advocates no changes in the vertical stress owing to the assumption of uniaxial rock deformation. Therefore, the stress path predicted by this model would be = ∕(1 − ), where is Poisson's ratio of the rock. Because of its simplicity, Equation (1) or variations thereof have been widely used in the petroleum geomechanics literature. [14][15][16] Despite the popularity and common use, field evidence 17,18 indicates that the true reservoir stress path depends on location and may be significantly different from the predictions made by a uniaxial model. More involved solutions based on time-independent, decoupled theory of poroelasticity and axisymmetric models of depleting reservoirs in a homogenous model of the subsurface are available in the literature. 8,19 The former solution is recently extended to the case of viscoelastic rock behavior in Ref. 20 The corresponding time-dependent and fully coupled solution in terms of axisymmetric Green's functions for a semi-infinite poroelastic medium is published in Ref. 21 Homogeneity of the subsurface rock layers is a common, yet, oversimplifying assumption in the existing analytical solutions for large-scale geomechanics problems. Due to their mathematical complexity, numerical approaches are often adopted to obtain the corresponding displacement and stress solutions to the pertinent heterogenous problems. 22 However, a limited number of analytical solutions to such problems have been reported on the subject. An example is the stress solution to axisymmetric reservoirs above a rigid basement layer. 23 A one-dimensional, poroelastic stress solution to the problem of reservoir fluid injection via a three-layer model of reservoir, base rock, and caprock is reported in Ref. 24 Through a different approach, the theory of ellipsoidal inhomogeneities 25 is successfully adopted in ref. 4,26,27 to estimate the change in the stress field of ellipsoidal reservoirs upon depletion. These solutions predict the same stress path as the one obtained by the uniaxial model when the aspect ratio of the ellipsoidal inhomogeneity inclusion grows sufficiently large. A major shortcoming of ellipsoidal inhomogeneity solutions is that the predicted stress change field within the ellipsoidal volume is uniform. This outcome is in contradiction with field evidence, for example, stress arching 28 effect which advocates nonuniform redistribution of the in-situ stress within depleted reservoirs.
The solution to the problem of stress path inside and outside an axisymmetric volume of pore fluid pressure within a homogenous and infinite model of the subsurface is recently reported in Ref. 29 This paper presents a generalization of the work by Su et al. 29 to the case of heterogeneous rock layers. For this purpose, an analytical solution is developed for displacement and stress inside and outside an axisymmetric volume of prescribed pore fluid pressure change within a porous rock layer. The porous layer is confined by two semi-infinite burden layers of rock with dissimilar mechanical properties. These burden layers have much smaller permeabilities than the confined porous layer.
Biot's theory of poroelasticity 30 is used to describe the constitutive behavior of the considered rock layers. The solution is based on a general representation of the displacement and stress fields within each rock layer in the Hankel transform space. The solution to the considered problem is obtained by securing the continuity of displacement and stress at the interface between the permeable layer and surrounding layers. The proposed solution is verified against the numerical solution of COMSOL R Multiphysics. 31 Parametric analysis of the solution is conducted for the stiffness ratio of the impermeable and permeable layers, = 2 ∕ 1 , and the thickness-to-diameter ratio ℎ∕(2 ) of perturbed pore pressure volume. Further, distributions of permeable layer compaction and stress change inside and outside the permeable layer are presented. Lastly, the changes in magnitude and orientation of in-situ principal stress components are thoroughly visualized and discussed. Figure 1 shows a three-layer model of the subsurface which consists of a permeable middle layer being confined by two semi-infinite rock formations from the top and bottom. All rock layers are porous and fluid-saturated. The permeable layer is assumed to be sufficiently deep to allow for neglecting the effect of traction-free ground surface on the developed stresses in its general vicinity. A quantified assessment of this assumption is later made in the model verification section. The top and bottom layers' permeability is sufficiently small to render an undrained rock behavior upon mechanical deformation.

ANALYTICAL SOLUTION
F I G U R E 1 Schematics of depleted rock layer and confining rock strata with dissimilar mechanical properties.
The permeable middle layer has different mechanical properties compared to the impermeable layers and undergoes a prescribed axisymmetric distribution of pore fluid pressure change. As such, the considered model is time-independent, therefore, does not capture the relevant coupled effects of transient poroelasticity, for example, Noordbergum effect 32 or dilative intake effect. 33 However, linear poroelastic constitutive relations 34 are herein used to derive the solution for the mechanical response of rock layers to the arbitrary axisymmetric profile of pore fluid pressure change, ( , ), within the permeable middle layer. The basic set of equations for the considered model is classified and outlined, as follows.
• The relation between the strain tensor, , and displacement vector, : • Navier's equilibrium equation of static stress tensor, : • The mixed-stiffness poroelastic constitutive relations for porous rocks: In Equations (1)-(5), is the pore pressure and is the poroelastic kinematic variable for fluid content increment. is the rock shear modulus, is Poisson's ratio, is Biot's coefficient, and is Biot's modulus. Biot's modulus can be written in terms of , , and undrained Poisson's ratio as, By substituting the Equation (4) into Equation (3) and considering the strain-displacement relation from Equation (1), Navier's equations of static equilibrium for porous solids are obtained, as follows.
where, = is known as Geertsma's uniaxial compaction coefficient. 8 Given the axisymmetric nature of the problem at hand, Hankel integral transform is used to derive the displacement and stress field of the rock layers in Figure 1. The th order Hankel transform of function and the corresponding inverse transform are recalled, as follows, 35̂( By using the condition of axial symmetry and applying Equation (9) on Equation (8), the volumetric strain of the permeable layer can be obtained in Hankel transform domain as, Here, (sub-)subscript 1 denotes a property of the permeable middle layer, and 1, ( ) denotes the th arbitrary function of integration pertaining to the solution for the same. By substituting Equation (11) into the -component of Equation (7), the radial displacement component, , of permeable layer can be obtained in Hankel transform domain as, . The volumetric strain of a cylindrical coordinate system can be expressed as, Consequently, by substituting Equations (11) and (12) into Equation (13), the vertical displacement of permeable layer in Hankel transform domain can be derived as, The displacement and strain solutions of the permeable layer can then be substituted into constitutive equation Equation (4) to obtain the stress components. The resulting stress expressions of the permeable layer are outlined, as follows.
1 + 1 1 1,1 1 − 1 1 1,2 − 1 + 1,3 1 + 1,4 The solution for the impermeable layer is obtained, similarly. Here, the pure-stiffness form of the poroelastic constitutive equations, 34 along with the following expression for Navier's equations of static equilibrium, are more convenient to use. , Similarly, the Beltrami-Michell compatibility equation can be expressed in terms of strain tensor and fluid content change, , as, where = . Equations (19)- (21) are the counterparts of Equations (4), (7), and (8) in which the pore fluid pressure, , is replaced by the pore fluid content increment, . With using similar derivations, the general displacement and stress components of the impermeable layer can be obtained as, where 2 = 1∕[2(1 − 2 2 )]. 2, ( ) denotes the th arbitrary function of integration of the solution for impermeable layers. Here subscript 2 indicates a property of the impermeable layer. As such, there are a total of eight unknown functions of integration for the permeable and impermeable layers solutions, that is, 1,1 , 1,2 , 1,3 , 1,4 , 2,1 , 2,2 , 2,3 , 2,4 . Since the considered solution domains of the impermeable layers are infinite, 2,1 and 2,3 must be zero so that the displacement and stress values remain finite at the far field, that is, where 2 → ∞. In addition, the low permeability of the impermeable rock layers justifies consideration of an undrained response in the same layers. This condition may be mathematically imposed by letting = 0 in Equations (22)- (27).
The continuity of displacements, as well as the and components of stress tensor, across the interface between the permeable and impermeable layers ( 1 = ℎ 2 ) are expressed through the following equations.
1. Radial displacement is continuous along the interface, .
The resulting displacement and stress solutions for the permeable and impermeable layers are written, as follows.
• Complete solution for the permeable layer: • Complete solution for the impermeable layer: where = 2 ∕ 1 . The constants 1 − 8 and 1 − 6 appearing in Equations (34)-(45) are given in Appendix B. Due to the complexity of the solution listed through Equations (34)-(45), numerical integration is required to obtain the resulting displacement or stress terms. In the following sections, the use of NIntegrate function of Wolfram Mathematica 36 for this purpose is shown to return satisfactory results. In particular, the integrals were successfully taken using NIntegrate function without a requirement for separate asymptotic expansions of the integrands for small or large values of the involved Bessel function arguments.

SOLUTION VERIFICATION
The analytical solution that was derived in Section 2 is herein verified against the numerical solution to the same problem. COMSOL R Multiphysics 31 is used for this purpose. The axisymmetric numerical model and discretization of solution domain into finite elements are shown in Figure 2. The burden layers are undrained while the permeable layer has a negative uniform pore pressure change, that is, pore fluid pressure reduction within a cylindrical volume of 500 ft radius and 50 ft thickness. The input parameters for each layer, along with the dimensions of the numerical solution domain, are given in Table 1. Since same values of Poisson's ratio and undrained Poisson's ratio are assigned to both the permeable layer and impermeable layer in the following sections, shear modulus ratio is same as Young's modulus ratio of the impermeable to permeable layer, that is, = 2 ∕ 1 = 2 ∕ 1 . To eliminate the far field boundary effects, the dimensions of the overall numerical solution domain are set to be 60 and 600 times larger than the radial and vertical dimensions of the disturbed pore fluid pressure volume, respectively. The numerical solution domain is discretized using an adaptive mesh of unstructured, triangular finite elements with quadratic shape functions. Partial views of the infinite subsurface and permeable layer are presented in Figure 2B,C. Consistent with the numerical model, a uniform pore pressure reduction 1 = Δ H( − ) is used in the analytical solution. Expression of the same pressure decline profile in Hankel transform space may be made by substitution in  Table 2 for radial and tangential stress components and in Table 3 for vertical and shear stress components. Positive values imply tension or counterclockwise rotation tendency. The differences of all four stress components predicted by the two solutions are less than 1% with all three stiffness ratios = 0.5, 1, 2 at all three points.
To further evaluate the effect of neglecting the free ground surface of the subsurface in the analytical solution, the numerical model is next modified to incorporate a traction-free boundary on the top surface of the numerical model shown in Figure 2A. This new numerical model is shown in Figure 3A where the permeable layer is located at depth below the traction-free surface. A comparison of the obtained stress components at point P1 between the numerical and analytical solutions is demonstrated in Figure 2C. Here, the stress component obtained from COMSOL simulation is denoted by whereas the same stress component obtained from the analytical solution is referred to as, . The horizontal axis in Figure 3B shows the relative difference between the two solutions, that is, ( − )∕ . The vertical downward axis of the same plot is the depth of porous layer measured from the traction-free boundary of the numerical TA B L E 2 Comparison of radial and tangential stresses between analytical and numerical solutions   model representing the ground surface. The shear stress component shows the least convergence rate between the two solutions with rates of relative difference that are smaller than 5% for 2 > 1.6. Convergence of other stress components between the two solutions occurs at substantially shallower depths.

DISCUSSION OF THE SOLUTION RESULTS
Parametric analysis of the developed analytical solution against the solution key parameters is presented in this section. The studied parameters include the ratio of the impermeable layer stiffness to that of the permeable middle layer, = 2 ∕ 1 , as well as the aspect ratio of the disturbed pore fluid pressure volume, ℎ∕ (2 ). These parameters are shown in Figure 1. The values in Table 1 are used as the input parameters of the solution. A negative uniform pore fluid change, that is, pore pressure reduction, is considered in the permeable layer to conduct and discuss this parametric analysis. However, the results and discussions in this section may be conveniently converted to the fluid injection case by assuming a positive pore pressure change in the permeable layer.

In-situ stress change
Variations of stress components with stiffness ratio, , are presented in Figure 4 for the point ( ∕ , 2 ∕ℎ) = (0.5, 0.1) within the impermeable layer. The stresses on vertical axes of the plots are scaled by Δ and presented in the tensilepositive convention. The gray strips in Figure 4 denote the typical range of stiffness ratios, 0.1 < = 2 ∕ 1 < 10, that may occur in practice. However, for the sake of completeness of the pertinent demonstrations, variations of stress components are plotted against a much larger range of values. An extremely small would imply that the impermeable layers are much softer than the permeable layer. As a result, the induced stress change of the impermeable layer should diminish when → 0. With increasing the stiffness ratio, , the normalized radial and tangential stress changes become more compressive until ≅ 2. This can be explained by less compaction at the boundary of disturbed pressure volume compared to the center axis of the same. In this case, the permeable layer is compressing the surrounding layers. As increases, the surrounding layers become stiffer. In such a case, the described compression effect would be negated by the depletion-induced Poisson's effect in the undrained surrounding layers rendering tensile stress changes in the impermeable layers. Consequently, the normalized radial and tangential stress changes convert to tensile mode and asymptotically approach a constant value as the value of grows. For extremely small aspect ratios of the disturbed pore volume, that is, ℎ∕(2 ) ≪ 1, this asymptotic value would be −(1 + 2 )∕2. This result can be directly produced from the elastic solution of a semi-infinite space subjected to a uniformly distributed vertical load on the surface. The derivation for this special case is presented in Appendix C. Figure 4C shows that the vertical stress change is always tensile and asymptotically converges to − Δ when grows large. Figure 4D presents the shear stress changes with stiffness ratio . Indeed, smaller aspect ratios would yield smaller shear stress magnitudes. Furthermore, the induced shear stress is generally smaller than the normal components of induced stress by one or two orders of magnitude. However, variations of shear stress with are more complex and not as intuitively explicable.
Variations of dimensionless stress change ∕( Δ ) with radial distance from the axis of symmetry of the impermeable layer at 2 ∕ℎ = 0.1 are presented in Figure 5. The pore fluid change induces compressive radial and tangential stress changes in impermeable layers for the cases of = 0.5 and 1. Due to Poisson's effect, a larger stiffness ratio yields tensile, radial, and tangential stress changes in the impermeable layer. The induced radial and tangential stress changes diminish when grows large. The vertical stress change is tensile with < . It converts to compressive changes when crossing the lateral boundary of perturbed pore pressure volume. In either region, the magnitude of vertical stress change is larger for a larger stiffness ratio. The shear stress is relatively small if the point of interest is not close to the lateral boundary of disturbed pore pressure volume where high gradients with a finite peak value at the boundary are observed. The stress profiles along the vertical coordinates of the layers, ( 1 , 2 ), at different radii are presented in Figure 6 where axes of all plots denote the vertical distance from the symmetry plane, 1 = 0. Discontinuities in the radial and tangential components of stress are observed when crossing the interface between the impermeable and permeable layers in Figure 6A,B. As expected, the vertical stress change is tensile at all radii and reaches a maximum value at a depth close to the center plane of the permeable layer.

Deformation of the porous layer
Compaction or bulging of the permeable layer refers to the difference between the vertical displacements of the top and bottom surfaces of the permeable layer. Because of the model symmetry about the center plane of the permeable layer, compaction or bulging of permeable layer, ( ) can be calculated as, (46) F I G U R E 7 Compaction at the center of the permeable layer with variations of stiffness ratio, , for different thickness-to-diameter ratios, ∕(2 ).

F I G U R E 8
Compaction profiles of the permeable layer at various stiffness ratios, = 2 ∕ 1 . Figure 7 shows compaction of the permeable layer for Δ < 0 at the center of depleted pore fluid pressure volume, (0), versus stiffness ratio = 2 ∕ 1 for four aspect ratios, ℎ 2 = 0.1, 0.05, 0.01, and 0.001. Compaction is normalized by ,hm = 1 Δ ℎ which is the compaction of a disk-shaped volume with extremely thin aspect ratio, ℎ∕ ≪ 1, of pore pressure disturbance within a homogeneous model ( = 1) of the subsurface. 8 Figure 7 shows that ,hm = 1 Δ ℎ is the maximum possible compaction for any stiffness ratio or aspect ratio. An asymptotic value of = (1 − ) 1 Δ ℎ is observed for → 0 which is the case when the permeable depleted layer has a much larger stiffness than the surrounding impermeable strata. Same conclusions apply to the case of fluid injection resulting in pore fluid pressure increase (Δ > 0) and the permeable layer bulging instead of compaction. Figure 8 presents the compaction profiles of the permeable layer, ( ), with radial distance from the vertical axis of symmetry at different ratios, = 2 ∕ 1 . For all stiffness ratios, compaction is maximum at the axis of symmetry of the permeable layer. Compaction decreases significantly when crossing the lateral boundary of pore pressure change volume and gradually diminishes when moving further away from the boundary.

In-situ principal stress change and rotation
Rotation and change in magnitudes of principal stress components because of disturbances in the pore fluid pressure change within the permeable layer are discussed in this section. For this purpose, a case study including normal faulting regime with initial gradients of overburden stress, ∕ = −0.85 psi∕f t, horizontal stresses, ∕ = −0.78 psi∕f t, and ℎ ∕ = −0.72 psi∕f t, is considered. The initial pore fluid pressure gradient is ∕ = 0.65 psi∕f t. The permeable layer depth, , is assumed to be 16,700 ft. A uniform pore pressure decline, Δ = −5427.5 psi, is considered. The in-situ stress changes both inside and outside the permeable layer are obtained from the analytical solutions. A Cartesian coordinates system ( , , ) is considered with = 0 plane residing on the plane of symmetry at the middle of the permeable layer in Figure 1 and the axis being aligned with the direction of . The stress components in Section 2 are derived in cylindrical coordinate system ( , , ). The corresponding stress components in Cartesian coordinates system can be obtained as follows, The polar angle in Equation (48) is measured from the positive axis. The stress components induced by reservoir depletion are rotated onto the Cartesian coordinate system ( , , ), and then superposed on the initial in-situ stresses to give, The magnitude and direction of post-depletion principal stresses can be obtained from the eigenvalues and eigenvectors of . Consequently, the corresponding changes in magnitudes of three principal stresses can be found as, ) and Δ 3 = ( 3 − ℎ ), respectively. Figure 9 presents the changes in the magnitudes of three principal stresses at 2 ∕ℎ = 0.05 within the impermeable layer for the cases of = 1∕3 and = 3. The radial direction of the polar charts represents the radial distance from = 0 axis while the tangential values of the polar charts correspond to the azimuth that follows the geographic convention of being measured clockwise from the positive axis. For both = 1∕3 and = 3 cases, the first principal stress (formerly overburden or vertical stress) would reduce except for a region close to the reservoir lateral boundary where it becomes more compressive. The magnitude of the first principal stress change is larger with stiffer burden layers. With soft burden layers, that is, = 1∕3, the second and third principal stresses become more compressive due to the reservoir depletion. However, because of Poisson's effect, a stiffer burden would cause larger changes of the second and third principal stresses, particularly, in the direction of = 0 • . Figure 10 shows the variations of first in-situ principal stress which used to be the overburden (vertical) stress before depleting the porous layer. The cases of = 1∕3 and = 3 are presented at the cross section of = 0 • . The color plots show the change in magnitude of first principal stress scaled by depletion pressure, Δ , while the arrows show the new orientation of the first principal stress. In general, the first principal stress becomes more tensile at < while a reverse effect is observed when > . This observation holds for both the permeable and impermeable layers and is somewhat augmented at the region close to the lateral boundary of the disturbed volume in the permeable layer. In the same region, the first principal stress rotates, therefore, would be no longer vertical.

CONCLUSIONS
This paper presented a three-dimensional analytical solution for stress and displacement inside and outside an axisymmetric volume of pore pressure change within a permeable layer of rock that is sandwiched by two impermeable layers with dissimilar mechanical properties. Findings indicate the rather nontrivial effect of mechanical inhomogeneity in between the layers on the re-equilibrated state of the in-situ stress. The magnitudes and signs of changes in the in-situ principal stress components vary depending on location inside either of the layers, contrast between the mechanical properties of the layers, as well as the geometric aspect ratio of the perturbed pore pressure volume. The orientation of two horizontal stress components may switch while substantial tilt of the first principal stress may occur in both the permeable and the impermeable layers at a region close to the lateral boundary of disturbed pore pressure. Maximum compaction of permeable layer occurs in the homogenous case while any contrast in mechanical moduli would yield smaller compaction.

A C K N O W L E D G M E N T S
The Petroleum Research Fund of the American Chemical Society under grant number 59894-DNI9 is gratefully acknowledged.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The authors confirm that the data supporting the findings of this study are available within the article. Consequently, the six unknown constants can be obtained as, ] .