# Effect of Natural Fractures on Stress Evolution of Unconventional Reservoirs Using a Finite Element Method and a Fracture Continuum Method

Authors: Qiang Wang; Yongquan Hu; Jinzhou Zhao
Journal: Geofluids
 (2019)

---

## Abstract

Refracturing, temporary plugging, and infilling well design play an important role in the development of reservoirs. The prediction of stress distribution can provide the basic guiding theory for the design and implementation of these techniques. In this paper, a fully-coupled three-dimensional production model based on the finite element method (FEM) and fracture continuum method (FCM) for naturally fractured reservoirs is presented to study the effects of fluid consumption on the reservoir stress. Furthermore, the effects of natural fractures on the stress re-distribution and stress re-orientation are also studied. The model also considers the influence of natural fractures on the permeability, and the effect of the effective stress on natural fracture openings, pore-elastic deformation, and fluid consumption. An analytical solution model and Eclipse were used for the comparison, which verifies the accuracy of the model results. Based on two cases of one cluster of fractures and three clusters of non-planar fractures, the research results revealed that natural fractures have a significant influence on the surrounding drainage, stress distribution, and stress re-orientation during the development. Under the influence of natural fractures, the production of the fluid along the direction of natural fractures is significantly easier, and it is highly probably that the insufficient consumption area is perpendicular to the direction of natural fractures. Compared with the conventional model, the stress distribution in the proposed model is deflected to a certain extent under the flow mode dominated by natural fractures, which is significantly prominent in the non-planar fracture model. Due to the effect of natural cracks, the absolute values of the stress, displacement, and stress difference in this model are relatively larger than those in the conventional model. Moreover, the re-orientation angles of the maximum principal stress are significantly different. After considering the natural cracks, there was an increase in the change in re-orientation and the re-orientation range. The research findings reported in this paper can be used to predict the initiation, extension, and steering process of temporary plugging fracturing fractures and refracturing fractures in fractured reservoirs.


---

## Body

## 1. Introduction

In oil and gas exploitation, the pore pressure in reservoirs gradually decreases in accordance with the development. At this point, the pore pressure gradient due to the heterogeneity of the reservoir pore pressure distribution induces a change in the initial reservoir stress distribution [1]. In conventional reservoir research, the reservoir is generally considered as a porous medium elastomer. When evaluating the coupling law between fluid flow and geological stress in the porous elastic model, the fluid flow is generally considered to occur in three main directions. However, the fluid flow is due to the pressure gradient. For reservoirs such as shale, wherein natural fractures develop and dominate the reservoir fluid flow, the fluid flow is mainly along the direction of natural fractures. In this flow mode, the conventional permeability representation cannot accurately describe the flow of fluid. Therefore, the relationship between the natural fracture permeability and permeability tensor [2], the relationship between the permeability and fracture width [3], and the relationship between the effective stress and natural fracture width were established [4]. The permeability tensor can then be used to describe the flow of the fluid in other directions, and the variation in the width of the natural fracture can be used to describe the permeability evolution during production. Hence, a three-dimensional fluid-solid coupling mathematical model in porous media was developed as an equivalent simulator for the analysis and evaluation of the stress evolution and stress steering in fractured reservoirs during depletion.Several studies have been conducted on the fluid-solid coupling in porous media; however, the effect of natural fractures on the reservoir stress was considered in a few of these studies. The study on the coupling between fluid and geomechanics in porous media can be divided into two parts. The first aspect is the coupling of fluid and geomechanics during injection [5–10]. The other aspect is the coupling of the fluid and geological stress in the depletion [11–13]. The change in the reservoir stress due to the fluid loss during hydraulic fracturing is the result of the coupling between the flow and geomechanics, which is a process that occurs within a short-time period. Ghassemi et al. established a three-dimensional fluid-solid coupling model for hydraulic fracturing and then solved the model using the finite element and displacement discontinuity methods. The distribution of the pore pressure and stress in the vicinity of the hydrofractures was obtained, and the failure of the rocks in the vicinity of the hydrofractures was analyzed using the Mohr-Coulomb failure criterion [14]. Salimzadeh et al. proposed a fully coupled three-dimensional (3D) finite element model for the hydraulic fracturing of permeable rocks and then evaluated the applicable range of classical analytic solutions in the limit case using this model [15]. Gao and Ghassemi established a pore elastic model in the process of full 3D heterogeneous hydraulic fracturing using the coupling method of fluid flow and geomechanics and then solved it using the finite element method. The induced stress distribution and steering around hydraulic fractures in heterogeneous reservoirs were confirmed as functions of time and space [16].In addition to changes in the reservoir stress during fluid injection, changes were observed when the reservoir was depleted. Oil and gas production is a significantly long-term reservoir consumption process [17–19]. Gupta et al. evaluated the influence of production consumption and completion migration on the reservoir stress, and the findings of the study indicated that the stress is more easily diverted under low-stress difference conditions. Moreover, the fracture characteristics, pressure exhaustion, and interwell interference have a significant influence on the infill well design and refracturing design in the latter stages [20]. Using a state-of-the-art coupled symmetric Galerkin boundary-element method (SGBEM) and the finite element method (FEM) model, Roussel and Sharma simulated the influence of prior production from the offset wells on the propagation direction of fractures initiated in an infill horizontal well. A stress reorientation zone was observed near the offset well after several years of production, which exhibited a significant influence on the fracture extension in the infill well [11]. Safari and Ghassemi proposed a numerical evaluation model for the changes in stress due to the depletion of the infill wells and then analyzed the bending of hydraulic fractures in the disturbed stress field and the connectivity of the hydraulic fractures between wells [21]. The failure effect and critical parameters that determine the fracture bending were systematically evaluated [21, 22]. The results revealed that the depletion perturbed the stress tensor of the stratum around the fractured horizontal well. The perturbation stress field is related to stress/formation anisotropy, fluid fluidity, pore pressure, bottom hole working pressure, and Biot’s constant. Sangnimnuan et al. evaluated the influence of the hydrofracture morphology on the stress evolution in the production process using the embedded discrete fracture method (EDFM) and the finite volume method (FVM) [1]. The study revealed that the fracture morphology has a significant influence on the stress distribution, surrounding drainage, and stress reorientation. Moreover, with a decrease in the stress difference, the steering is easier to achieve. The propagation path and radius of curvature of the hydrofractures can be predicted using the fracture propagation model combined with the varying stress field. In this paper, reservoir simulation was integrated with rock geomechanics to predict the well poststimulation productivities [23]. Several numerical and field studies have been presented to study the impact of the reservoir depletion on the pore pressure and the stress variation in the infill well zones [24]. Rezaei et al. presented refracturing analysis using a 2D poroelastic plain strain displacement discontinuity model built on the work of Zhang, Ghassemi and Zhang, and Chun [25–28]. Kumar and Ghassemi presented a geomechanical perspective to better understand the problem of “frac-hits” in horizontal well refracturing and to design solutions for it using geomechanics analysis and modeling [13]. The numerical analysis was based on a fully coupled 3D model “GeoFrac3D” with the capabilities to simulate multistage fracturing of multiple horizontal wells.However, for reservoirs with natural fractures, the evolution of the reservoir-induced stress is influenced by the pressure consumption, hydrofracture morphology, and natural fractures, among others. Moreover, the influence of natural fractures is significant. Most researchers know that natural fractures have a significant impact on production during oil and gas production. However, when natural fractures and fluid-solid coupling are considered at the same time, the treatment of natural fractures becomes a difficult task. In summary, no research was found on the effect of the natural fracture heterogeneity on the induced stress in reservoir production. This paper summarizes and presents the results of the analysis conducted on the effect of natural fractures on the stress and stress diversion under the condition of plane fractures and nonplane fractures. This study was a combination and expansion of the fluid flow/geomechanics coupling theory and natural fracture heterogeneity in the pore elastic model. Compared with previous studies, our model is closer to the actual geological situation of the reservoir, so the results are more reliable and accurate. The research findings reported in this paper can be used to predict the initiation, extension, and steering process of temporary plugging fracturing fractures and refracturing fractures in fractured reservoirs.


## 2. Mathematic Model

### 2.1. Natural Fracture Continuum Method

In naturally fractured reservoirs such as shale and tight sandstone, the permeability of the reservoir matrix is significantly small and it is difficult to extract the fluid based only on the matrix [23, 29]. In engineering, multicluster fracturing is generally employed in horizontal wells to activate reservoirs with natural fractures. The natural fractures then become the main channel of the oil and gas production and then dominate the flow mode of the reservoir fluid. Given that the permeability of natural fractures is significantly greater than that of the matrix, the contribution of the matrix permeability to fluid flow is neglected in this model and only the influence of natural fractures is considered [29]. The distribution angle of natural fractures developed in shale reservoirs is typically consistent, i.e., the inclination angle and approximation angle are generally distributed within a certain range of the main direction [30]. Therefore, it was assumed that the approximation angle and inclination angle of all-natural fractures in the model are consistent. In this model, the fracture continuum method proposed by McKenna and Reeves is used to describe the characteristics of natural fractures. The permeability of the element grid is characterized by the natural fracture inclination angle, approximation angle, opening angle, and fracture spacing [31]. The advantage of this method is that the natural fracture grid cannot be processed separately, and the results were obtained under the condition that the element is sufficiently small to meet the precision requirements [3, 32]. Based on previous research [3], the natural fracture permeability was dispersed to the continuous unit. The permeability tensor of natural fractures defines the permeability for each grid cell in the model domain. For one fracture set, the permeability tensor can be expressed as follows:
(1)k=knfn22+n32−n1n2−n3n1−n1n2n22+n32−n2n3−n3n1−n2n3n22+n32,where k is the permeability tensor; n1, n2, and n3 are the units normal to the fracture plane in the x, y, and z directions, respectively; and knf is the permeability of the natural fracture. Moreover, n1, n2, and n3 can be expressed as follows:
(2)n1=cosςπ180sinξπ180,n2=cosςπ180cosξπ180,n3=−sinςπ180,where ς=90°−dip and ξ=strike−90°. Using Equations (1) and (2), the permeability tensor can be used to represent the heterogeneity caused by natural fractures. This heterogeneity is determined by the opening, spacing, inclination, and approach angle of natural fractures. In the model, the fluid-solid coupling of the depletion is considered and the stress distribution of the reservoir changes with respect to time, thus leading to changes in the natural fracture permeability. The natural fracture permeability is mainly determined by the fracture opening as follows [3]:
(3)knf=δ312d,where δ is the reduction of the fracture aperture (fracture closure) and d is the fracture spacing.To accurately describe the natural fracture permeability as a function of time, a coupling model of the natural fracture opening and effective stress was introduced. Bandis et al. evaluated the crack opening and effective normal stress and then proposed a semi-logarithmic relationship [33]:
(4)logσn=ω+χδ,where σn is the normal stress on the crack surface and ω and χ represent the constants in the semilog fracture closure/stress relationship.The effective normal stress that acts on the crack surface can be expressed as follows:(5)σn′=12σH′+σh′+12σH′−σh′cos2θ,where θ is the angle between the normal to the fracture and the σH direction and σH′ and σh′, respectively, represent the maximum effective principal stress and the minimum effective principal stress, MPa.By combining Equations (4) and (5), the following was obtained [4]:
(6)δ=1χlogn+12n+n−12ncos2θ+logσH′χ−ωχ,where n is the effective stress ratio, which can be expressed as n=σH′/σh′.


### 2.2. Mechanism of Fluid-Solid Coupling Mechanics

#### 2.2.1. Principle of Effective Stress

Rock is a porous medium composed of a solid skeleton and pores. In addition, there are fluids (oil, gas, and water) stored in the pores of the skeleton. The fluid in the rock can bear or transfer pressure, which is defined as the pore pressure. Moreover, the stress transmitted through the contact surface between the rock particles is the effective stress. The Biot effective stress can be expressed as follows [34, 35]:
(7)σij′=σij+αδijpf,where σij′ is the effective stress, σij is the total stress, δij is the Kronecker symbol, and α represents the Biot constant, which can be defined as follows [36]:
(8)α=1−KvKs,where Kv and Ks represent the volume compression modulus of the rock and the compression modulus of the solid particles, respectively. Moreover, α is typically close to 1; thus, α=1 was used.


#### 2.2.2. Stress Balance Equation

The theory of saturated-unsaturated seepage evaluates the movement of two or more fluids through pores. This study addressed the flow of a single-phase fluid in a porous medium. Based on the principle of virtual work, it can be known that at a certain moment, the virtual work of the rock mass of a unit is equal to the virtual work generated by the force that acts on the unit (physical force and surface force); thus,(9)∫VδεTdσdV−∫VδuTdfdV−∫SδuTdFds=0,where F is the surface force, f is the body force, and δε and δu are the virtual strain and virtual displacement, respectively.The constitutive relationship expressed in an incremental form is(10)dσ′=Depdε−dεl,where Dep is the elastoplastic matrix and dεl is the particle compression due to the pore fluid pressure, which only deforms in the positive direction without shear direction. The specific expression is as follows:
(11)dεl=−mdpf3Ks,where m=1,1,1,0,0,0T.Formula (7) can be expressed as follows:
(12)σ′=σ+αmpf.Forα=1, the combination of Formula (12) and the sum Formula (9) can be expressed as follows:
(13)∫VδεTDepdε+mdpf3KsdV−∫VδεTmdpfdV−∫VδuTdfdV−∫SδuTdfdV=0.Given that the seepage continuity equation contains a time-related term, to couple the stress with the seepage, the time derivative of the virtual work Equation (13) is required. The specific expression is as follows:
(14)∫VδεTDep∂ε∂t+m13Ks∂pf∂tdV−∫VδεTm∂pf∂tdV−∫VδεTm∂pf∂tdV−∫SδuT∂F∂tdS=0.


#### 2.2.3. Continuity Equation

By combining Darcy’s law with the conservation of fluid mass, a matrix flow model can be established to characterize the flow of porous media. The amount of fluid that flows into this volume in timedt should be equal to the increase in its internal water storage. Through derivation, the continuity equation of seepage flow can be obtained as follows:
(15)ϕcf+α−ϕKs∂pf∂t+α∂ε∂t−∇kμf∇pf+ρfg=0,where g is the acceleration vector of gravity, cf is the compressibility of the fluid, and ϕ is the matrix porosity.In the process of oil and gas development, hydraulic fracturing is generally used to generate hydrofractures to increase the productive drainage area. As shown in the Figure1, a hydrofracture was located at the y-z section. Hence, the continuity equation can be expressed as follows:
(16)ϕcf+α−ϕKs∂pf∂t+α∂ε∂t−∇kμf∇pf+ρfg+qΓ=0,where qΓ represents the interchange term between the fracture and matrix.Figure 1
Single-phase solution validation model.


## 2.1. Natural Fracture Continuum Method

In naturally fractured reservoirs such as shale and tight sandstone, the permeability of the reservoir matrix is significantly small and it is difficult to extract the fluid based only on the matrix [23, 29]. In engineering, multicluster fracturing is generally employed in horizontal wells to activate reservoirs with natural fractures. The natural fractures then become the main channel of the oil and gas production and then dominate the flow mode of the reservoir fluid. Given that the permeability of natural fractures is significantly greater than that of the matrix, the contribution of the matrix permeability to fluid flow is neglected in this model and only the influence of natural fractures is considered [29]. The distribution angle of natural fractures developed in shale reservoirs is typically consistent, i.e., the inclination angle and approximation angle are generally distributed within a certain range of the main direction [30]. Therefore, it was assumed that the approximation angle and inclination angle of all-natural fractures in the model are consistent. In this model, the fracture continuum method proposed by McKenna and Reeves is used to describe the characteristics of natural fractures. The permeability of the element grid is characterized by the natural fracture inclination angle, approximation angle, opening angle, and fracture spacing [31]. The advantage of this method is that the natural fracture grid cannot be processed separately, and the results were obtained under the condition that the element is sufficiently small to meet the precision requirements [3, 32]. Based on previous research [3], the natural fracture permeability was dispersed to the continuous unit. The permeability tensor of natural fractures defines the permeability for each grid cell in the model domain. For one fracture set, the permeability tensor can be expressed as follows:
(1)k=knfn22+n32−n1n2−n3n1−n1n2n22+n32−n2n3−n3n1−n2n3n22+n32,where k is the permeability tensor; n1, n2, and n3 are the units normal to the fracture plane in the x, y, and z directions, respectively; and knf is the permeability of the natural fracture. Moreover, n1, n2, and n3 can be expressed as follows:
(2)n1=cosςπ180sinξπ180,n2=cosςπ180cosξπ180,n3=−sinςπ180,where ς=90°−dip and ξ=strike−90°. Using Equations (1) and (2), the permeability tensor can be used to represent the heterogeneity caused by natural fractures. This heterogeneity is determined by the opening, spacing, inclination, and approach angle of natural fractures. In the model, the fluid-solid coupling of the depletion is considered and the stress distribution of the reservoir changes with respect to time, thus leading to changes in the natural fracture permeability. The natural fracture permeability is mainly determined by the fracture opening as follows [3]:
(3)knf=δ312d,where δ is the reduction of the fracture aperture (fracture closure) and d is the fracture spacing.To accurately describe the natural fracture permeability as a function of time, a coupling model of the natural fracture opening and effective stress was introduced. Bandis et al. evaluated the crack opening and effective normal stress and then proposed a semi-logarithmic relationship [33]:
(4)logσn=ω+χδ,where σn is the normal stress on the crack surface and ω and χ represent the constants in the semilog fracture closure/stress relationship.The effective normal stress that acts on the crack surface can be expressed as follows:(5)σn′=12σH′+σh′+12σH′−σh′cos2θ,where θ is the angle between the normal to the fracture and the σH direction and σH′ and σh′, respectively, represent the maximum effective principal stress and the minimum effective principal stress, MPa.By combining Equations (4) and (5), the following was obtained [4]:
(6)δ=1χlogn+12n+n−12ncos2θ+logσH′χ−ωχ,where n is the effective stress ratio, which can be expressed as n=σH′/σh′.


## 2.2. Mechanism of Fluid-Solid Coupling Mechanics

### 2.2.1. Principle of Effective Stress

Rock is a porous medium composed of a solid skeleton and pores. In addition, there are fluids (oil, gas, and water) stored in the pores of the skeleton. The fluid in the rock can bear or transfer pressure, which is defined as the pore pressure. Moreover, the stress transmitted through the contact surface between the rock particles is the effective stress. The Biot effective stress can be expressed as follows [34, 35]:
(7)σij′=σij+αδijpf,where σij′ is the effective stress, σij is the total stress, δij is the Kronecker symbol, and α represents the Biot constant, which can be defined as follows [36]:
(8)α=1−KvKs,where Kv and Ks represent the volume compression modulus of the rock and the compression modulus of the solid particles, respectively. Moreover, α is typically close to 1; thus, α=1 was used.


### 2.2.2. Stress Balance Equation

The theory of saturated-unsaturated seepage evaluates the movement of two or more fluids through pores. This study addressed the flow of a single-phase fluid in a porous medium. Based on the principle of virtual work, it can be known that at a certain moment, the virtual work of the rock mass of a unit is equal to the virtual work generated by the force that acts on the unit (physical force and surface force); thus,(9)∫VδεTdσdV−∫VδuTdfdV−∫SδuTdFds=0,where F is the surface force, f is the body force, and δε and δu are the virtual strain and virtual displacement, respectively.The constitutive relationship expressed in an incremental form is(10)dσ′=Depdε−dεl,where Dep is the elastoplastic matrix and dεl is the particle compression due to the pore fluid pressure, which only deforms in the positive direction without shear direction. The specific expression is as follows:
(11)dεl=−mdpf3Ks,where m=1,1,1,0,0,0T.Formula (7) can be expressed as follows:
(12)σ′=σ+αmpf.Forα=1, the combination of Formula (12) and the sum Formula (9) can be expressed as follows:
(13)∫VδεTDepdε+mdpf3KsdV−∫VδεTmdpfdV−∫VδuTdfdV−∫SδuTdfdV=0.Given that the seepage continuity equation contains a time-related term, to couple the stress with the seepage, the time derivative of the virtual work Equation (13) is required. The specific expression is as follows:
(14)∫VδεTDep∂ε∂t+m13Ks∂pf∂tdV−∫VδεTm∂pf∂tdV−∫VδεTm∂pf∂tdV−∫SδuT∂F∂tdS=0.


### 2.2.3. Continuity Equation

By combining Darcy’s law with the conservation of fluid mass, a matrix flow model can be established to characterize the flow of porous media. The amount of fluid that flows into this volume in timedt should be equal to the increase in its internal water storage. Through derivation, the continuity equation of seepage flow can be obtained as follows:
(15)ϕcf+α−ϕKs∂pf∂t+α∂ε∂t−∇kμf∇pf+ρfg=0,where g is the acceleration vector of gravity, cf is the compressibility of the fluid, and ϕ is the matrix porosity.In the process of oil and gas development, hydraulic fracturing is generally used to generate hydrofractures to increase the productive drainage area. As shown in the Figure1, a hydrofracture was located at the y-z section. Hence, the continuity equation can be expressed as follows:
(16)ϕcf+α−ϕKs∂pf∂t+α∂ε∂t−∇kμf∇pf+ρfg+qΓ=0,where qΓ represents the interchange term between the fracture and matrix.Figure 1
Single-phase solution validation model.


## 2.2.1. Principle of Effective Stress

Rock is a porous medium composed of a solid skeleton and pores. In addition, there are fluids (oil, gas, and water) stored in the pores of the skeleton. The fluid in the rock can bear or transfer pressure, which is defined as the pore pressure. Moreover, the stress transmitted through the contact surface between the rock particles is the effective stress. The Biot effective stress can be expressed as follows [34, 35]:
(7)σij′=σij+αδijpf,where σij′ is the effective stress, σij is the total stress, δij is the Kronecker symbol, and α represents the Biot constant, which can be defined as follows [36]:
(8)α=1−KvKs,where Kv and Ks represent the volume compression modulus of the rock and the compression modulus of the solid particles, respectively. Moreover, α is typically close to 1; thus, α=1 was used.


## 2.2.2. Stress Balance Equation

The theory of saturated-unsaturated seepage evaluates the movement of two or more fluids through pores. This study addressed the flow of a single-phase fluid in a porous medium. Based on the principle of virtual work, it can be known that at a certain moment, the virtual work of the rock mass of a unit is equal to the virtual work generated by the force that acts on the unit (physical force and surface force); thus,(9)∫VδεTdσdV−∫VδuTdfdV−∫SδuTdFds=0,where F is the surface force, f is the body force, and δε and δu are the virtual strain and virtual displacement, respectively.The constitutive relationship expressed in an incremental form is(10)dσ′=Depdε−dεl,where Dep is the elastoplastic matrix and dεl is the particle compression due to the pore fluid pressure, which only deforms in the positive direction without shear direction. The specific expression is as follows:
(11)dεl=−mdpf3Ks,where m=1,1,1,0,0,0T.Formula (7) can be expressed as follows:
(12)σ′=σ+αmpf.Forα=1, the combination of Formula (12) and the sum Formula (9) can be expressed as follows:
(13)∫VδεTDepdε+mdpf3KsdV−∫VδεTmdpfdV−∫VδuTdfdV−∫SδuTdfdV=0.Given that the seepage continuity equation contains a time-related term, to couple the stress with the seepage, the time derivative of the virtual work Equation (13) is required. The specific expression is as follows:
(14)∫VδεTDep∂ε∂t+m13Ks∂pf∂tdV−∫VδεTm∂pf∂tdV−∫VδεTm∂pf∂tdV−∫SδuT∂F∂tdS=0.


## 2.2.3. Continuity Equation

By combining Darcy’s law with the conservation of fluid mass, a matrix flow model can be established to characterize the flow of porous media. The amount of fluid that flows into this volume in timedt should be equal to the increase in its internal water storage. Through derivation, the continuity equation of seepage flow can be obtained as follows:
(15)ϕcf+α−ϕKs∂pf∂t+α∂ε∂t−∇kμf∇pf+ρfg=0,where g is the acceleration vector of gravity, cf is the compressibility of the fluid, and ϕ is the matrix porosity.In the process of oil and gas development, hydraulic fracturing is generally used to generate hydrofractures to increase the productive drainage area. As shown in the Figure1, a hydrofracture was located at the y-z section. Hence, the continuity equation can be expressed as follows:
(16)ϕcf+α−ϕKs∂pf∂t+α∂ε∂t−∇kμf∇pf+ρfg+qΓ=0,where qΓ represents the interchange term between the fracture and matrix.Figure 1
Single-phase solution validation model.


## 3. Model Calculation and Validation

### 3.1. Finite Element Dispersion

The shape function can be defined as follows:(17)u=Nuu,¯ε=Bu,¯pf=Npp¯f,where B and Nu represent the shape function of the displacement, Np is the shape function of pore pressure, u¯ is the nodal displacement, and pw¯ is the pore pressure of the element node.By substituting Equation (17) into Equation (15), the solid phase finite element formula can be obtained after simplification, as follows:
(18)Kdu¯dt+Cdp¯wdt=dfdt,where
(19)K=∫VBTDepBdV,C=∫VBTDepm3KsNpdV−∫VBTmNpdV,df=∫VNuTdfdV+∫SNuTdtdS.In the analysis of the seepage field, there are two types of boundary conditions. The first is the flow boundary condition, and the second is the pore pressure boundary condition. The flow boundary conditions can be expressed as follows:(20)qΓ=kμf∂pf∂nc,where nc is the normal unit of the flow boundary. This boundary condition can be used as the outer boundary condition. When there are production fractures in the reservoir after fracturing, the flow boundary condition can also be considered as the internal boundary condition or fracture boundary condition.The pore pressure boundary conditions can be expressed as follows:(21)pf=pfb,where pfb is the pore pressure value at the known boundary.The Galerkin method can be used, which can be expressed as follows:(22)∫VaTA¯dV+∫SbTB¯dS=0,where a and b are arbitrary functions, A¯ is the governing equation, and B¯ is the continuity equation through the boundary.By substituting Equation (15) as A¯, Equation (20) as B¯, and the formal function expression (17) into Equation (22) and by setting a=−b, after the simplification, the following can be obtained:
(23)Edu¯dt+Fp¯w+Gdp¯wdt=f^,where
(24)E=∫VNpTαBdV,F=∫V∇NpTkμf∇NpdV,G=∫VNpTϕcf+α−ϕKsNpdV,f^=∫sNpTqΓdS−∫V∇NpTkμfρfgdV.By combining Equation (18) with Equation (23), the geomechanic-seepage coupling equation can be expressed as follows:
(25)KCEGddtu¯p¯w+000Fu¯p¯w=dfdtf^.


### 3.2. Model Validation

#### 3.2.1. Single-Phase Flow Solver

To verify the oil-water flow model, a single fracture seepage model was established (as shown in Figure1) and the same model was established by Eclipse software for a simulation comparison. The permeability, porosity, and reservoir thickness of the model were 5 md, 0.1, and 50 m, respectively. The boundary of the model was closed; thus, the reserves in the model were constant. The other main parameters were the same as those shown in Table 1. In the proposed model, the flow model was established without considering natural fractures. Fluid-solid coupling is not considered in Eclipse software. The daily output and cumulative oil production after one year of production were obtained by adopting the method of constant pressure production, as shown in Figure 2. The results of the model were approximately consistent with those of Eclipse. Due to the different factors considered by the models, the change process was slightly different. Therefore, the model in this paper demonstrated a high accuracy.Table 1
Primary parameters of calculation.
Parameter
Value
Unit
Young’s modulus
35.5
GPa
Poisson’s ratio
0.2
—
Biot coefficient
1
—
Initial porosity
0.25
—
Fluid compressibility
8×10−4
MPa-1
Fluid viscosity
1.1
MPa s
The principal stress inx-direction
50
MPa
The principal stress iny-direction
54
MPa
The principal stress inz-direction
57
MPa
The initial pore pressure
30
MPa
The initial opening of natural fracture
0.2
mm
Average natural fracture spacing
20
m
Natural fracture dip
80
°
Natural fracture strike
40
°
Well radius
0.0762
mFigure 2
Flow rate comparison between the proposed model and Eclipse for one year of production.


## 3.1. Finite Element Dispersion

The shape function can be defined as follows:(17)u=Nuu,¯ε=Bu,¯pf=Npp¯f,where B and Nu represent the shape function of the displacement, Np is the shape function of pore pressure, u¯ is the nodal displacement, and pw¯ is the pore pressure of the element node.By substituting Equation (17) into Equation (15), the solid phase finite element formula can be obtained after simplification, as follows:
(18)Kdu¯dt+Cdp¯wdt=dfdt,where
(19)K=∫VBTDepBdV,C=∫VBTDepm3KsNpdV−∫VBTmNpdV,df=∫VNuTdfdV+∫SNuTdtdS.In the analysis of the seepage field, there are two types of boundary conditions. The first is the flow boundary condition, and the second is the pore pressure boundary condition. The flow boundary conditions can be expressed as follows:(20)qΓ=kμf∂pf∂nc,where nc is the normal unit of the flow boundary. This boundary condition can be used as the outer boundary condition. When there are production fractures in the reservoir after fracturing, the flow boundary condition can also be considered as the internal boundary condition or fracture boundary condition.The pore pressure boundary conditions can be expressed as follows:(21)pf=pfb,where pfb is the pore pressure value at the known boundary.The Galerkin method can be used, which can be expressed as follows:(22)∫VaTA¯dV+∫SbTB¯dS=0,where a and b are arbitrary functions, A¯ is the governing equation, and B¯ is the continuity equation through the boundary.By substituting Equation (15) as A¯, Equation (20) as B¯, and the formal function expression (17) into Equation (22) and by setting a=−b, after the simplification, the following can be obtained:
(23)Edu¯dt+Fp¯w+Gdp¯wdt=f^,where
(24)E=∫VNpTαBdV,F=∫V∇NpTkμf∇NpdV,G=∫VNpTϕcf+α−ϕKsNpdV,f^=∫sNpTqΓdS−∫V∇NpTkμfρfgdV.By combining Equation (18) with Equation (23), the geomechanic-seepage coupling equation can be expressed as follows:
(25)KCEGddtu¯p¯w+000Fu¯p¯w=dfdtf^.


## 3.2. Model Validation

### 3.2.1. Single-Phase Flow Solver

To verify the oil-water flow model, a single fracture seepage model was established (as shown in Figure1) and the same model was established by Eclipse software for a simulation comparison. The permeability, porosity, and reservoir thickness of the model were 5 md, 0.1, and 50 m, respectively. The boundary of the model was closed; thus, the reserves in the model were constant. The other main parameters were the same as those shown in Table 1. In the proposed model, the flow model was established without considering natural fractures. Fluid-solid coupling is not considered in Eclipse software. The daily output and cumulative oil production after one year of production were obtained by adopting the method of constant pressure production, as shown in Figure 2. The results of the model were approximately consistent with those of Eclipse. Due to the different factors considered by the models, the change process was slightly different. Therefore, the model in this paper demonstrated a high accuracy.Table 1
Primary parameters of calculation.
Parameter
Value
Unit
Young’s modulus
35.5
GPa
Poisson’s ratio
0.2
—
Biot coefficient
1
—
Initial porosity
0.25
—
Fluid compressibility
8×10−4
MPa-1
Fluid viscosity
1.1
MPa s
The principal stress inx-direction
50
MPa
The principal stress iny-direction
54
MPa
The principal stress inz-direction
57
MPa
The initial pore pressure
30
MPa
The initial opening of natural fracture
0.2
mm
Average natural fracture spacing
20
m
Natural fracture dip
80
°
Natural fracture strike
40
°
Well radius
0.0762
mFigure 2
Flow rate comparison between the proposed model and Eclipse for one year of production.


## 3.2.1. Single-Phase Flow Solver

To verify the oil-water flow model, a single fracture seepage model was established (as shown in Figure1) and the same model was established by Eclipse software for a simulation comparison. The permeability, porosity, and reservoir thickness of the model were 5 md, 0.1, and 50 m, respectively. The boundary of the model was closed; thus, the reserves in the model were constant. The other main parameters were the same as those shown in Table 1. In the proposed model, the flow model was established without considering natural fractures. Fluid-solid coupling is not considered in Eclipse software. The daily output and cumulative oil production after one year of production were obtained by adopting the method of constant pressure production, as shown in Figure 2. The results of the model were approximately consistent with those of Eclipse. Due to the different factors considered by the models, the change process was slightly different. Therefore, the model in this paper demonstrated a high accuracy.Table 1
Primary parameters of calculation.
Parameter
Value
Unit
Young’s modulus
35.5
GPa
Poisson’s ratio
0.2
—
Biot coefficient
1
—
Initial porosity
0.25
—
Fluid compressibility
8×10−4
MPa-1
Fluid viscosity
1.1
MPa s
The principal stress inx-direction
50
MPa
The principal stress iny-direction
54
MPa
The principal stress inz-direction
57
MPa
The initial pore pressure
30
MPa
The initial opening of natural fracture
0.2
mm
Average natural fracture spacing
20
m
Natural fracture dip
80
°
Natural fracture strike
40
°
Well radius
0.0762
mFigure 2
Flow rate comparison between the proposed model and Eclipse for one year of production.


## 4. Case Studies

### 4.1. Basic Parameters

As discussed in this section, a three-dimensional reservoir physical model was established to evaluate the influence of anisotropy due to natural fractures on the reservoir stress during the oil and gas development. To fully demonstrate the significant influence of natural fractures on the generation of induced stress, a model was established that does not consider the influence ofkxz, kyz, and kxy. However, the influence of the permeability in the main direction was considered for comparison. As shown in Figure 3, the length, width, and height of the model were 400 m, 400 m, and 100 m, respectively, and the hydraulic fracture was located at the cross section of the y-z-axis. The domain was discretized to 50 four-node cells in the x- and y-directions and 20 four-node cells in the z-direction. For the treatment of the reservoir boundary, the nodes on the left and right surfaces of the model were free to move along the x-axis and not the y-axis. All the nodes at the top and bottom outer surfaces were free to move along the y-axis and not the x-axis. Only the vertical displacement of all the nodes on the outer surface of the top and bottom could be achieved. Moreover, the six boundaries were nonflow boundaries, i.e., the pore pressure in the reservoir was continuously reduced in the process of production. In all the calculations below, the downhole flow velocity was calculated using the Peaceman equation [37]. The other parameters in the calculation are shown in Table 1. Figures 4(a) and 4(b) present the locations of the single fractures and three clusters of nonplanar fractures, respectively.Figure 3
A schematic diagram of the 3D physical model.Figure 4
A schematic diagram of the hydrofracture layout ((a) presents the planar single fracture, and (b) presents the nonplanar multiple fractures).


### 4.2. Single Fracture Study

Based on the model and physical parameters established in Section4.1, in this section, the focus is on the influence of natural fractures on the stress distribution and steering in the production process of a single fracture in the y-z section. To fully reflect the effect of natural fractures on the stress, a conventional model was established without considering the effect of natural fractures, for comparison. Compared with the boundary conditions in the model, the calculated parameters were consistent with the proposed model.During the production, the pore pressure decreases continuously due to the continuous fluid extraction and the pressure is the lowest near the fracture. Figure5 compares the distribution of the reservoir pore pressure after 5 years of production with and without considering the influence of natural fractures. As can be seen from the figure, due to the influence of natural fractures, there was a change in the fluid flow pattern. Compared with the conventional model, the pore pressure dissipation area was greater and it was easier to reduce the pore pressure along the direction of natural fractures for the formation of a clockwise deflecting pressure drop funnel. Therefore, under the natural fracture distribution condition, the fluid in the upper right and lower left corners could be easily extracted, whereas the fluid in the upper left and lower right corners could not be easily extracted. This can serve as an essential guide to the design of infill wells.Figure 5
Panels (a) and (b), respectively, represent the pore pressure distribution without the influence of the natural fractures and the pore pressure under the influence of the natural fractures.
(a)
(b)According to the principle of effective stress, the reservoir stress changes under the influence of the constant dissipation of the pore pressure. For the conventional model, the distribution of the stressσxx, σyy, and σzz in the reservoir after 5 years of production is shown in Figures 6(a), 6(b), and 6(c). Due to fluid consumption, the reservoir stress distribution exhibited a high complexity. Moreover, there was a decrease in σxx near the crack as a result of the pressure dissipation, whereas σxx increased on the left and right parts of the domain to support the pressure depletion in the x-direction. In addition, there was a decrease in σyy near the crack as a result of the pressure dissipation, whereas σyy increased at the top and bottom parts of the domain to support the pressure depletion in the y-direction. Moreover, σzz mainly exhibited a decrease in stress in the central region of the domain. Considering the influence of natural fractures on the stress evolution, the results are shown in Figures 7(a), 7(b), and 7(c). Compared with the conventional model, the stress distribution of the proposed model is similar to that of the conventional model. However, considering the influence of tensors, the pressure dissipation of the proposed model was more rapid and the changes in stress were greater. The most significant difference is that the distributions of σxx, σyy, and σzz exhibited a degree of deflection. The stress distribution deflection was found to be significantly different from the stress distribution near the fractures in the conventional model, and this difference has a significant influence on the fracture initiation direction and fracture extension of refracturing. As can be seen in Figures 8(a) and 8(b), the distribution of the shear stress on the x-y section exhibited a clockwise deflection when compared with the conventional model, due to the influence of natural fractures. From a comparison between Figures 9(c) and 9(d) and Figures 9(a) and 9(b), the distribution of the displacement was found to be deflected to a certain extent under the influence of natural fractures. Figures 10(a) and 10(b) present the displacement profiles on Paths A and B, respectively. In all the figures, the solid lines represent the proposed model and the dotted lines represent the conventional model. The displacement trend of the two models was similar. The main difference was that there was a more significant change in the absolute value of the displacement after the influence of natural fractures was considered.Figure 6
(a), (b), and (c) present the principal stress in thex-direction, y-direction, and z-direction without considering the influence of the permeability tensor of natural fractures.
(a)
(b)
(c)Figure 7
(a), (b), and (c) present the principal stress in thex-direction, y-direction, and z-direction considering the influence of the permeability tensor of natural fractures.
(a)
(b)
(c)Figure 8
(a) and (b) present the induced shear stress in thex-y section with and without considering natural fractures, respectively.
(a)
(b)Figure 9
(a) and (b) presentux and uy on the x-y section without considering natural fractures, respectively, and (c) and (d) present ux and uy on the x-y section considering natural fractures, respectively.
(a)
(b)
(c)
(d)Figure 10
(a) and (b) present the displacement profiles on Paths A and B.
(a)
(b)Figures11(a), 11(b), and 11(c) present the σxx, σyy, and σzx stress distributions of the proposed model and conventional model after 5 years of production on the z-x section, respectively. As can be seen from the figure, the stress distribution of the conventional model exhibited a symmetrical distribution. Compared with the conventional model, the overall trend of the stress distribution in the proposed model was in good agreement. However, the stress distribution was symmetrically distributed under the influence of natural fractures, and there was a relative increase in the variation range of the stress value.Figure 11
(a), (b), and (c) present the comparison of the principal stress along thex-direction, y-direction, and the induced shear stress on the z-x section, respectively.
(a)
(b)
(c)Significant attention has been directed toward the reservoir stress evolution and surrounding drainage in the process of fluid recovery, given their significance for the latter refracturing and completion of infill wells. The change in the direction of the horizontal principal stress direction has a significant influence on the design of refracturing. Due to the heterogeneity of the permeability and natural fractures, the heterogeneity consumption of the pore pressure may result in stress reorientation and changes in the direction of the horizontal maximum principal stress. As can be seen from Figures12(a) and 12(b), a region with a stress difference less than 0 was observed near the fracture after one year of production. In addition, the stress in this region was reversed. It can also be found that the time when the stress difference changes the most is not the 20 years with the longest production life but the 10th year. This is because with an increase in the production life, most of the fluid in the easily flowing region was produced, whereas the fluid in the previously relatively noneasily flowing region was gradually produced. Therefore, there was a decrease in the difference between the pore pressures in the regions, which led to a decrease in the stress differences and an increase in the influence range. Figures 12(c) and 12(d) present the shear stress profiles on Paths A and B. Compared with the conventional model, the shear stress absolute value of the proposed model is influenced more significantly by natural cracks.Figure 12
(a) and (b) are present the stress difference profiles on Paths A and B of different production years, and (c) and (d) present the shear stress profiles on Paths A and B of different production years.
(a)
(b)
(c)
(d)Refracturing is one of the most effective methods for latter reservoir development. The distribution of the reservoir stress after the production and the direction of the horizontal maximum principal stress are two critical factors that determine the fracture pressure, fracture initiation direction, and fracture extension trajectory of refracturing. The stress evolution in the production process was previously analyzed. For the solution of the change in reorientation of the maximum horizontal principal stress (σHmax) on the x-y section, the action of the vertical principal stress (σzz) was neglected; thus, the original direction of the maximum horizontal principal stress is shown in Figure 13(a). After 5 and 10 years of production, there was a significant change in the direction of σHmax. In the conventional model, the change area of the principal stress direction was in the shape of a vertical ellipse, whereas the change area of the σHmax direction in the proposed model was in the shape of an inclined ellipse. The degree of inclination is related to the location of natural fractures. The area with the largest reorientation change of σHmax was near the hydrofracture. Compared with the conventional model, the change in reorientation of σHmax in the proposed model was influenced more significantly by natural fractures (θ1′>θ1, θ2′>θ2). In the conventional model, the change in reorientation of σHmax was 90° (Figure 14(a)); however, the change in reorientation of σHmax in the proposed model was 110° under the influence of natural cracks (Figure 14(b)). Similar to the change in the stress difference, the change in the reorientation was the largest at the 10th year, not the 20th year. Moreover, the range of the reorientation areas is proportional to production time.Figure 13
(a) presents the direction of the maximum principal stress on thex-y section in the initial state. (b) and (c) present the direction of the maximum horizontal principal stress after 5 years and 10 years of production without considering natural fractures. (d) and (e) present the direction of the maximum horizontal principal stress after 5 years and 10 years of production considering the influence of natural fractures.Figure 14
(a) and (b) present the steering angle of the maximum horizontal principal stress along Path B without and with considering the influence of natural fractures in different production years, respectively.
(a)
(b)


### 4.3. Nonplanar Multifracture Study

In the actual segmented multicluster fracturing process, there is stress interference between the fractures when multiple fractures are simultaneously extended, which leads to the case wherein the middle fractures are restrained and the fractures on both sides are extended in a nonplanar manner. Hence, the displacement discontinuity method was used to simulate the extended trajectories of three clusters of hydrofractures with spacings of 50 m. The three hydraulic fractures are considered as production fractures. It was assumed that all three clusters of fractures pass through the reservoir in thez-direction. The other calculation parameters are the same as those in Table 1. Under this condition, the reservoir stress distribution after 1, 5, and 10 years of production was simulated. Figure 15 presents the pore pressure distribution after 1, 5, and 10 years of nonplanar fracture production using the conventional model. The pore pressure distribution in the domain was axially symmetric. However, under the influence of natural fractures, there was an initial decrease in the pore pressure along the strike direction of the natural fractures, and the distribution exhibited central symmetry rather than axial symmetry (as shown in Figure 16). The upper left part and the lower right part of the domain are underexploited areas under such reservoir geological conditions and distributed fractures.Figure 15
(a), (b), and (c) present the pore pressure distribution after 1, 5, and 10 years of nonplanar fracture production without considering the influence of natural fractures.
(a)
(b)
(c)Figure 16
(a), (b), and (c) present the pore pressure distribution after 1, 5, and 10 years of nonplanar fracture production considering the influence of natural fractures.
(a)
(b)
(c)Due to the influence of the nonplanar fracture geometry, the stress distribution after production exhibited a high complexity. Without considering the effect of natural fractures, the stress distribution of the conventional model was axially symmetric. Compared with a single fracture, there was an increase in the rates of the pore pressure dissipation and stress decrease and there was an increase in the size of the influence range, as shown in Figures17(a), 17(b), and 17(c). The stress (σxx) on the left and right of the domain and the stress (σyy) on the top and bottom of the domain increased to support the pressure depletion of the middle area. When the influence of natural fractures was considered, the stress distribution in the proposed model exhibited a significant degree of deflection, as shown in Figures 18(a), 18(b), and 18(c). The deflection of the stress distribution confirms that natural fractures have a significant influence on the reservoir fluid flow and production.Figure 17
(a), (b), and (c) present theσx, σy, and σz stress distributions after 10 years of nonplanar fracture production without considering the influence of natural fractures.
(a)
(b)
(c)Figure 18
(a), (b), and (c) present theσx, σy, σz stress distribution after 10 years of nonplanar fracture production considering the influence of natural fractures.
(a)
(b)
(c)Compared with a single crack, the distribution of the principal stress was more complex. The change in reorientation ofσHmax was θ1′>θ1, θ2′>θ2, as shown in Figure 19. Due to the influence of the nonplanar crack geometry, the rotation angle of the maximum principal stress direction in the conventional model reached 135°. In the proposed model, the rotation angle of the principal stress direction near the crack was concentrated at 140°. The reorientation region affected by the natural fractures was larger than those of the conventional models. Therefore, under the action of the nonplanar fracture geometry and natural fractures, the evolution of the reservoir stress tends to be complex in the production process. From a comparison of Figures 20(a) and Figure 18(a), the profile of the change in reorientation on Path B was found to be significantly different due to the influence of the hydraulic fracture geometry. During the first year of production, the central region did not turn; however, it turned in accordance with an increase in the production time. From a comparison of Figures 20 and Figure 18, there was a similarity in that the greatest change in the reorientation of σHmax occurred in the 10th year and the change in reorientation was correlated with the reservoir depletion time.Figure 19
(a) presents the direction of the initial maximum horizontal principal stress on thex-y section. (b) and (c) present the direction of the maximum horizontal principal stress after 3 years of depletion without considering the natural cracks. (d) and (e) present the direction of the maximum horizontal principal stress after 10 years of production considering the natural cracks.Figure 20
(a) and (b) present the change in reorientation ofσHmax along Path B without and with considering the influence of natural fractures and under different depletion years, respectively.
(a)
(b)


## 4.1. Basic Parameters

As discussed in this section, a three-dimensional reservoir physical model was established to evaluate the influence of anisotropy due to natural fractures on the reservoir stress during the oil and gas development. To fully demonstrate the significant influence of natural fractures on the generation of induced stress, a model was established that does not consider the influence ofkxz, kyz, and kxy. However, the influence of the permeability in the main direction was considered for comparison. As shown in Figure 3, the length, width, and height of the model were 400 m, 400 m, and 100 m, respectively, and the hydraulic fracture was located at the cross section of the y-z-axis. The domain was discretized to 50 four-node cells in the x- and y-directions and 20 four-node cells in the z-direction. For the treatment of the reservoir boundary, the nodes on the left and right surfaces of the model were free to move along the x-axis and not the y-axis. All the nodes at the top and bottom outer surfaces were free to move along the y-axis and not the x-axis. Only the vertical displacement of all the nodes on the outer surface of the top and bottom could be achieved. Moreover, the six boundaries were nonflow boundaries, i.e., the pore pressure in the reservoir was continuously reduced in the process of production. In all the calculations below, the downhole flow velocity was calculated using the Peaceman equation [37]. The other parameters in the calculation are shown in Table 1. Figures 4(a) and 4(b) present the locations of the single fractures and three clusters of nonplanar fractures, respectively.Figure 3
A schematic diagram of the 3D physical model.Figure 4
A schematic diagram of the hydrofracture layout ((a) presents the planar single fracture, and (b) presents the nonplanar multiple fractures).


## 4.2. Single Fracture Study

Based on the model and physical parameters established in Section4.1, in this section, the focus is on the influence of natural fractures on the stress distribution and steering in the production process of a single fracture in the y-z section. To fully reflect the effect of natural fractures on the stress, a conventional model was established without considering the effect of natural fractures, for comparison. Compared with the boundary conditions in the model, the calculated parameters were consistent with the proposed model.During the production, the pore pressure decreases continuously due to the continuous fluid extraction and the pressure is the lowest near the fracture. Figure5 compares the distribution of the reservoir pore pressure after 5 years of production with and without considering the influence of natural fractures. As can be seen from the figure, due to the influence of natural fractures, there was a change in the fluid flow pattern. Compared with the conventional model, the pore pressure dissipation area was greater and it was easier to reduce the pore pressure along the direction of natural fractures for the formation of a clockwise deflecting pressure drop funnel. Therefore, under the natural fracture distribution condition, the fluid in the upper right and lower left corners could be easily extracted, whereas the fluid in the upper left and lower right corners could not be easily extracted. This can serve as an essential guide to the design of infill wells.Figure 5
Panels (a) and (b), respectively, represent the pore pressure distribution without the influence of the natural fractures and the pore pressure under the influence of the natural fractures.
(a)
(b)According to the principle of effective stress, the reservoir stress changes under the influence of the constant dissipation of the pore pressure. For the conventional model, the distribution of the stressσxx, σyy, and σzz in the reservoir after 5 years of production is shown in Figures 6(a), 6(b), and 6(c). Due to fluid consumption, the reservoir stress distribution exhibited a high complexity. Moreover, there was a decrease in σxx near the crack as a result of the pressure dissipation, whereas σxx increased on the left and right parts of the domain to support the pressure depletion in the x-direction. In addition, there was a decrease in σyy near the crack as a result of the pressure dissipation, whereas σyy increased at the top and bottom parts of the domain to support the pressure depletion in the y-direction. Moreover, σzz mainly exhibited a decrease in stress in the central region of the domain. Considering the influence of natural fractures on the stress evolution, the results are shown in Figures 7(a), 7(b), and 7(c). Compared with the conventional model, the stress distribution of the proposed model is similar to that of the conventional model. However, considering the influence of tensors, the pressure dissipation of the proposed model was more rapid and the changes in stress were greater. The most significant difference is that the distributions of σxx, σyy, and σzz exhibited a degree of deflection. The stress distribution deflection was found to be significantly different from the stress distribution near the fractures in the conventional model, and this difference has a significant influence on the fracture initiation direction and fracture extension of refracturing. As can be seen in Figures 8(a) and 8(b), the distribution of the shear stress on the x-y section exhibited a clockwise deflection when compared with the conventional model, due to the influence of natural fractures. From a comparison between Figures 9(c) and 9(d) and Figures 9(a) and 9(b), the distribution of the displacement was found to be deflected to a certain extent under the influence of natural fractures. Figures 10(a) and 10(b) present the displacement profiles on Paths A and B, respectively. In all the figures, the solid lines represent the proposed model and the dotted lines represent the conventional model. The displacement trend of the two models was similar. The main difference was that there was a more significant change in the absolute value of the displacement after the influence of natural fractures was considered.Figure 6
(a), (b), and (c) present the principal stress in thex-direction, y-direction, and z-direction without considering the influence of the permeability tensor of natural fractures.
(a)
(b)
(c)Figure 7
(a), (b), and (c) present the principal stress in thex-direction, y-direction, and z-direction considering the influence of the permeability tensor of natural fractures.
(a)
(b)
(c)Figure 8
(a) and (b) present the induced shear stress in thex-y section with and without considering natural fractures, respectively.
(a)
(b)Figure 9
(a) and (b) presentux and uy on the x-y section without considering natural fractures, respectively, and (c) and (d) present ux and uy on the x-y section considering natural fractures, respectively.
(a)
(b)
(c)
(d)Figure 10
(a) and (b) present the displacement profiles on Paths A and B.
(a)
(b)Figures11(a), 11(b), and 11(c) present the σxx, σyy, and σzx stress distributions of the proposed model and conventional model after 5 years of production on the z-x section, respectively. As can be seen from the figure, the stress distribution of the conventional model exhibited a symmetrical distribution. Compared with the conventional model, the overall trend of the stress distribution in the proposed model was in good agreement. However, the stress distribution was symmetrically distributed under the influence of natural fractures, and there was a relative increase in the variation range of the stress value.Figure 11
(a), (b), and (c) present the comparison of the principal stress along thex-direction, y-direction, and the induced shear stress on the z-x section, respectively.
(a)
(b)
(c)Significant attention has been directed toward the reservoir stress evolution and surrounding drainage in the process of fluid recovery, given their significance for the latter refracturing and completion of infill wells. The change in the direction of the horizontal principal stress direction has a significant influence on the design of refracturing. Due to the heterogeneity of the permeability and natural fractures, the heterogeneity consumption of the pore pressure may result in stress reorientation and changes in the direction of the horizontal maximum principal stress. As can be seen from Figures12(a) and 12(b), a region with a stress difference less than 0 was observed near the fracture after one year of production. In addition, the stress in this region was reversed. It can also be found that the time when the stress difference changes the most is not the 20 years with the longest production life but the 10th year. This is because with an increase in the production life, most of the fluid in the easily flowing region was produced, whereas the fluid in the previously relatively noneasily flowing region was gradually produced. Therefore, there was a decrease in the difference between the pore pressures in the regions, which led to a decrease in the stress differences and an increase in the influence range. Figures 12(c) and 12(d) present the shear stress profiles on Paths A and B. Compared with the conventional model, the shear stress absolute value of the proposed model is influenced more significantly by natural cracks.Figure 12
(a) and (b) are present the stress difference profiles on Paths A and B of different production years, and (c) and (d) present the shear stress profiles on Paths A and B of different production years.
(a)
(b)
(c)
(d)Refracturing is one of the most effective methods for latter reservoir development. The distribution of the reservoir stress after the production and the direction of the horizontal maximum principal stress are two critical factors that determine the fracture pressure, fracture initiation direction, and fracture extension trajectory of refracturing. The stress evolution in the production process was previously analyzed. For the solution of the change in reorientation of the maximum horizontal principal stress (σHmax) on the x-y section, the action of the vertical principal stress (σzz) was neglected; thus, the original direction of the maximum horizontal principal stress is shown in Figure 13(a). After 5 and 10 years of production, there was a significant change in the direction of σHmax. In the conventional model, the change area of the principal stress direction was in the shape of a vertical ellipse, whereas the change area of the σHmax direction in the proposed model was in the shape of an inclined ellipse. The degree of inclination is related to the location of natural fractures. The area with the largest reorientation change of σHmax was near the hydrofracture. Compared with the conventional model, the change in reorientation of σHmax in the proposed model was influenced more significantly by natural fractures (θ1′>θ1, θ2′>θ2). In the conventional model, the change in reorientation of σHmax was 90° (Figure 14(a)); however, the change in reorientation of σHmax in the proposed model was 110° under the influence of natural cracks (Figure 14(b)). Similar to the change in the stress difference, the change in the reorientation was the largest at the 10th year, not the 20th year. Moreover, the range of the reorientation areas is proportional to production time.Figure 13
(a) presents the direction of the maximum principal stress on thex-y section in the initial state. (b) and (c) present the direction of the maximum horizontal principal stress after 5 years and 10 years of production without considering natural fractures. (d) and (e) present the direction of the maximum horizontal principal stress after 5 years and 10 years of production considering the influence of natural fractures.Figure 14
(a) and (b) present the steering angle of the maximum horizontal principal stress along Path B without and with considering the influence of natural fractures in different production years, respectively.
(a)
(b)


## 4.3. Nonplanar Multifracture Study

In the actual segmented multicluster fracturing process, there is stress interference between the fractures when multiple fractures are simultaneously extended, which leads to the case wherein the middle fractures are restrained and the fractures on both sides are extended in a nonplanar manner. Hence, the displacement discontinuity method was used to simulate the extended trajectories of three clusters of hydrofractures with spacings of 50 m. The three hydraulic fractures are considered as production fractures. It was assumed that all three clusters of fractures pass through the reservoir in thez-direction. The other calculation parameters are the same as those in Table 1. Under this condition, the reservoir stress distribution after 1, 5, and 10 years of production was simulated. Figure 15 presents the pore pressure distribution after 1, 5, and 10 years of nonplanar fracture production using the conventional model. The pore pressure distribution in the domain was axially symmetric. However, under the influence of natural fractures, there was an initial decrease in the pore pressure along the strike direction of the natural fractures, and the distribution exhibited central symmetry rather than axial symmetry (as shown in Figure 16). The upper left part and the lower right part of the domain are underexploited areas under such reservoir geological conditions and distributed fractures.Figure 15
(a), (b), and (c) present the pore pressure distribution after 1, 5, and 10 years of nonplanar fracture production without considering the influence of natural fractures.
(a)
(b)
(c)Figure 16
(a), (b), and (c) present the pore pressure distribution after 1, 5, and 10 years of nonplanar fracture production considering the influence of natural fractures.
(a)
(b)
(c)Due to the influence of the nonplanar fracture geometry, the stress distribution after production exhibited a high complexity. Without considering the effect of natural fractures, the stress distribution of the conventional model was axially symmetric. Compared with a single fracture, there was an increase in the rates of the pore pressure dissipation and stress decrease and there was an increase in the size of the influence range, as shown in Figures17(a), 17(b), and 17(c). The stress (σxx) on the left and right of the domain and the stress (σyy) on the top and bottom of the domain increased to support the pressure depletion of the middle area. When the influence of natural fractures was considered, the stress distribution in the proposed model exhibited a significant degree of deflection, as shown in Figures 18(a), 18(b), and 18(c). The deflection of the stress distribution confirms that natural fractures have a significant influence on the reservoir fluid flow and production.Figure 17
(a), (b), and (c) present theσx, σy, and σz stress distributions after 10 years of nonplanar fracture production without considering the influence of natural fractures.
(a)
(b)
(c)Figure 18
(a), (b), and (c) present theσx, σy, σz stress distribution after 10 years of nonplanar fracture production considering the influence of natural fractures.
(a)
(b)
(c)Compared with a single crack, the distribution of the principal stress was more complex. The change in reorientation ofσHmax was θ1′>θ1, θ2′>θ2, as shown in Figure 19. Due to the influence of the nonplanar crack geometry, the rotation angle of the maximum principal stress direction in the conventional model reached 135°. In the proposed model, the rotation angle of the principal stress direction near the crack was concentrated at 140°. The reorientation region affected by the natural fractures was larger than those of the conventional models. Therefore, under the action of the nonplanar fracture geometry and natural fractures, the evolution of the reservoir stress tends to be complex in the production process. From a comparison of Figures 20(a) and Figure 18(a), the profile of the change in reorientation on Path B was found to be significantly different due to the influence of the hydraulic fracture geometry. During the first year of production, the central region did not turn; however, it turned in accordance with an increase in the production time. From a comparison of Figures 20 and Figure 18, there was a similarity in that the greatest change in the reorientation of σHmax occurred in the 10th year and the change in reorientation was correlated with the reservoir depletion time.Figure 19
(a) presents the direction of the initial maximum horizontal principal stress on thex-y section. (b) and (c) present the direction of the maximum horizontal principal stress after 3 years of depletion without considering the natural cracks. (d) and (e) present the direction of the maximum horizontal principal stress after 10 years of production considering the natural cracks.Figure 20
(a) and (b) present the change in reorientation ofσHmax along Path B without and with considering the influence of natural fractures and under different depletion years, respectively.
(a)
(b)


## 5. Conclusion

A fully coupled 3D pore elastic model based on the FEM was developed to analyze the stress evolution and stress reorientation in heterogeneous porous media in the reservoir depletion process. In the model, the continuous natural fracture model is used to characterize the effect of natural fractures on the reservoir fluid flow. The relationship between the natural fracture width and effective stress was used to establish a fluid cluster single fracture model and a three-cluster nonplane fracture model for an example analysis. There was a good agreement between the analytical solutions and numerical results. By comparing the conventional models, it was confirmed that natural cracks have a significant influence on the stress evolution, surrounding drainage, and change in reorientation in the depletion process. Under the influence of natural fractures, the fluids in reservoirs along the direction of the natural fractures are easier to recover, and areas that are perpendicular to the natural fracture surface are prone to underconsumption. In addition, under the flow mode dominated by natural fractures, the stress distribution in the proposed model is deflected to a certain extent, especially in the nonplanar fracture model. Due to the action of natural cracks, the absolute values of the stress, displacement, and stress difference in the proposed model are relatively greater and the change in the reorientation of the maximum principal stress steering is significantly different. After considering the natural cracks, the steering angle and steering range increased. Moreover, this directly results in significant differences in the design of the late-stage fracturing temporary plugging steering and refracturing technology for conventional reservoirs. The simulation results revealed that it is very important to consider the effects of natural fractures on the stress redistribution, surrounding drainage, and reorientation in production. The findings of this study can serve as a theoretical guide for the design of refracturing, temporary plugging steering, and infilling wells and help predict the direction of new fractures and the trajectory of fracture extensions.



---

*Source: 1481960-2019-12-27.xml*