# Coupled Geomechanics and Flow Modeling of Fractured Reservoirs considering Matrix Permeability Anisotropy

Authors: Dawei Wu; Yuan Di; Yu-Shu Wu
Journal: Geofluids
 (2022)

---

## Abstract

The effect of geomechanics is crucial in the modeling of fractured reservoirs since fractures can be more stress-sensitive than the rock matrix. In fractured reservoirs, the microscale fractures are often homogenized into the matrix continuum leading to matrix permeability anisotropy. The embedded discrete fracture model (EDFM) is becoming increasingly popular for numerical simulation of fractured reservoirs and has been successfully employed in the simulation of coupled geomechanics and flow systems. However, the anisotropic permeability of matrix cannot be considered in traditional EDFM. An approach of coupled geomechanics and flow modeling is proposed for fractured reservoirs considering anisotropic matrix permeability. In order to calculate the fluid exchange between fractures and rock matrix in an anisotropic formation, an integrally embedded discrete fracture model (IEDFM) is used. The geomechanics model of the proposed approach uses an equivalent continuum approach, which introduces an equivalent material to represent the overall deformation of the fractured rock under normal and shear stresses. The constitutive relations of the equivalent continuum are derived rigorously from stress-strain analysis, where the stress-dependent moduli of natural fractures are included. The coupled geomechanics and flow system is solved using the fixed-stress split iterative coupling strategy with the dynamic hydraulic parameters of matrix and fractures updated separately. Several examples are performed to demonstrate the applicability of the proposed approach for modeling the coupled geomechanics and flow system in fractured reservoirs considering anisotropic permeability. The effect of anisotropy is investigated, which indicates that the dynamic behavior of a fracture is highly orientation-related due to initial stress anisotropy and matrix permeability anisotropy. Simulations also show that the anisotropic matrix permeability affects the compaction in the reservoir domain, which reflects on the performance of production.


---

## Body

## 1. Introduction

During the exploitation of hydrocarbon reservoirs, pore pressure changes caused by production can lead to rock compaction, which in return affects the hydraulic properties such as porosity and permeability. In fractured reservoirs, the effect of geomechanics is more crucial since fractures can be more stress-sensitive than the rock matrix. In order to capture the effect of geomechanics on the production of fractured reservoirs, it is necessary to couple flow and geomechanics.The fundamentals of coupled geomechanics and flow problems date back to the concept of effective stress first proposed by Terzaghi [1]. Biot [2, 3] later extended this concept and derived the quasistatic constitutive equations that relate the strain tensor to the stress tensor and fluid pressure, which is generalized as the three-dimensional poroelastic theory. In the poroelastic theory, rock blocks are composed of compressible soil skeleton and fluid that fills the pore space. However, in a fractured porous media, the single porosity poroelastic theory is no longer applicable. According to the dual continuum concept introduced by Barenblatt et al. [4], the matrix and fractures can be envisioned as two separate continua that interact with each other. Based on the dual continuum concept, a dual porosity poroelastic theory (DPP) is proposed [5–7], where the effective stress law is modified to capture different extents of deformation of the rock matrix and fractures. Then, Wittke [8] rigorously derived the constitutive model for a fractured rock using the concept of equivalent continuum. Later, series of attempts to simulate coupled geomechanics and flow problems in fractured reservoirs using the dual continuum concept and the dual porosity poroelastic theory are performed [9–15].Though dual continuum models provide efficient ways of modeling fractured reservoirs, these models cannot handle reservoirs with complex fracture networks composed of large-scale fractures since the real fracture geometries cannot be represented [16, 17]. The embedded discrete fracture model (EDFM) is proposed by Li and Lee [18] to improve this limitation and later successfully implemented for various fracture geometries [19–21]. In the EDFM, fractures are treated explicitly to capture the complex geometry of fracture networks. The rock matrix is discretized using Cartesian grids to keep the computational efficiency, and the fracture elements are represented by individual control volumes that are connected to the parent matrix blocks. The EDFM has become an increasingly popular method for modeling fractured reservoirs and has been successfully adopted in coupled geomechanics and flow modeling. Moinfar et al. [22] proposed a coupled flow and geomechanics model for fractured reservoirs using the EDFM. Sangnimnuan et al. [23] presented a coupled flow and geomechanics model using the EDFM to predict stress evolution in unconventional reservoirs during production. Liu et al. [24] investigated the effects of hydraulic fracture and natural fracture properties on production behavior as well as pressure and stress evolution of shale gas reservoirs using the EDFM and FEM coupled simulations. Bai et al. [25] developed a coupled compositional flow and geomechanics simulator using the EDFM and investigated the effect of confined phase behavior in shale oil reservoirs. Ren et al. [26] proposed a coupled XFEM-EDFM hybrid model for geomechanics and flow in fractured reservoirs, which uses the extended finite element method (XFEM) to model the geomechanics and the EDFM to model the flow equations. This hybrid model is later applied in the modeling of the fractured vuggy porous media [27], hydraulic fracture propagation process [28], and enhanced geothermal system [29].Eventhough extensive research works have been done on the coupled geomechanics and flow modeling of fractured reservoirs using the EDFM, the effect of matrix permeability anisotropy is often overlooked, which might significantly affect the performance of the reservoir in coupled geomechanics and flow simulations. Naturally fractured reservoirs contain multiple scaled fractures, including a limited number of macroscale fractures and a large number of microscale fractures. Though the macroscale fractures can be treated explicitly and efficiently using the EDFM, individually treating all the large amount of microscale fractures as control volumes using the EDFM often causes large computational burden. One practical method to treat these microscale fractures is to homogenize them into the matrix continuum, which results in matrix permeability anisotropy. Traditional EDFM cannot precisely model the anisotropic flow between fracture and matrix due to the oversimplified linear pressure distribution assumption applied in the EDFM [30]. Shao et al. [31, 32] developed an integrally embedded discrete fracture model (IEDFM), where fractures are envisioned as series of point sinks to derive the semi-analytical pressure distribution around the fractures. The IEDFM allows for the calculation of matrix-fracture connection transmissibilities considering the anisotropy of rock matrix, which is used to calculate the fluid exchange between the rock matrix and fracture.A coupled geomechanics and flow model for the simulation of fractured reservoirs considering matrix anisotropy is proposed in this study. The flow simulation uses the IEDFM to account for the fluid exchange between fracture control volumes and the rock matrix with anisotropic permeability. In the geomechanics modeling, an equivalent continuum approach is used to derive the stiffness matrix required in the finite element method, where an equivalent material is introduced to represent the overall deformation of fractured rock. An empirical joint model is employed to capture the stress-dependent fracture moduli. Several examples are performed to demonstrate the applicability of the proposed approach for coupled geomechanics and flow modeling of fractured reservoirs with anisotropic permeability, and the effect of anisotropy in coupled geomechanics and flow systems is investigated.


## 2. Governing Equations

The governing equations of a coupled geomechanics and flow system consist of two parts: the mass conservation equation and the force balance equation.


### 2.1. Mass Conservation Equations

The mass conservation equations in fractured media are given by(1)∂∂tϕSβρβ+∇⋅ρβvβ=qβ,where β represents the fluid phases, ϕ is the porosity, Sβ is the saturation of phase β, ρβ is the density of phase β, vβ is the velocity of phase β, and qβ is the source term of phase β. This mass balance equation applies to both matrix and fracture control volumes.The velocity of phaseβ is defined by Darcy’s law:(2)vβ=−krβμβK⋅∇Pβ−ρβg∇Z,where K is the absolute permeability tensor, krβ is the relative permeability of phase β, μβ is the viscosity of phase β, Pβ is the pressure of phase β, g is the gravitational acceleration, and Z is the depth.Discretizing equation (1) using the finite volume method in space and first-order scheme in time gives(3)ϕSβρβin+1−ϕSβρβinViΔt=∑j∈ηiFβ,ijn+1+Qβ,in+1,where subscript i denotes the values of element i, superscript n+1 represents the current time level, superscript n represents the previous time level, Fβ,ij is the flow term for phase β between element i and element j, ηi is a set that contains all the neighbor elements connected to element i, and Qβ,i is the source term for phase β.The fluid flow through matrix-matrix, fracture-fracture, and matrix-fracture connections is universally described by the flow termFβ,ij, which can be expressed as(4)Fβ,ij=TijρβkrβμβijPβ,j−Pβ,i−ρβ,ijgZj−Zi,where Tij is the transmissibility of the connection, which is the harmonic average of the half-transmissibility of element i and element j:(5)Tij=TiTjTi+Tj,(6)Ti=Aijkidi,Tj=Aijkjdj,where Ti and Tj are the half-transmissibilities, ki and kj are the absolute permeabilities along each element center to the interface, di and dj are the distances from the centers to the interface, and Aij is the interface area.


### 2.2. Force Balance Equations

The vectors and tensors are expressed in components. The force balance equation can be obtained from linear momentum conservation law, which gives(7)ρu¨is=∂σij∂xj+ρfi,where ρ=1−ϕρs+ϕ∑ρβSβ is the density of the fluid-solid mixture, ρs is the density of rock skeleton, u¨is is the acceleration vector of the rock skeleton, σij is the total stress tensor, and fi is the body force vector.In a fractured porous media, the rock matrix and fracture are envisioned as two separate continua for the flow system; thus, the traditional effective stress law is not appropriate. The effective stress law for fractured media can be derived by extending the single porosity effective stress law, which gives(8)σ′ij=σij−αMPMδij−αFPFδij,where σij′ is the effective stress tensor, PM=∑PMβSβ is the average pressure of matrix blocks, PF=∑PFβSβ is the average pressure of fractures, δij is the Kronecker delta, αM=1−Kb/Kr is the Biot coefficient of matrix, αF=Kb/Kr−Kb/Ks is the Biot coefficient of fracture, Kb is the bulk modulus of porous media, Ks is the modulus of solid grain, and Kr is the drained bulk modulus of porous media without fractures.Since the loadings caused by fluid depletion do not change drastically, deformation in the reservoir can be assumed to be quasistatic. In the case of static analysis, the dynamic term in equation (7) can be subtracted. Then, substituting equation (8) into equation (7) gives(9)∂σij′∂xj+αM∂PM∂xi+αF∂PF∂xi+ρfi=0.In the case of an elastic tensor, the constitutive relations are given by(10)σij′=Dijklεkl,where Dijkl is the fourth-order modulus tensor, and εkl is the strain tensor.Spatially discretizing equation (9) in incremental form using the finite element method gives the discretized force balance equation:(11)KmΔuN+αMKvΔPM+αFKvΔPF=ΔF,where Km=∫VBTDBdV is the nodal stiffness matrix, B is the strain-nodal displacement matrix, D is the elastic stiffness matrix, which is the matrix form of the fourth-order modulus tensor Dijkl, ΔuN is the nodal displacement increment vector, Kv is the volumetric stiffness vector, ΔPM and ΔPF are the pressure increments in matrix and fracture, and ΔF is the external loading vector.


## 2.1. Mass Conservation Equations

The mass conservation equations in fractured media are given by(1)∂∂tϕSβρβ+∇⋅ρβvβ=qβ,where β represents the fluid phases, ϕ is the porosity, Sβ is the saturation of phase β, ρβ is the density of phase β, vβ is the velocity of phase β, and qβ is the source term of phase β. This mass balance equation applies to both matrix and fracture control volumes.The velocity of phaseβ is defined by Darcy’s law:(2)vβ=−krβμβK⋅∇Pβ−ρβg∇Z,where K is the absolute permeability tensor, krβ is the relative permeability of phase β, μβ is the viscosity of phase β, Pβ is the pressure of phase β, g is the gravitational acceleration, and Z is the depth.Discretizing equation (1) using the finite volume method in space and first-order scheme in time gives(3)ϕSβρβin+1−ϕSβρβinViΔt=∑j∈ηiFβ,ijn+1+Qβ,in+1,where subscript i denotes the values of element i, superscript n+1 represents the current time level, superscript n represents the previous time level, Fβ,ij is the flow term for phase β between element i and element j, ηi is a set that contains all the neighbor elements connected to element i, and Qβ,i is the source term for phase β.The fluid flow through matrix-matrix, fracture-fracture, and matrix-fracture connections is universally described by the flow termFβ,ij, which can be expressed as(4)Fβ,ij=TijρβkrβμβijPβ,j−Pβ,i−ρβ,ijgZj−Zi,where Tij is the transmissibility of the connection, which is the harmonic average of the half-transmissibility of element i and element j:(5)Tij=TiTjTi+Tj,(6)Ti=Aijkidi,Tj=Aijkjdj,where Ti and Tj are the half-transmissibilities, ki and kj are the absolute permeabilities along each element center to the interface, di and dj are the distances from the centers to the interface, and Aij is the interface area.


## 2.2. Force Balance Equations

The vectors and tensors are expressed in components. The force balance equation can be obtained from linear momentum conservation law, which gives(7)ρu¨is=∂σij∂xj+ρfi,where ρ=1−ϕρs+ϕ∑ρβSβ is the density of the fluid-solid mixture, ρs is the density of rock skeleton, u¨is is the acceleration vector of the rock skeleton, σij is the total stress tensor, and fi is the body force vector.In a fractured porous media, the rock matrix and fracture are envisioned as two separate continua for the flow system; thus, the traditional effective stress law is not appropriate. The effective stress law for fractured media can be derived by extending the single porosity effective stress law, which gives(8)σ′ij=σij−αMPMδij−αFPFδij,where σij′ is the effective stress tensor, PM=∑PMβSβ is the average pressure of matrix blocks, PF=∑PFβSβ is the average pressure of fractures, δij is the Kronecker delta, αM=1−Kb/Kr is the Biot coefficient of matrix, αF=Kb/Kr−Kb/Ks is the Biot coefficient of fracture, Kb is the bulk modulus of porous media, Ks is the modulus of solid grain, and Kr is the drained bulk modulus of porous media without fractures.Since the loadings caused by fluid depletion do not change drastically, deformation in the reservoir can be assumed to be quasistatic. In the case of static analysis, the dynamic term in equation (7) can be subtracted. Then, substituting equation (8) into equation (7) gives(9)∂σij′∂xj+αM∂PM∂xi+αF∂PF∂xi+ρfi=0.In the case of an elastic tensor, the constitutive relations are given by(10)σij′=Dijklεkl,where Dijkl is the fourth-order modulus tensor, and εkl is the strain tensor.Spatially discretizing equation (9) in incremental form using the finite element method gives the discretized force balance equation:(11)KmΔuN+αMKvΔPM+αFKvΔPF=ΔF,where Km=∫VBTDBdV is the nodal stiffness matrix, B is the strain-nodal displacement matrix, D is the elastic stiffness matrix, which is the matrix form of the fourth-order modulus tensor Dijkl, ΔuN is the nodal displacement increment vector, Kv is the volumetric stiffness vector, ΔPM and ΔPF are the pressure increments in matrix and fracture, and ΔF is the external loading vector.


## 3. Modeling of Flow in Fractured Reservoirs considering Matrix Permeability Anisotropy

The IEDFM is adopted to model the fluid flow in fractured reservoirs considering matrix permeability anisotropy. In the IEDFM, rock matrix is discretized by Cartesian grids, and additional elements are added to represent the fracture control volumes. There are three kinds of connections between control volumes in the IEDFM: matrix-matrix (M-M) connections, fracture-fracture (F-F) connections, and matrix-fracture (M-F) connections.For M-M connections, the multipoint flux approximation (MPFA) method [33] is used to establish the connections between matrix elements in an anisotropic reservoir. Since the permeability difference between two connected matrix elements is relatively small, the pressure distribution inside matrix elements can be assumed to be linear. Thus, the parameters in equation (6) have clear physical meanings: di and dj are the distances from the matrix element centers to the interface, and ki and kj are the matrix permeabilities that can be calculated from the absolute permeability tensor K using MPFA.For F-F connections, the linear pressure distribution is also appropriate since the permeability difference is small. The parameters are as follows:di and dj are the distances from the fracture element centers to the common interface, and ki and kj are the fracture permeabilities kF.For M-F connections, due to the anisotropy of matrix permeability, simplified linear assumptions cannot be used for the pressure distribution around a fracture. In the IEDFM, the M-F transmissibilities are calculated semi-analytically to address the anisotropic matrix permeability. A fracture is considered as a series of point sinksSi, as shown in Figure 1. In a 2D anisotropic medium, if a point sink locates at the origin 0,0, the pressure at position x,y is given by(12)Px,y=qμ2πhkelnre+c,(13)ke=kxxkyy−kxy2,(14)re=kyykex2+kxxkey2+2kxykexy,where q is the point sink flow rate, h is the strata height, ke is the equivalent permeability of the anisotropic matrix, re is the equivalent distance, and c is a constant number.Figure 1
Illustration of IEDFM.In a matrix element, the pressure of pointX can be derived as the superposition of all the pressure drops caused by point sinks:(15)PXx,y=∑i=1NSqSiμ2πhkelnre,iX+CP,where NS is the number of point sinks, qSi is the flow rate of sink Si, re,iX is the equivalent distance between point X and sink Si, and CP is a constant related to the fracture pressure.The fractures are assumed to be equipotential; thus, several pointsFj on the fracture surface can be selected as the reference points, and these points have the same pressure PF. Applying equation (15) on every reference point forms a linear equation system:(16)PF=∑i=1NSqSiμ2πhkelnre,ij+CP,j=1,⋯,NS,where re,ij is the equivalent distance between reference point Fj and sink Si.The linear equation system can be rewritten as(17)∑i=1NSξi2πhlnre,ij=1,j=1,⋯,NS,where ξi is defined by(18)ξi=μqSikePF−CP.By solving the linear equation system, the exact expression of pressure at pointX can be obtained:(19)PX=∑i=1NSqSiμ2πhkelnre,iX−1PM−PFλ−1+PF,(20)PM=1VM∬MPXdVX,(21)λ=1VM∬∑i=1NSξilnre,iX2πhdVX,where PM is the matrix pressure which is defined as the average of all point pressures, and VM is the matrix volume.The transmissibility between the fracture elementF and the neighboring matrix element M can be derived using the matrix pressure, fracture pressure, and the flow term definition given in equation (4), which gives(22)TM−F=keλ−1∑ξi.The transmissibility given in equation (22) is calculated semi-analytically, which improves the accuracy and allows for anisotropic matrix permeability. Validation and detailed discussion of the IEDFM in anisotropic formations can be found in the work of Shao et al. [32].


## 4. Modeling of Fractured Rock Deformation

### 4.1. Equivalent Continuum Model

To derive the nodal stiffness matrixKm in equation (11), the elastic stiffness matrix D is required, which describes the deformation of a fractured rock.A fractured rock consists of two separate continua, rock matrix and fracture, which behave differently under stress. The concept of equivalent continuum is introduced to establish the stiffness matrix for the fractured rock. In the equivalent continuum model, an equivalent material is introduced to represent the overall deformation of a fractured rock. The deformation behavior of the introduced equivalent material is the combination of rock matrix deformation and the displacement of fracture planes.The constitutive relations of the equivalent material can be rigorously derived from stress-strain analysis. In case of a fractured rock crossed by a single fracture, a coordinate transformation is first performed to simplify the derivation. After the coordinate transformation, theZ-axis of the coordinate system is set to align with the normal direction of the surface of the fracture. All kinds of normal and shear loadings are applied on the rock block to obtain the stress and strain tensors, as shown in Figure 2. Young’s modulus and shear modulus of the fracture are assumed to be much smaller than the rock matrix, and the volume fraction of fracture in a fractured rock is much smaller than the rock matrix. The rock matrix is assumed to be elastic isotropic. In the local coordinate system of fracture, the compliance matrix can be written as(23)C=1ζEM−νMζEM−νMEM−νMζEM1ζEM−νMEM−νMEM−νMEMζEM+ηECF1ζGMζGM+1GF1GM+1GF,where ζ is the volume fraction of the rock matrix in the fractured rock, η is the volume fraction of the fracture in the fractured rock, EM is Young’s modulus of rock matrix, GM is the shear modulus of rock matrix, νM is Poisson’s ratio of rock matrix, EcF is the compression modulus of the fracture, and GF is the shear modulus of the fracture. The detailed derivation of equation (23) from stress-strain analysis is given in the Appendix.Figure 2
Illustration of stress-strain analysis.The compression modulusEcF describes the modulus of the laterally confined fracture medium, which is given by(24)EcF=EF1−νF1−νF−2νF2,where EF is Young’s modulus of the fracture and νF is Poisson’s ratio of the fracture.Since the volume fraction of the rock matrix is much larger than the fracture, it can be assumed thatζ≈1. The compliance matrix of the equivalent material C can be rewritten as(25)C=CM+ηCF,where CM and CF are the compliance matrixes for the rock matrix and the fracture given by(26)CM=1EM−νMEM−νMEM−νMEM1EM−νMEM−νMEM−νMEM1EM1GM1GM1GM,(27)CF=001EcF01GF1GF.The compliance matrix given in equation (25) is derived in the local coordinate system where Z-axis aligns with the normal direction of the fracture. To obtain the compliance matrix in the global coordinate system Cglobal, a coordinate transformation needs to be performed:(28)Cglobal=TTCT,where T is the coordinate transformation matrix from the local coordinate system to the global coordinate system, which is determined by the orientation and dip angle of the fracture.Substituting equation (25) into equation (28) gives(29)Cglobal=CM+ηTTCFT,where CM remains unchanged under the coordinate transformation due to the assumption of elastic isotropic rock matrix.When a fractured rock is crossed by more than one fractures, the compliance matrix of the equivalent material is derived by sequentially accumulating the contribution of each individual fracture. The overall compliance matrix of a fractured rock containing multiple fractures is given by(30)Coverall=CM+∑i=1NfηiTiTCFiTi,where Nf is the number of fractures, ηi is the volume fraction of fracture i, CFi is the compliance matrix of fracture i, and Ti is the coordinate transformation matrix from the local coordinate system of fracture i to the global coordinate system.The elastic stiffness matrixD in equation (11) is the inversion of the overall compliance matrix Coverall:(31)D=Coverall−1.


### 4.2. Constitutive Model of Fractures

The deformation of fractures under stress is significantly different from intact rock matrix, which is more complex and strongly nonlinear. Thus, a fracture constitutive model is needed to calculate the dynamic fracture modulus used in equation (27). For natural fractures, an empirical model can be employed, which is based on laboratory experimental data on interlocked rock samples with natural unfilled joints [34].For a fracture under normal stress, Young’s modulus is given by(32)EF=kniw1−σn/Δumaxkni+σn2,where kni is the initial fracture normal stiffness, w is the fracture aperture, σn is the effective normal stress, and Δumax is the maximum normal closure of fracture.The initial fracture normal stiffnesskni can be obtained from the following equation:(33)kni=−7.15+1.75JRC+0.02JCSw0,where JRC is the joint roughness coefficient, JCS is the joint compressive strength, and w0 is the zero-stress fracture aperture. In this work, the term “zero-stress” refers to the state when no stress is applied, and the term “initial” refers to the state under the initial stress field before the simulation begins.The maximum normal closure of fractureΔumax can be estimated from the following empirical relation:(34)Δumax=−0.1032−0.0074JRC+0.02JCSw0−0.251.The fracture aperturew is stress-dependent, and the fracture aperture change under normal stress is given by(35)Δw=σnkni+σn/Δumax.Fracture behavior under shear stress is more sophisticated and less well understood compared to deformation under normal stress. For simplicity, a nonlinear model that only applies to the prepeak fracture shear deformation is used to describe fracture behavior. The shear modulus of fracture is obtained from the empirical relation:(36)GF=Kjσnnj1−τRfτpeak2,where Kj is the stiffness number, nj is the stiffness exponent, τ is the shear stress, τpeak is the peak shear stress, and Rf is the failure ratio.The stiffness numberKj can be estimated using the empirical relation:(37)Kj=−17.19+3.86JRC.The peak shear stressτpeak is given by(38)τpeak=σntanJRClog10JCSσn+φr,where φr is the fracture residue angle.


### 4.3. Validation of the Fractured Rock Geomechanics Model

The proposed geomechanics model of the fractured rock is compared with the experimental data from the work of Bandis et al. [34]. The experimental data are acquired from applying gradually increasing loading on the top of single-jointed rectangular rock samples. The curves are drawn from experiments on a limestone sample with a fully interlocked joint in the first loading cycle. Young’s modulus for limestone matrix is 49.0 GPa. The joint roughness coefficient (JRC) is 7.6, and the joint compressive strength (JCS) is 157 MPa. The zero-stress aperture of the fracture is 0.2 mm.Figure3 shows the comparison of the normal deformation of solid rock, fracture, and fractured rock between experimental data and calculated results. The good agreement between experimental and calculated data validates the equivalent continuum model and the fracture constitutive model. The deformation of solid rock is approximately linear, indicating the rock matrix can be assumed to be elastic. The deformation of the fracture shows nonlinear behavior and almost stops increasing at high stress level, since the maximum normal closure is about to be reached. The deformation of fractured rock is a combination of fracture and solid rock deformation, which verifies the assumption of the equivalent continuum model.Figure 3
Validation of the fractured rock geomechanics model.


## 4.1. Equivalent Continuum Model

To derive the nodal stiffness matrixKm in equation (11), the elastic stiffness matrix D is required, which describes the deformation of a fractured rock.A fractured rock consists of two separate continua, rock matrix and fracture, which behave differently under stress. The concept of equivalent continuum is introduced to establish the stiffness matrix for the fractured rock. In the equivalent continuum model, an equivalent material is introduced to represent the overall deformation of a fractured rock. The deformation behavior of the introduced equivalent material is the combination of rock matrix deformation and the displacement of fracture planes.The constitutive relations of the equivalent material can be rigorously derived from stress-strain analysis. In case of a fractured rock crossed by a single fracture, a coordinate transformation is first performed to simplify the derivation. After the coordinate transformation, theZ-axis of the coordinate system is set to align with the normal direction of the surface of the fracture. All kinds of normal and shear loadings are applied on the rock block to obtain the stress and strain tensors, as shown in Figure 2. Young’s modulus and shear modulus of the fracture are assumed to be much smaller than the rock matrix, and the volume fraction of fracture in a fractured rock is much smaller than the rock matrix. The rock matrix is assumed to be elastic isotropic. In the local coordinate system of fracture, the compliance matrix can be written as(23)C=1ζEM−νMζEM−νMEM−νMζEM1ζEM−νMEM−νMEM−νMEMζEM+ηECF1ζGMζGM+1GF1GM+1GF,where ζ is the volume fraction of the rock matrix in the fractured rock, η is the volume fraction of the fracture in the fractured rock, EM is Young’s modulus of rock matrix, GM is the shear modulus of rock matrix, νM is Poisson’s ratio of rock matrix, EcF is the compression modulus of the fracture, and GF is the shear modulus of the fracture. The detailed derivation of equation (23) from stress-strain analysis is given in the Appendix.Figure 2
Illustration of stress-strain analysis.The compression modulusEcF describes the modulus of the laterally confined fracture medium, which is given by(24)EcF=EF1−νF1−νF−2νF2,where EF is Young’s modulus of the fracture and νF is Poisson’s ratio of the fracture.Since the volume fraction of the rock matrix is much larger than the fracture, it can be assumed thatζ≈1. The compliance matrix of the equivalent material C can be rewritten as(25)C=CM+ηCF,where CM and CF are the compliance matrixes for the rock matrix and the fracture given by(26)CM=1EM−νMEM−νMEM−νMEM1EM−νMEM−νMEM−νMEM1EM1GM1GM1GM,(27)CF=001EcF01GF1GF.The compliance matrix given in equation (25) is derived in the local coordinate system where Z-axis aligns with the normal direction of the fracture. To obtain the compliance matrix in the global coordinate system Cglobal, a coordinate transformation needs to be performed:(28)Cglobal=TTCT,where T is the coordinate transformation matrix from the local coordinate system to the global coordinate system, which is determined by the orientation and dip angle of the fracture.Substituting equation (25) into equation (28) gives(29)Cglobal=CM+ηTTCFT,where CM remains unchanged under the coordinate transformation due to the assumption of elastic isotropic rock matrix.When a fractured rock is crossed by more than one fractures, the compliance matrix of the equivalent material is derived by sequentially accumulating the contribution of each individual fracture. The overall compliance matrix of a fractured rock containing multiple fractures is given by(30)Coverall=CM+∑i=1NfηiTiTCFiTi,where Nf is the number of fractures, ηi is the volume fraction of fracture i, CFi is the compliance matrix of fracture i, and Ti is the coordinate transformation matrix from the local coordinate system of fracture i to the global coordinate system.The elastic stiffness matrixD in equation (11) is the inversion of the overall compliance matrix Coverall:(31)D=Coverall−1.


## 4.2. Constitutive Model of Fractures

The deformation of fractures under stress is significantly different from intact rock matrix, which is more complex and strongly nonlinear. Thus, a fracture constitutive model is needed to calculate the dynamic fracture modulus used in equation (27). For natural fractures, an empirical model can be employed, which is based on laboratory experimental data on interlocked rock samples with natural unfilled joints [34].For a fracture under normal stress, Young’s modulus is given by(32)EF=kniw1−σn/Δumaxkni+σn2,where kni is the initial fracture normal stiffness, w is the fracture aperture, σn is the effective normal stress, and Δumax is the maximum normal closure of fracture.The initial fracture normal stiffnesskni can be obtained from the following equation:(33)kni=−7.15+1.75JRC+0.02JCSw0,where JRC is the joint roughness coefficient, JCS is the joint compressive strength, and w0 is the zero-stress fracture aperture. In this work, the term “zero-stress” refers to the state when no stress is applied, and the term “initial” refers to the state under the initial stress field before the simulation begins.The maximum normal closure of fractureΔumax can be estimated from the following empirical relation:(34)Δumax=−0.1032−0.0074JRC+0.02JCSw0−0.251.The fracture aperturew is stress-dependent, and the fracture aperture change under normal stress is given by(35)Δw=σnkni+σn/Δumax.Fracture behavior under shear stress is more sophisticated and less well understood compared to deformation under normal stress. For simplicity, a nonlinear model that only applies to the prepeak fracture shear deformation is used to describe fracture behavior. The shear modulus of fracture is obtained from the empirical relation:(36)GF=Kjσnnj1−τRfτpeak2,where Kj is the stiffness number, nj is the stiffness exponent, τ is the shear stress, τpeak is the peak shear stress, and Rf is the failure ratio.The stiffness numberKj can be estimated using the empirical relation:(37)Kj=−17.19+3.86JRC.The peak shear stressτpeak is given by(38)τpeak=σntanJRClog10JCSσn+φr,where φr is the fracture residue angle.


## 4.3. Validation of the Fractured Rock Geomechanics Model

The proposed geomechanics model of the fractured rock is compared with the experimental data from the work of Bandis et al. [34]. The experimental data are acquired from applying gradually increasing loading on the top of single-jointed rectangular rock samples. The curves are drawn from experiments on a limestone sample with a fully interlocked joint in the first loading cycle. Young’s modulus for limestone matrix is 49.0 GPa. The joint roughness coefficient (JRC) is 7.6, and the joint compressive strength (JCS) is 157 MPa. The zero-stress aperture of the fracture is 0.2 mm.Figure3 shows the comparison of the normal deformation of solid rock, fracture, and fractured rock between experimental data and calculated results. The good agreement between experimental and calculated data validates the equivalent continuum model and the fracture constitutive model. The deformation of solid rock is approximately linear, indicating the rock matrix can be assumed to be elastic. The deformation of the fracture shows nonlinear behavior and almost stops increasing at high stress level, since the maximum normal closure is about to be reached. The deformation of fractured rock is a combination of fracture and solid rock deformation, which verifies the assumption of the equivalent continuum model.Figure 3
Validation of the fractured rock geomechanics model.


## 5. Coupling Strategy and Parameter Update

The coupled geomechanics and flow problem is discretized using a hybrid finite element-finite volume method (FE-FVM). The FVM is applied for the flow model with the primary variables located at the grid center, which is given in equations (3) and (4). The FEM is applied for the geomechanical model with the displacement vectors located at the grid nodes, which is given in equation (11). The coupled geomechanics and flow problem is solved using the iterative coupling strategy.


### 5.1. Fixed-Stress Split Iterative Coupling Strategy

In the iterative coupling strategy, the coupled problem is split into flow and mechanics subproblems; then, flow and geomechanics modules can be invoked sequentially to solve the coupled problem. At each iteration step, either flow or mechanics problem is solved first, and then, the other problem is solved using the solution information of the previous problem. This sequential process is repeated until the solution meets the convergence criteria, and the simulation moves on to the next time step.There are several schemes to develop an iterative coupling simulator based on the operator splitting strategies for the coupled flow and mechanics problem. The fixed-stress split scheme is employed in this work, which is recommended by Kim et al. [35]. In the fixed-stress split scheme, the flow problem is solved first while freezing the total mean stress. This scheme is unconditionally stable, and the convergence is proved by many previous works [36, 37]. This scheme is successfully applied in various attempts of solving coupled flow and mechanics problems [38, 39]. Using the fixed-stress split iterative coupling strategy, a flowchart showing the framework of solving the coupled problem is shown in Figure 4.Figure 4
Flowchart of fixed-stress split scheme.


### 5.2. Parameter Update after an Iteration Step

In every iteration step, the flow simulation is first performed, and then, the geomechanics module is called after the flow simulation converges to calculate changes in stress and strain fields using the information from the flow problem solution. According to the fixed-stress split scheme, changes in matrix parameters, such as porosity and permeability, are updated by freezing the total mean stress before the next iterative step.The variation of true matrix porosity in a deformable porous media can be derived as [40](39)dϕM=1−ϕMKb−1KsdPM−dσm,where σm is the mean total stress.By freezing the mean total stress, the reservoir porosity of rock matrix can be updated by [41](40)ϕM∗,k+1=1−εvk+1ϕMk+1−εvk+11−ϕMkKb−1KsPMk+1−PMk,where the superscript k denotes the previous iteration step, the superscript k+1 denotes the current iteration step, εv is the volume strain, and ϕM∗=1−εvϕM is the reservoir porosity, which is the ratio of the pore volume in the deformed configuration to the total volume in the undeformed configuration.The matrix permeability tensor is updated using an empirical relation [42]:(41)Kk+1=K0expψϕM∗,k+1ϕM∗,0−1,where the superscript 0 denotes the zero-stress condition and ψ is a parameter for permeability update.


### 5.3. Parameter Update after a Time Step

After the solution meets the convergence criteria, the hydraulic parameters such as porosity and permeability for matrix and fractures are updated for the next time step. Since the fracture moduli are stress-dependent, the fracture stiffness is also updated to build the equivalent constitutive matrix in the next time step.According to equation (39), the updation of matrix reservoir porosity is given by(42)ϕM∗,n+1=ϕM∗,n+1−ϕMnKb−1Ks1−ενk+1+ϕMn1Kb−1KsPMn+1−PMn−1−ϕMnKb−1Ks1−ενk+1+ϕMnKbσmn+1−σmn,where the superscript n denotes the previous iteration step, and the superscript n+1 denotes the current iteration step.The matrix permeability tensor is updated by(43)Kn+1=K0expψϕM∗,n+1ϕM∗,0−1.Normal effective stress on a fracture is calculated from the stress field and fracture orientation and then used to calculate fracture aperture change using equation (35). Porosity and permeability of fractures are updated using the updated fracture apertures:(44)wn+1=w0−Δw,(45)ϕFn+1=wn+1w0ϕF0,(46)kFn+1=kF0wn+1w02.Using the updated matrix permeability tensor, the M-F transmissibilities are updated by(47)TM−Fn+1=ken+1ke0TM−F0.The fracture moduli are updated using equations (32) and (36) before the next time loop.


### 5.4. Validation of Coupled Geomechanics and Flow Model

The accuracy of the discretized coupled geomechanics and flow model is examined by the 1D Terzaghi’s consolidation problem [1].Terzaghi’s problem consists of a soil column fully saturated with water, as shown in Figure5. A uniformly distributed load is applied on the top surface of the soil column. Drainage is only allowed on the upper boundary of the column. The lateral boundaries are fixed on horizontal direction, and the bottom boundary is fixed on both horizontal and vertical directions. In this case, the computational domain is 5 m in length and 10 m in height. The parameters used in Terzaghi’s problem are given in Table 1. The grids are shown in Figure 6.Figure 5
Diagram of Terzaghi’s problem.Table 1
Parameters of Terzaghi’s problem.
ParameterValueYoung’s modulus,E20.0 GPaPoisson ratio,ν0.2Porosity,ϕ0.3Permeability,k100.0 mDRock density,ρs1700.0 kg/m3Fluid density,ρf1000.0 kg/m3Fluid viscosity,μ1.0 cPDistributed load,T20.0 MPaFigure 6
Grids used in Terzaghi’s problem.The calculated results of pore pressure and vertical displacement distributions along the central column of grids at several times are compared with the analytical solutions [43], as shown in Figure 7. A great consistency between the analytical solutions and numerical solutions can be found, which validates the proposed coupled geomechanics and flow model.Figure 7
Comparisons of (a) vertical displacement and (b) pore pressure distribution between analytical solution and numerical solution in Terzaghi’s problem.
(a)(b)


## 5.1. Fixed-Stress Split Iterative Coupling Strategy

In the iterative coupling strategy, the coupled problem is split into flow and mechanics subproblems; then, flow and geomechanics modules can be invoked sequentially to solve the coupled problem. At each iteration step, either flow or mechanics problem is solved first, and then, the other problem is solved using the solution information of the previous problem. This sequential process is repeated until the solution meets the convergence criteria, and the simulation moves on to the next time step.There are several schemes to develop an iterative coupling simulator based on the operator splitting strategies for the coupled flow and mechanics problem. The fixed-stress split scheme is employed in this work, which is recommended by Kim et al. [35]. In the fixed-stress split scheme, the flow problem is solved first while freezing the total mean stress. This scheme is unconditionally stable, and the convergence is proved by many previous works [36, 37]. This scheme is successfully applied in various attempts of solving coupled flow and mechanics problems [38, 39]. Using the fixed-stress split iterative coupling strategy, a flowchart showing the framework of solving the coupled problem is shown in Figure 4.Figure 4
Flowchart of fixed-stress split scheme.


## 5.2. Parameter Update after an Iteration Step

In every iteration step, the flow simulation is first performed, and then, the geomechanics module is called after the flow simulation converges to calculate changes in stress and strain fields using the information from the flow problem solution. According to the fixed-stress split scheme, changes in matrix parameters, such as porosity and permeability, are updated by freezing the total mean stress before the next iterative step.The variation of true matrix porosity in a deformable porous media can be derived as [40](39)dϕM=1−ϕMKb−1KsdPM−dσm,where σm is the mean total stress.By freezing the mean total stress, the reservoir porosity of rock matrix can be updated by [41](40)ϕM∗,k+1=1−εvk+1ϕMk+1−εvk+11−ϕMkKb−1KsPMk+1−PMk,where the superscript k denotes the previous iteration step, the superscript k+1 denotes the current iteration step, εv is the volume strain, and ϕM∗=1−εvϕM is the reservoir porosity, which is the ratio of the pore volume in the deformed configuration to the total volume in the undeformed configuration.The matrix permeability tensor is updated using an empirical relation [42]:(41)Kk+1=K0expψϕM∗,k+1ϕM∗,0−1,where the superscript 0 denotes the zero-stress condition and ψ is a parameter for permeability update.


## 5.3. Parameter Update after a Time Step

After the solution meets the convergence criteria, the hydraulic parameters such as porosity and permeability for matrix and fractures are updated for the next time step. Since the fracture moduli are stress-dependent, the fracture stiffness is also updated to build the equivalent constitutive matrix in the next time step.According to equation (39), the updation of matrix reservoir porosity is given by(42)ϕM∗,n+1=ϕM∗,n+1−ϕMnKb−1Ks1−ενk+1+ϕMn1Kb−1KsPMn+1−PMn−1−ϕMnKb−1Ks1−ενk+1+ϕMnKbσmn+1−σmn,where the superscript n denotes the previous iteration step, and the superscript n+1 denotes the current iteration step.The matrix permeability tensor is updated by(43)Kn+1=K0expψϕM∗,n+1ϕM∗,0−1.Normal effective stress on a fracture is calculated from the stress field and fracture orientation and then used to calculate fracture aperture change using equation (35). Porosity and permeability of fractures are updated using the updated fracture apertures:(44)wn+1=w0−Δw,(45)ϕFn+1=wn+1w0ϕF0,(46)kFn+1=kF0wn+1w02.Using the updated matrix permeability tensor, the M-F transmissibilities are updated by(47)TM−Fn+1=ken+1ke0TM−F0.The fracture moduli are updated using equations (32) and (36) before the next time loop.


## 5.4. Validation of Coupled Geomechanics and Flow Model

The accuracy of the discretized coupled geomechanics and flow model is examined by the 1D Terzaghi’s consolidation problem [1].Terzaghi’s problem consists of a soil column fully saturated with water, as shown in Figure5. A uniformly distributed load is applied on the top surface of the soil column. Drainage is only allowed on the upper boundary of the column. The lateral boundaries are fixed on horizontal direction, and the bottom boundary is fixed on both horizontal and vertical directions. In this case, the computational domain is 5 m in length and 10 m in height. The parameters used in Terzaghi’s problem are given in Table 1. The grids are shown in Figure 6.Figure 5
Diagram of Terzaghi’s problem.Table 1
Parameters of Terzaghi’s problem.
ParameterValueYoung’s modulus,E20.0 GPaPoisson ratio,ν0.2Porosity,ϕ0.3Permeability,k100.0 mDRock density,ρs1700.0 kg/m3Fluid density,ρf1000.0 kg/m3Fluid viscosity,μ1.0 cPDistributed load,T20.0 MPaFigure 6
Grids used in Terzaghi’s problem.The calculated results of pore pressure and vertical displacement distributions along the central column of grids at several times are compared with the analytical solutions [43], as shown in Figure 7. A great consistency between the analytical solutions and numerical solutions can be found, which validates the proposed coupled geomechanics and flow model.Figure 7
Comparisons of (a) vertical displacement and (b) pore pressure distribution between analytical solution and numerical solution in Terzaghi’s problem.
(a)(b)


## 6. Example Simulations

A coupled geomechanics and flow simulator is developed to demonstrate the feasibility of the proposed model. The main programs of the flow module and mechanics module are developed on FORTRAN platform. Several examples are performed to investigate the effect of reservoir anisotropy in the modeling of the coupled geomechanics and flow system in fractured reservoirs.


### 6.1. Example 1: Effect of Anisotropy on the Closure of Fractures

During the depletion, the dynamic behaviors of fractures are largely orientation-related. The initial stress field in the reservoir can be anisotropic since it is affected by the geological conditions of the reservoir. Besides, the matrix permeability of the reservoir can be anisotropic due to the microscale fractures that are homogenized into the matrix continuum. To investigate these effects, a simple fractured reservoir is designed with three natural fractures, as is shown in Figure8.Figure 8
Configuration of the reservoir.The reservoir is 420 m in length, 420 m in width, and 40 m in depth. The reservoir domain is discretized by Cartesian grids, as shown in Figure9. The three fractures have the same length and are equally distant from the production well. Fracture 1 is parallel to the X-axis, fracture 2 is parallel to Y-axis, and fracture 3 is 45° inclined. The parameters of fractures are given in Table 1. The zero-stress aperture of fractures is 0.2 mm.Figure 9
Mesh in the reservoir.The porosity in the reservoir is 0.30. The reservoir is saturated with oil with the density of 800.0 kg/m3 and viscosity of 10.0 cP. The geomechanical configuration is shown in Figure 10. The displacement boundary conditions are as follows: the lateral boundaries are fixed on horizontal direction and the bottom boundary is fixed on all directions. The mechanical parameters are given in Table 2. The initial pressure in the reservoir is 10.0 MPa. A production well located at the center of the reservoir drains at fixed bottom hole pressure of 1.0 MPa. Constant pressure boundary condition is considered in this case; thus, the oil pressure at the boundary of the reservoir remains 10.0 MPa throughout the simulation process. The total simulation time is 3000 days.Figure 10
Geomechanical configuration.Table 2
Mechanical parameters.
ParameterValueJoint roughness coefficient, JRC6.0Joint compressive strength, JCS70.0 MPaFailure ratio,Rf0.67Stiffness exponent,nj0.75Fracture residue angle,φr30.0Matrix residue porosity,ϕr0.27Matrix permeability update parameter,ψ20.0Biot coefficient of matrix,αM0.875Biot coefficient of fracture,αF0.025The simulation is first run using isotropic initial stress field and matrix permeability as control. The initial stress is assumed to be 4.0 MPa, and the zero-stress permeability is 0.1 D. Under the initial stress, the initial matrix permeability is 0.088 D. Figure11 shows the fracture aperture evolution histories of the three fractures. The curves are in perfect agreement, since the distances from the fracture centers to the production well are the same.Figure 11
Fracture aperture histories for isotropic case.


#### 6.1.1. Effect of Initial Stress Anisotropy

The effect of initial stress anisotropy is investigated. Figure12 shows the two anisotropic initial stress fields considered in this part, and the matrix permeability is kept isotropic.Figure 12
Initial stress fields.Figure13 shows the fracture aperture histories of three fractures in the two anisotropic cases. Due to the existence of the initial stress, the initial aperture of the three fractures are smaller than the zero-stress fracture aperture. For example, when σx/σy=3, the largest closure occurs on fracture 1, which is perpendicular to the maximum initial stress direction. Fracture 1 remains only 52.6% of the zero-stress aperture before the simulation. Fracture 2, which is perpendicular to minimum initial stress direction, undergoes the smallest closure and remains 64.3% of the zero-stress aperture. Fracture 3 remains 56.2% of the zero-stress aperture. At the end of the simulation, fracture 1 remains 51.7% of the zero-stress aperture, fracture 2 remains 60.0% of the zero-stress aperture, and fracture 3 remains 54.4% of the zero-stress aperture. Because of the symmetry in geometry, the fractures behave oppositely when σx/σy=1/3.Figure 13
Effect of initial stress anisotropy on fracture aperture: (a)σx/σy=3; (b)σx/σy=1/3.
(a)(b)Therefore, fractures perpendicular to minimum initial stress direction have larger permeabilities before the simulation, yet are more prone to closure during the depletion, such as fracture 2 whenσx/σy=3 and fracture 1 when σx/σy=1/3.


#### 6.1.2. Effect of Matrix Permeability Anisotropy

Similarly, the effect of matrix permeability anisotropy is investigated. The initial stress field is assumed to be isotropic with the value of 4.0 MPa. The two anisotropic permeability tensors under initial stress are shown in Figure14.Figure 14
Initial permeability tensors.Figure15 shows the fracture aperture histories of three fractures using two anisotropic permeability tensors. Since the isotropic initial stress tensor is used, the initial apertures of the three fractures are the same. It can be found that the closure of fractures perpendicular to the maximum permeability direction is the largest, and the closure of fractures perpendicular to the minimum permeability direction is the smallest.Figure 15
Effect of matrix permeability anisotropy on fracture aperture: (a)kx/ky=4; (b)kx/ky=1/4.
(a)(b)


#### 6.1.3. The Combined Effect of Initial Stress Anisotropy and Matrix Permeability Anisotropy

Figure16 shows the comparison of the aperture evolution of fracture 1 with all the combinations of initial stress anisotropy and matrix permeability anisotropy discussed in this section. For ease of comparison, the ratio of fracture aperture to the initial aperture, w/wi, is used to represent the closure of fracture 1.Figure 16
Fracture aperture history for fracture 1 with all the combinations of initial stress anisotropy and matrix permeability anisotropy.As shown in Figure16, it can be concluded that a fracture undergoes the largest closure when it is perpendicular to the minimum initial stress direction, and the minimum initial stress direction aligns with the maximum permeability direction. Conversely, when a fracture is perpendicular to the maximum initial stress direction and the maximum initial stress direction aligns with the minimum permeability direction, the fracture closure is the smallest.


### 6.2. Example 2: A Synthetic Naturally Fractured Reservoir

Another 45 natural fractures with random length and orientations are added in the reservoir of Example 1, as shown in Figure17. The fractures are assumed to be vertical and may only penetrate specific layers of the reservoir. The lower-middle and upper-middle layers are assumed to be penetrated by more fractures than the bottom and top layers.Figure 17
Configuration of the fractured reservoir.With the additional fractures, the production well is now perforated on the fracture elements in the center of each layer, and the equivalent permeability of the reservoir is reduced to 0.1 mD. For the comparison of matrix permeability anisotropy, the initial stress field is fixed withx-direction initial stress at 6.0 MPa and y-direction initial stress at 2.0 MPa.Figures18–20 show the aperture change histories of fracture 1, fracture 2, and fracture 3, which indicate that the conclusions in Example 1 still hold for more complex fracture geometry. The closure of fractures perpendicular to the maximum permeability direction is the largest, and the closure of fractures perpendicular to the minimum permeability direction is the smallest. The aperture change of fracture 3 is only modestly influenced by the anisotropic permeability tensor since it aligns with the angular bisector of the principle permeability tensor directions.Figure 18
Fracture aperture history for fracture 1.Figure 19
Fracture aperture history for fracture 2.Figure 20
Fracture aperture history for fracture 3.The oil pressure profiles using three permeability tensors are shown in Figure21. It can be found that permeability anisotropy affects the pressure distribution in the reservoir. In the direction of the maximum permeability tensor principle axis, the pressure drop is the fastest, which induces the largest aperture change in the fractures that are perpendicular to this direction.Figure 21
Oil pressure profiles when (a)kx/ky=1; (b)kx/ky=4; (c)kx/ky=1/4.
(a)(b)(c)Since the oil pressures are used to calculate the stress, strain, and displacement in the reservoir, the permeability anisotropy directly influences the results of geomechanics simulation. The distributions of mean effective stress in the reservoir using the three permeability tensors are compared, as shown in Figure22. Stress is larger in the vicinity of fractures, which can be attributed to both the larger pressure drop around the fractures and the stiffness reduction of a fractured rock compared to intact rock.Figure 22
Mean effective stress distributions when (a)kx/ky=1, (b)kx/ky=4, and (c)kx/ky=1/4.
(a)(b)(c)The influence of matrix permeability anisotropy on fracture deformation and the geomechanics behavior of the reservoir is eventually reflected on the performance of the production well. Figure23 shows the cumulative oil produced using three permeability tensors. For comparison, the results of uncoupled simulations are also included in the figure. At 1000 days, the cumulative production is decreased by 24.4%, 25.0%, and 25.3% for kx/ky=1, kx/ky=4, and kx/ky=1/4 in coupled simulations.Figure 23
Cumulative oil production histories.


### 6.3. Example 3: Water Flooding in a Five-Spot Well System

A five-spot well system water flooding case is considered. Figure24 shows the configuration of a synthetic 2D reservoir. The primary natural fractures (marked in red) are modified from reported seismic survey data [44], and the secondary natural fractures (marked in blue) have random lengths and orientations. The parameters of primary and secondary fractures are given in Table 3. The zero-stress aperture of primary and secondary fractures are 0.3 mm and 0.15 mm. All the fractures are assumed to penetrate the entire depth of the reservoir.Figure 24
(a) Configuration of the reservoir. (b) Grids in the reservoir.
(a)(b)Table 3
Fracture parameters.
ParameterPrimary fracturesSecondary fracturesJoint roughness coefficient, JRC15.07.5Joint compressive strength, JCS80.0 MPa60.0 MPaFailure ratio,Rf0.70.75Stiffness exponent,nj0.80.65Fracture residue angle,φr30.010.0Zero-stress aperture,w00.3 mm0.15 mmThe reservoir is 500 m in length, 500 m in width, and 40 m in depth. The porosity of the reservoir is 0.20, and the initial oil saturation in the reservoir is 0.8. The initial pressure in the reservoir is 15.0 MPa. In this example, an isotropic initial stress field is applied. The initial stress is 4.0 MPa in thex-direction and y-direction and 19.0 MPa in the z-direction.Four cases with different permeability tensors are investigated. The zero-stress permeabilities tensors for four cases are given in Table4. The four tensors have the same equivalent permeability of 10.0 mD. Under the initial stress, the components of permeability tensors are reduced, and the equivalent permeability is reduced to 8.23 mD. Figure 25 shows the permeability tensors under initial stress.Table 4
Permeability tensors (zero-stress).
kxx (mD)kyy (mD)kxy (mD)ke (mD)Case 15.020.00.010.0Case 220.05.00.010.0Case 312.512.57.510.0Case 412.512.5-7.510.0Figure 25
Initial permeability tensors for four cases.Water is injected from an injection well in the center of the reservoir at the fixed pressure of 15.0 MPa. Four wells located at four corners of the reservoir produce the fixed bottom hole pressure of 5.0 MPa. The water density is 1000.0 kg/m3, and the oil density is 800.0 kg/m3. The water viscosity is 0.9 cP, and the oil viscosity is 3.0 cP. The relative permeability curves are assumed to be linear. A closed boundary is considered in this example. Therefore, all the boundaries of the reservoir are impermeable. The geomechanical configuration in this example is the same with that of Example 1. The simulation is run for 10000 days.Figure26 shows the oil saturation profiles for four anisotropic cases after 10000 days of simulation, which clearly demonstrate the effect of matrix permeability anisotropic on the water displacement process. A two-dimensional permeability tensor has two principle axes, and the direction of the maximum permeability tensor principle axis becomes a favored direction for water flooding in anisotropic reservoirs.Figure 26
Oil saturation profiles for (a) Case 1, (b) Case 2, (c) Case 3, and (d) Case 4 after 10000 days of production.
(a)(b)(c)(d)The saturation distribution shows that the existence of fractures causes strong heterogeneity across the reservoir. Due to the aperture difference between primary and secondary fractures, the primary fractures have larger impacts on the oil saturation distribution. The primary fractures become major flow channels in the reservoir, and the secondary fractures increase the heterogeneity of the reservoir. Fractures whose orientation aligns with the direction of the maximum permeability tensor principle axis create fast water displacement channels.Figure27 shows the comparison of the oil and water producing rate of four cases in four production wells. The differences in producing rates result both from the anisotropic matrix permeability tensors and the heterogeneity induced by fractures. Production wells in the direction of the maximum permeability tensor principle axis to the injection well have the highest producing rate, such as well 1 and well 3 in Case 3 and well 2 and well 4 in Case 4. Production wells in the direction of the maximum permeability tensor principle axis to the injection well have the highest producing rate, and those in the direction of the minimum permeability tensor principle axis to the injection well have the lowest producing rate.Figure 27
Histories of producing rate from (a) well 1, (b) well 2, (c) well 3, and (d) well 4.
(a)(b)(c)(d)Figure28 shows the oil pressure distribution after 10000 days. Pressure drop in the direction of maximum permeability tensor principle axis is the slowest, and pressure drop in the direction of minimum permeability tensor principle axis is the fastest. The symmetry of the pressure fields is broken by the fractures. The geometry of fractures also influences the overall pressure depletion in the reservoir since the preferential flow path created by fractures induces larger water inflow to compensate for the pressure drop. The average oil pressures in the reservoir are 10.19 MPa for Case 1, 10.13 MPa for Case 2, 10.28 MPa for Case 3, and 10.38 MPa for Case 4.Figure 28
Oil pressure profiles for (a) Case 1, (b) Case 2, (c) Case 3, and (d) Case 4 after 10000 days.
(a)(b)(c)(d)Due to the differences in pressure profiles, the deformation distributions in the reservoir in four cases are also different, which is illustrated in the volume strain profiles shown in Figure29. The average volume strains are 6.41E – 4 for Case 1, 6.49E – 4 for Case 2, 6.30E – 4 for Case 3, and 6.17E − 4 for Case 4.Figure 29
Volume strain distributions for (a) Case 1, (b) Case 2, (c) Case 3, and (d) Case 4 after 10000 days.
(a)(b)(c)(d)


## 6.1. Example 1: Effect of Anisotropy on the Closure of Fractures

During the depletion, the dynamic behaviors of fractures are largely orientation-related. The initial stress field in the reservoir can be anisotropic since it is affected by the geological conditions of the reservoir. Besides, the matrix permeability of the reservoir can be anisotropic due to the microscale fractures that are homogenized into the matrix continuum. To investigate these effects, a simple fractured reservoir is designed with three natural fractures, as is shown in Figure8.Figure 8
Configuration of the reservoir.The reservoir is 420 m in length, 420 m in width, and 40 m in depth. The reservoir domain is discretized by Cartesian grids, as shown in Figure9. The three fractures have the same length and are equally distant from the production well. Fracture 1 is parallel to the X-axis, fracture 2 is parallel to Y-axis, and fracture 3 is 45° inclined. The parameters of fractures are given in Table 1. The zero-stress aperture of fractures is 0.2 mm.Figure 9
Mesh in the reservoir.The porosity in the reservoir is 0.30. The reservoir is saturated with oil with the density of 800.0 kg/m3 and viscosity of 10.0 cP. The geomechanical configuration is shown in Figure 10. The displacement boundary conditions are as follows: the lateral boundaries are fixed on horizontal direction and the bottom boundary is fixed on all directions. The mechanical parameters are given in Table 2. The initial pressure in the reservoir is 10.0 MPa. A production well located at the center of the reservoir drains at fixed bottom hole pressure of 1.0 MPa. Constant pressure boundary condition is considered in this case; thus, the oil pressure at the boundary of the reservoir remains 10.0 MPa throughout the simulation process. The total simulation time is 3000 days.Figure 10
Geomechanical configuration.Table 2
Mechanical parameters.
ParameterValueJoint roughness coefficient, JRC6.0Joint compressive strength, JCS70.0 MPaFailure ratio,Rf0.67Stiffness exponent,nj0.75Fracture residue angle,φr30.0Matrix residue porosity,ϕr0.27Matrix permeability update parameter,ψ20.0Biot coefficient of matrix,αM0.875Biot coefficient of fracture,αF0.025The simulation is first run using isotropic initial stress field and matrix permeability as control. The initial stress is assumed to be 4.0 MPa, and the zero-stress permeability is 0.1 D. Under the initial stress, the initial matrix permeability is 0.088 D. Figure11 shows the fracture aperture evolution histories of the three fractures. The curves are in perfect agreement, since the distances from the fracture centers to the production well are the same.Figure 11
Fracture aperture histories for isotropic case.


### 6.1.1. Effect of Initial Stress Anisotropy

The effect of initial stress anisotropy is investigated. Figure12 shows the two anisotropic initial stress fields considered in this part, and the matrix permeability is kept isotropic.Figure 12
Initial stress fields.Figure13 shows the fracture aperture histories of three fractures in the two anisotropic cases. Due to the existence of the initial stress, the initial aperture of the three fractures are smaller than the zero-stress fracture aperture. For example, when σx/σy=3, the largest closure occurs on fracture 1, which is perpendicular to the maximum initial stress direction. Fracture 1 remains only 52.6% of the zero-stress aperture before the simulation. Fracture 2, which is perpendicular to minimum initial stress direction, undergoes the smallest closure and remains 64.3% of the zero-stress aperture. Fracture 3 remains 56.2% of the zero-stress aperture. At the end of the simulation, fracture 1 remains 51.7% of the zero-stress aperture, fracture 2 remains 60.0% of the zero-stress aperture, and fracture 3 remains 54.4% of the zero-stress aperture. Because of the symmetry in geometry, the fractures behave oppositely when σx/σy=1/3.Figure 13
Effect of initial stress anisotropy on fracture aperture: (a)σx/σy=3; (b)σx/σy=1/3.
(a)(b)Therefore, fractures perpendicular to minimum initial stress direction have larger permeabilities before the simulation, yet are more prone to closure during the depletion, such as fracture 2 whenσx/σy=3 and fracture 1 when σx/σy=1/3.


### 6.1.2. Effect of Matrix Permeability Anisotropy

Similarly, the effect of matrix permeability anisotropy is investigated. The initial stress field is assumed to be isotropic with the value of 4.0 MPa. The two anisotropic permeability tensors under initial stress are shown in Figure14.Figure 14
Initial permeability tensors.Figure15 shows the fracture aperture histories of three fractures using two anisotropic permeability tensors. Since the isotropic initial stress tensor is used, the initial apertures of the three fractures are the same. It can be found that the closure of fractures perpendicular to the maximum permeability direction is the largest, and the closure of fractures perpendicular to the minimum permeability direction is the smallest.Figure 15
Effect of matrix permeability anisotropy on fracture aperture: (a)kx/ky=4; (b)kx/ky=1/4.
(a)(b)


### 6.1.3. The Combined Effect of Initial Stress Anisotropy and Matrix Permeability Anisotropy

Figure16 shows the comparison of the aperture evolution of fracture 1 with all the combinations of initial stress anisotropy and matrix permeability anisotropy discussed in this section. For ease of comparison, the ratio of fracture aperture to the initial aperture, w/wi, is used to represent the closure of fracture 1.Figure 16
Fracture aperture history for fracture 1 with all the combinations of initial stress anisotropy and matrix permeability anisotropy.As shown in Figure16, it can be concluded that a fracture undergoes the largest closure when it is perpendicular to the minimum initial stress direction, and the minimum initial stress direction aligns with the maximum permeability direction. Conversely, when a fracture is perpendicular to the maximum initial stress direction and the maximum initial stress direction aligns with the minimum permeability direction, the fracture closure is the smallest.


## 6.1.1. Effect of Initial Stress Anisotropy

The effect of initial stress anisotropy is investigated. Figure12 shows the two anisotropic initial stress fields considered in this part, and the matrix permeability is kept isotropic.Figure 12
Initial stress fields.Figure13 shows the fracture aperture histories of three fractures in the two anisotropic cases. Due to the existence of the initial stress, the initial aperture of the three fractures are smaller than the zero-stress fracture aperture. For example, when σx/σy=3, the largest closure occurs on fracture 1, which is perpendicular to the maximum initial stress direction. Fracture 1 remains only 52.6% of the zero-stress aperture before the simulation. Fracture 2, which is perpendicular to minimum initial stress direction, undergoes the smallest closure and remains 64.3% of the zero-stress aperture. Fracture 3 remains 56.2% of the zero-stress aperture. At the end of the simulation, fracture 1 remains 51.7% of the zero-stress aperture, fracture 2 remains 60.0% of the zero-stress aperture, and fracture 3 remains 54.4% of the zero-stress aperture. Because of the symmetry in geometry, the fractures behave oppositely when σx/σy=1/3.Figure 13
Effect of initial stress anisotropy on fracture aperture: (a)σx/σy=3; (b)σx/σy=1/3.
(a)(b)Therefore, fractures perpendicular to minimum initial stress direction have larger permeabilities before the simulation, yet are more prone to closure during the depletion, such as fracture 2 whenσx/σy=3 and fracture 1 when σx/σy=1/3.


## 6.1.2. Effect of Matrix Permeability Anisotropy

Similarly, the effect of matrix permeability anisotropy is investigated. The initial stress field is assumed to be isotropic with the value of 4.0 MPa. The two anisotropic permeability tensors under initial stress are shown in Figure14.Figure 14
Initial permeability tensors.Figure15 shows the fracture aperture histories of three fractures using two anisotropic permeability tensors. Since the isotropic initial stress tensor is used, the initial apertures of the three fractures are the same. It can be found that the closure of fractures perpendicular to the maximum permeability direction is the largest, and the closure of fractures perpendicular to the minimum permeability direction is the smallest.Figure 15
Effect of matrix permeability anisotropy on fracture aperture: (a)kx/ky=4; (b)kx/ky=1/4.
(a)(b)


## 6.1.3. The Combined Effect of Initial Stress Anisotropy and Matrix Permeability Anisotropy

Figure16 shows the comparison of the aperture evolution of fracture 1 with all the combinations of initial stress anisotropy and matrix permeability anisotropy discussed in this section. For ease of comparison, the ratio of fracture aperture to the initial aperture, w/wi, is used to represent the closure of fracture 1.Figure 16
Fracture aperture history for fracture 1 with all the combinations of initial stress anisotropy and matrix permeability anisotropy.As shown in Figure16, it can be concluded that a fracture undergoes the largest closure when it is perpendicular to the minimum initial stress direction, and the minimum initial stress direction aligns with the maximum permeability direction. Conversely, when a fracture is perpendicular to the maximum initial stress direction and the maximum initial stress direction aligns with the minimum permeability direction, the fracture closure is the smallest.


## 6.2. Example 2: A Synthetic Naturally Fractured Reservoir

Another 45 natural fractures with random length and orientations are added in the reservoir of Example 1, as shown in Figure17. The fractures are assumed to be vertical and may only penetrate specific layers of the reservoir. The lower-middle and upper-middle layers are assumed to be penetrated by more fractures than the bottom and top layers.Figure 17
Configuration of the fractured reservoir.With the additional fractures, the production well is now perforated on the fracture elements in the center of each layer, and the equivalent permeability of the reservoir is reduced to 0.1 mD. For the comparison of matrix permeability anisotropy, the initial stress field is fixed withx-direction initial stress at 6.0 MPa and y-direction initial stress at 2.0 MPa.Figures18–20 show the aperture change histories of fracture 1, fracture 2, and fracture 3, which indicate that the conclusions in Example 1 still hold for more complex fracture geometry. The closure of fractures perpendicular to the maximum permeability direction is the largest, and the closure of fractures perpendicular to the minimum permeability direction is the smallest. The aperture change of fracture 3 is only modestly influenced by the anisotropic permeability tensor since it aligns with the angular bisector of the principle permeability tensor directions.Figure 18
Fracture aperture history for fracture 1.Figure 19
Fracture aperture history for fracture 2.Figure 20
Fracture aperture history for fracture 3.The oil pressure profiles using three permeability tensors are shown in Figure21. It can be found that permeability anisotropy affects the pressure distribution in the reservoir. In the direction of the maximum permeability tensor principle axis, the pressure drop is the fastest, which induces the largest aperture change in the fractures that are perpendicular to this direction.Figure 21
Oil pressure profiles when (a)kx/ky=1; (b)kx/ky=4; (c)kx/ky=1/4.
(a)(b)(c)Since the oil pressures are used to calculate the stress, strain, and displacement in the reservoir, the permeability anisotropy directly influences the results of geomechanics simulation. The distributions of mean effective stress in the reservoir using the three permeability tensors are compared, as shown in Figure22. Stress is larger in the vicinity of fractures, which can be attributed to both the larger pressure drop around the fractures and the stiffness reduction of a fractured rock compared to intact rock.Figure 22
Mean effective stress distributions when (a)kx/ky=1, (b)kx/ky=4, and (c)kx/ky=1/4.
(a)(b)(c)The influence of matrix permeability anisotropy on fracture deformation and the geomechanics behavior of the reservoir is eventually reflected on the performance of the production well. Figure23 shows the cumulative oil produced using three permeability tensors. For comparison, the results of uncoupled simulations are also included in the figure. At 1000 days, the cumulative production is decreased by 24.4%, 25.0%, and 25.3% for kx/ky=1, kx/ky=4, and kx/ky=1/4 in coupled simulations.Figure 23
Cumulative oil production histories.


## 6.3. Example 3: Water Flooding in a Five-Spot Well System

A five-spot well system water flooding case is considered. Figure24 shows the configuration of a synthetic 2D reservoir. The primary natural fractures (marked in red) are modified from reported seismic survey data [44], and the secondary natural fractures (marked in blue) have random lengths and orientations. The parameters of primary and secondary fractures are given in Table 3. The zero-stress aperture of primary and secondary fractures are 0.3 mm and 0.15 mm. All the fractures are assumed to penetrate the entire depth of the reservoir.Figure 24
(a) Configuration of the reservoir. (b) Grids in the reservoir.
(a)(b)Table 3
Fracture parameters.
ParameterPrimary fracturesSecondary fracturesJoint roughness coefficient, JRC15.07.5Joint compressive strength, JCS80.0 MPa60.0 MPaFailure ratio,Rf0.70.75Stiffness exponent,nj0.80.65Fracture residue angle,φr30.010.0Zero-stress aperture,w00.3 mm0.15 mmThe reservoir is 500 m in length, 500 m in width, and 40 m in depth. The porosity of the reservoir is 0.20, and the initial oil saturation in the reservoir is 0.8. The initial pressure in the reservoir is 15.0 MPa. In this example, an isotropic initial stress field is applied. The initial stress is 4.0 MPa in thex-direction and y-direction and 19.0 MPa in the z-direction.Four cases with different permeability tensors are investigated. The zero-stress permeabilities tensors for four cases are given in Table4. The four tensors have the same equivalent permeability of 10.0 mD. Under the initial stress, the components of permeability tensors are reduced, and the equivalent permeability is reduced to 8.23 mD. Figure 25 shows the permeability tensors under initial stress.Table 4
Permeability tensors (zero-stress).
kxx (mD)kyy (mD)kxy (mD)ke (mD)Case 15.020.00.010.0Case 220.05.00.010.0Case 312.512.57.510.0Case 412.512.5-7.510.0Figure 25
Initial permeability tensors for four cases.Water is injected from an injection well in the center of the reservoir at the fixed pressure of 15.0 MPa. Four wells located at four corners of the reservoir produce the fixed bottom hole pressure of 5.0 MPa. The water density is 1000.0 kg/m3, and the oil density is 800.0 kg/m3. The water viscosity is 0.9 cP, and the oil viscosity is 3.0 cP. The relative permeability curves are assumed to be linear. A closed boundary is considered in this example. Therefore, all the boundaries of the reservoir are impermeable. The geomechanical configuration in this example is the same with that of Example 1. The simulation is run for 10000 days.Figure26 shows the oil saturation profiles for four anisotropic cases after 10000 days of simulation, which clearly demonstrate the effect of matrix permeability anisotropic on the water displacement process. A two-dimensional permeability tensor has two principle axes, and the direction of the maximum permeability tensor principle axis becomes a favored direction for water flooding in anisotropic reservoirs.Figure 26
Oil saturation profiles for (a) Case 1, (b) Case 2, (c) Case 3, and (d) Case 4 after 10000 days of production.
(a)(b)(c)(d)The saturation distribution shows that the existence of fractures causes strong heterogeneity across the reservoir. Due to the aperture difference between primary and secondary fractures, the primary fractures have larger impacts on the oil saturation distribution. The primary fractures become major flow channels in the reservoir, and the secondary fractures increase the heterogeneity of the reservoir. Fractures whose orientation aligns with the direction of the maximum permeability tensor principle axis create fast water displacement channels.Figure27 shows the comparison of the oil and water producing rate of four cases in four production wells. The differences in producing rates result both from the anisotropic matrix permeability tensors and the heterogeneity induced by fractures. Production wells in the direction of the maximum permeability tensor principle axis to the injection well have the highest producing rate, such as well 1 and well 3 in Case 3 and well 2 and well 4 in Case 4. Production wells in the direction of the maximum permeability tensor principle axis to the injection well have the highest producing rate, and those in the direction of the minimum permeability tensor principle axis to the injection well have the lowest producing rate.Figure 27
Histories of producing rate from (a) well 1, (b) well 2, (c) well 3, and (d) well 4.
(a)(b)(c)(d)Figure28 shows the oil pressure distribution after 10000 days. Pressure drop in the direction of maximum permeability tensor principle axis is the slowest, and pressure drop in the direction of minimum permeability tensor principle axis is the fastest. The symmetry of the pressure fields is broken by the fractures. The geometry of fractures also influences the overall pressure depletion in the reservoir since the preferential flow path created by fractures induces larger water inflow to compensate for the pressure drop. The average oil pressures in the reservoir are 10.19 MPa for Case 1, 10.13 MPa for Case 2, 10.28 MPa for Case 3, and 10.38 MPa for Case 4.Figure 28
Oil pressure profiles for (a) Case 1, (b) Case 2, (c) Case 3, and (d) Case 4 after 10000 days.
(a)(b)(c)(d)Due to the differences in pressure profiles, the deformation distributions in the reservoir in four cases are also different, which is illustrated in the volume strain profiles shown in Figure29. The average volume strains are 6.41E – 4 for Case 1, 6.49E – 4 for Case 2, 6.30E – 4 for Case 3, and 6.17E − 4 for Case 4.Figure 29
Volume strain distributions for (a) Case 1, (b) Case 2, (c) Case 3, and (d) Case 4 after 10000 days.
(a)(b)(c)(d)


## 7. Conclusions

A coupled geomechanics and flow simulation approach is proposed for modeling fractured reservoirs considering matrix permeability anisotropy. The connections between matrix elements are established using the multipoint flux approximation (MPFA) to adapt the matrix permeability anisotropy. Simulation of the fluid exchange between the matrix and fractures uses the integrally embedded discrete fracture model (IEDFM) to account for the anisotropic permeability of the rock matrix. The geomechanical model of the approach uses an equivalent continuum approach. The compliance matrix is derived through stress-strain analysis, and the stress-dependent moduli of fractures are derived using an empirical model. In the proposed approach, flow simulation and geomechanics simulation are performed sequentially using the fixed-stress split iterative coupling strategy. In the iteration loops, the matrix porosity and permeability are updated using the results of pore pressure and strain field. The fracture porosity and permeability are updated using the aperture change calculated from normal stress.Several examples are presented to demonstrate the applicability of the proposed approach for simulating the coupled geomechanics and flow system in fractured reservoirs with anisotropic rock matrix. The effect of reservoir anisotropy in coupled geomechanics and flow modeling is investigated. The numerical results indicate that the dynamic behaviors of fractures are largely orientation-related due to the initial stress anisotropy and matrix permeability anisotropy. Fractures perpendicular to the minimum initial stress direction are more fluid conductive yet more prone to deformation due to depletion, and fractures perpendicular to the maximum matrix permeability direction undergo the largest closure due to the rapid pressure drop. The distributions of oil pressure and saturation show that the presence of fractures causes strong heterogeneity in fractured reservoirs. Fractures that align with the direction of the maximum permeability tensor principle axis create fast water displacement channels in the reservoir. Besides, the matrix permeability anisotropy affects the compaction in the reservoir domain, which eventually reflects on the performance of production.



---

*Source: 1200012-2022-02-22.xml*