# A Novel Continuous Fracture Network Model: Formation Mechanism, Numerical Simulation, and Field Application

Authors: Zixi Guo; Jinzhou Zhao; Xinhao Sun; Chengwang Wang; Dali Guo; Haoran Hu; Hongna Wang; Qinggang Zeng
Journal: Geofluids
 (2022)

---

## Abstract

After a shale gas reservoir is fractured, hydraulic fractures interweave with natural fractures to form a fracture network. Numerical simulation based on the continuous fracture network model is a relatively economical and convenient method to predict fracture network morphology and size in the field application. However, some important factors, such as fracture height variation and filtration loss, have not been considered in the past continuous fracture network models. Therefore, this paper is aimed at establishing a novel continuous fracture network model to improve simulation accuracy. Firstly, this paper established a method to judge whether natural fractures develop or not. Then, a novel continuous fracture network model considering fracture height variation and asymmetry, filtration loss, fluid flow, and other key factors was established, and the forward algorithm and inverse algorithm of the model were proposed. At last, this model was applied in a field case to verify accuracy, and the average accuracy is more than 90%. Compared with the traditional Meyer software, the average error of prediction was reduced by 7.86%.


---

## Body

## 1. Introduction

Hydraulic fracturing has been used commonly in the development of shale gas reservoirs to initiate several dominant hydraulic fractures from the wellbore [1–4]. The technique is aimed at opening the widely existing natural fractures in the shale formation and generating the largest possible hydraulic fracture network for transporting hydrocarbon to the wellbore [5–9]. Thus, accurate knowledge of the fracture initiation and propagation of hydraulic fractures is significant for optimizing the fracturing treatment design and, ultimately, improving well productivity [10–12].Many key parameters affect the complexity of a fracture network, including the existence of natural fractures, anisotropy of stress, and heterogeneity of rock properties [13, 14]. Moreover, simulating hydraulic fracture propagation in a formation with preexisting natural fractures is very complex and requires proper consideration of key physical elements, including fracture propagation [15], fluid flow in the fracture networks [16], fracture height growth [17], interaction between hydraulic and natural fractures [18], and proppant transport in the fracture networks [19]. Microseismic and other geophysical methods are the main techniques for monitoring a fracture network in a shale reservoir [20]. However, the uncertainty of microseismic events and the high cost of fracture monitoring make the implementation of these techniques difficult for multiple wells [21].Numerical simulation is a relatively economical and convenient method to predict fracture morphology and size [22]. Different modeling approaches have been developed recently to simulate fracture networks, mainly divided into the discrete fracture network model [23, 24] and the continuous fracture network model. The main difference between these two models is that the discrete model fully considers the influence of natural fractures, while the continuous model uses the orthogonal fracture model to approximately describe the natural fracture network. In theory, the discrete model is closer to reality. However, in field applications, natural fracture diagnosis is limited by technology and high cost. Therefore, both methods have their own research value. This paper is aimed at developing and applying a new continuous fracture network model.The important breakthroughs in the development of fracture network models are summarized below: In 2008, Olson presented a continuous fracture network model capable of predicting hydraulic fracture propagation and interaction with preexisting natural fractures [25]. However, the model is based on fracture mechanics alone and does not include fluid flow or proppant transport. From 2009 to 2010, Xu et al. proposed a wire-mesh model to estimate fracture network dimensions and proppant placement in a network [26, 27]. This model has the advantage of fast computation speed due to being a semianalytical solution. But its fracture network pattern (i.e., fracture spacing) is relatively simple in simulating the geometric shape of fractures. In addition, the wire-mesh model neglects changes in fracture height and filtration of fracturing fluid. After 2011, Meyer established a relatively complete continuous fracture network model and developed commercial software [28], but its core problem was that the algorithm of model was not fully published, and this model did not take into account asymmetry of fracture height. In recent years, the discrete fracture network (DFN) model proposed by Meyer has been improved and extended by many scholars [29, 30]. For example, in 2017, Nejadi et al. proposed history matching and uncertainty quantification of DFN models in fractured reservoirs [31]. In 2020, Yao et al. discussed the role of natural fracture characteristics on the performance of hydraulic fracturing for deep energy extraction using DFN [32]. In 2021, Hyman and Dentz discussed the transport upscaling under flow heterogeneity and matrix-diffusion in three-dimensional DFN [33]. The development of the fracture network model is summarized in Table 1.Table 1
Comparison of different continuous fracture network models.
ReferenceTimeBaseInnovationOlson [25]2008OriginalContinuous fracture network modelXu et al. [26, 27]2009CFNSemianalytical solution modelMeyer et al. [28]2011CFNRelatively complete model and softwareNejadi et al. [31]2017DFNHistory matching and uncertainty quantificationYao et al. [32]2020DFNRole of natural fracture characteristicsHyman and Dentz [33]2021DFNTransport upscaling under flow heterogeneity and matrix-diffusionFrom the above literature review, a number of remaining challenges need to be resolved. To properly simulate the propagation of continuous fractures developed during hydraulic fracturing and get more accurate simulation results, a new model is needed. Therefore, this work researches the mechanical mechanism associated with fracture network formation based on fracture mechanics and proposes a kind of extended model of a continuous fracture network. This study establishes a mathematical model of the fracture network, including equations of width, fracture length, and fracture height. The model solves a system of equations governing fracture deformation, height growth, fluid flow, and proppant transport so as to predict the size of a complicated fracture network more accurately.The paper is organized as follows. Firstly, it provides mechanical condition requirements related to fracture network formation to determine whether natural fractures have developed or not. Secondly, it establishes a continuous fracture network extension model according to the propagation law of continuous fracture networking and considers the change in fracturing fluid filtration and fracture height. Finally, it uses field experimental data to verify the accuracy of the model.


## 2. Mechanical Mechanism of Fracture Network Formation Based on Fracture Mechanics

Warpinski et al. [28, 34, 35] pointed out that fracture networks are the result of interaction between hydraulic fractures and natural fractures; researchers have held that a continuous fracture network consists of one main fracture and several branch fractures, as shown in Figure 1. The key to fracture network formation lies in the fact that branch fractures develop around the main fracture [36, 37]. Some natural fractures develop, while some do not. This work studied the mechanical mechanism of reservoir formation with and without natural fracture development.Figure 1
Schematic diagram of the formation of a fracture network.


### 2.1. Mechanical Condition Requirements for Reservoir Formation with Natural Fracture Development

For reservoirs with natural fracture development, when the net pressure in the main fracture of a certain length increases after its formation, a natural fracture or weak plane opens, and then, a fracture network is formed. At present, the linear criterion put forward by Warpinski and Teufel [38, 39] is the most widely used fracture propagation criterion. The mechanical condition requirements for the formation of branch fractures in a fracture network can be analyzed based on the fracture propagation of a naturally fractured reservoir. Figure 2 shows a diagram of a reservoir fracture network with natural fracture development.Figure 2
Schematic diagram of a naturally fractured reservoir fracture network.The closed natural fracture (or weak plane) is acted on by the maximum horizontal principal stress (σH) and minimum horizontal principal stress (σh) of the far-field, and its angle from the maximum horizontal principal stress is θ0<θ<π/2.According to the 2D linear elastic theory, shear stress and normal stress can be expressed as follows [40]:
(1)τ=−σH−σh2sin2θ,(2)σn=σH+σh2−σH−σh2cos2θ.According to the Warpinski criterion widely applied to fracture propagation of naturally fractured reservoirs, when the fluid pressure of a natural fracture is larger than its normal stress, fracture extension will occur [41] as follows:
(3)pt>σn+△pnft,(4)△pnft=4pt−p0π∑n=0∞12n+1∗exp−2n+12π2knft4ϕnfμCtLf2sin2n+1π2.When two fractures intersect, the hydraulic fracture end is connected with the natural fracture, and fracturing fluid floods into the natural fracture. The pore pressure of the natural fracture’s near wall is as follows:(5)pt=σh+pnett.By substituting Equations (2) and (5) into Equation (3), the net pressure required for the tension fracture of the natural fracture is as follows:
(6)pnett>σH−σh21−cos2θ+△pnft.Whenθ=π/2, the net pressure reaches the maximum:
(7)pmax=σH−σh+△pnft.According to the Warpinski criterion, when the pressure in the natural fracture is low, normal stress acting on the fracture surface is pressure stress, and the fracture is closed. When the shear stress acting on the natural fracture is high, the shear slip of the natural fracture will occur easily [42]:
(8)τ>τ0+Kfσn−pt.Effective shear stress (τe) acting on fracture surface is as follows:
(9)τe=τ−τ0−Kfσn−pt.The fracture criterion of a compression shear fracture is as follows:(10)KΠ=KΠc,(11)KΠ=τeπl.By substituting Equations (1), (2), and (11) into Equation (9), we get critical pressure (pc) for the occurrence of shear slip and branch fracture formation:
(12)pc=1KfKΠcπl+τ0+σH+σh2−σH−σh2cos2θ+sin2θKf.According to Equation (12), when θ=π/2, the net pressure is the largest, and pmax is as follows:
(13)pmax=△σ+τ0Kf.The natural fracture surface is unbounded (τ0=0). Net pressure required for shear slip is horizontal principal stress difference (Δσ).According to Equations (7) and (13), maximum net pressure is positively correlated with . Therefore, the reduction in initial horizontal stress difference can help expand the fracture network fully.The calculation formula for the difference between two horizontal principal stresses of a reservoir can be derived based on the maximum horizontal principal stress formula [(36)]:
(14)△σ=2σh−Pi−Pf+St.Therefore, according toΔσ, combined with the net pressure in the fracture propagation model, a judgment can be made on whether a fracture network can be formed.


### 2.2. Mechanical Condition Requirements for Reservoir Formation without Natural Fracture Development

For a reservoir formation without natural fracture development, branch fractures in the rock body are required for fracture network formation [43]. According to the mechanical model shown in Figure 3, there is an elliptical fracture in the infinite field with the major semiaxis Lf and the minor semiaxis w. At infinity, the major semiaxis is acted on by pressure stress (σH) and the minor semiaxis (σh); there is uniform pressure action in the fracture.Figure 3
Schematic diagram of a naturally fractured reservoir fracture network.Boundary conditions:Aty=0,x<lf,
(15)σy=−pnet,τxy=0.Atx2+y2⟶∞,
(16)σx⟶σH,σy⟶σh,τxy⟶0.Expressions ofσx, σy, and τxy can be obtained according to the solutions for plane problems of elastic mechanics. Stress distribution on the fracture surface is expressed as follows:
(17)σθ=−1−3m2+2mcos2θ1+m2−2mcos2θpnet+1−m2−2m+2cos2θ1+m2−2mcos2θσh+1−m2+2m+2cos2θ1+m2−2mcos2θσH,σp=−pnet,τpθ=0,where m=lf−w/lf+w.Becauself≫w, m≈1, and putting into Equation (8), we obtain the following:
(18)σθ=pnet−σh+σH.According to the elastic failure criterion, settingσθ=−St, we obtain the following:
(19)pnet=−σH−σh−St.Therefore, when the value of net pressure in a fracture is larger thanσH−σh+St, the fracture in the rock body breaks and forms branch fractures, and then, a fracture network is formed.To sum up, mechanical condition requirements for a fracture network in a naturally fractured reservoir lie in net pressure in the constructed fracture exceeding the horizontal principal stress difference of the reservoir. For a fracture network in a reservoir without natural fracture development, mechanical condition requirements lie in the fracture breaking in the rock body and the value of net pressure in the fracture being larger than the sum of horizontal principal stress difference and reservoir tensile strength.


## 2.1. Mechanical Condition Requirements for Reservoir Formation with Natural Fracture Development

For reservoirs with natural fracture development, when the net pressure in the main fracture of a certain length increases after its formation, a natural fracture or weak plane opens, and then, a fracture network is formed. At present, the linear criterion put forward by Warpinski and Teufel [38, 39] is the most widely used fracture propagation criterion. The mechanical condition requirements for the formation of branch fractures in a fracture network can be analyzed based on the fracture propagation of a naturally fractured reservoir. Figure 2 shows a diagram of a reservoir fracture network with natural fracture development.Figure 2
Schematic diagram of a naturally fractured reservoir fracture network.The closed natural fracture (or weak plane) is acted on by the maximum horizontal principal stress (σH) and minimum horizontal principal stress (σh) of the far-field, and its angle from the maximum horizontal principal stress is θ0<θ<π/2.According to the 2D linear elastic theory, shear stress and normal stress can be expressed as follows [40]:
(1)τ=−σH−σh2sin2θ,(2)σn=σH+σh2−σH−σh2cos2θ.According to the Warpinski criterion widely applied to fracture propagation of naturally fractured reservoirs, when the fluid pressure of a natural fracture is larger than its normal stress, fracture extension will occur [41] as follows:
(3)pt>σn+△pnft,(4)△pnft=4pt−p0π∑n=0∞12n+1∗exp−2n+12π2knft4ϕnfμCtLf2sin2n+1π2.When two fractures intersect, the hydraulic fracture end is connected with the natural fracture, and fracturing fluid floods into the natural fracture. The pore pressure of the natural fracture’s near wall is as follows:(5)pt=σh+pnett.By substituting Equations (2) and (5) into Equation (3), the net pressure required for the tension fracture of the natural fracture is as follows:
(6)pnett>σH−σh21−cos2θ+△pnft.Whenθ=π/2, the net pressure reaches the maximum:
(7)pmax=σH−σh+△pnft.According to the Warpinski criterion, when the pressure in the natural fracture is low, normal stress acting on the fracture surface is pressure stress, and the fracture is closed. When the shear stress acting on the natural fracture is high, the shear slip of the natural fracture will occur easily [42]:
(8)τ>τ0+Kfσn−pt.Effective shear stress (τe) acting on fracture surface is as follows:
(9)τe=τ−τ0−Kfσn−pt.The fracture criterion of a compression shear fracture is as follows:(10)KΠ=KΠc,(11)KΠ=τeπl.By substituting Equations (1), (2), and (11) into Equation (9), we get critical pressure (pc) for the occurrence of shear slip and branch fracture formation:
(12)pc=1KfKΠcπl+τ0+σH+σh2−σH−σh2cos2θ+sin2θKf.According to Equation (12), when θ=π/2, the net pressure is the largest, and pmax is as follows:
(13)pmax=△σ+τ0Kf.The natural fracture surface is unbounded (τ0=0). Net pressure required for shear slip is horizontal principal stress difference (Δσ).According to Equations (7) and (13), maximum net pressure is positively correlated with . Therefore, the reduction in initial horizontal stress difference can help expand the fracture network fully.The calculation formula for the difference between two horizontal principal stresses of a reservoir can be derived based on the maximum horizontal principal stress formula [(36)]:
(14)△σ=2σh−Pi−Pf+St.Therefore, according toΔσ, combined with the net pressure in the fracture propagation model, a judgment can be made on whether a fracture network can be formed.


## 2.2. Mechanical Condition Requirements for Reservoir Formation without Natural Fracture Development

For a reservoir formation without natural fracture development, branch fractures in the rock body are required for fracture network formation [43]. According to the mechanical model shown in Figure 3, there is an elliptical fracture in the infinite field with the major semiaxis Lf and the minor semiaxis w. At infinity, the major semiaxis is acted on by pressure stress (σH) and the minor semiaxis (σh); there is uniform pressure action in the fracture.Figure 3
Schematic diagram of a naturally fractured reservoir fracture network.Boundary conditions:Aty=0,x<lf,
(15)σy=−pnet,τxy=0.Atx2+y2⟶∞,
(16)σx⟶σH,σy⟶σh,τxy⟶0.Expressions ofσx, σy, and τxy can be obtained according to the solutions for plane problems of elastic mechanics. Stress distribution on the fracture surface is expressed as follows:
(17)σθ=−1−3m2+2mcos2θ1+m2−2mcos2θpnet+1−m2−2m+2cos2θ1+m2−2mcos2θσh+1−m2+2m+2cos2θ1+m2−2mcos2θσH,σp=−pnet,τpθ=0,where m=lf−w/lf+w.Becauself≫w, m≈1, and putting into Equation (8), we obtain the following:
(18)σθ=pnet−σh+σH.According to the elastic failure criterion, settingσθ=−St, we obtain the following:
(19)pnet=−σH−σh−St.Therefore, when the value of net pressure in a fracture is larger thanσH−σh+St, the fracture in the rock body breaks and forms branch fractures, and then, a fracture network is formed.To sum up, mechanical condition requirements for a fracture network in a naturally fractured reservoir lie in net pressure in the constructed fracture exceeding the horizontal principal stress difference of the reservoir. For a fracture network in a reservoir without natural fracture development, mechanical condition requirements lie in the fracture breaking in the rock body and the value of net pressure in the fracture being larger than the sum of horizontal principal stress difference and reservoir tensile strength.


## 3. Continuous Fracture Network Modeling

The physical model of network fractures is shown in Figure4, which is composed of two clusters of orthogonal equally spaced vertical fractures. They are elliptical in the transverse direction, but the fracture height is not constant in the longitudinal direction, so they can be asymmetric in the upper and lower directions. The advantage of this physical model is that when the distribution of natural fractures is not clear, natural fractures and other factors are contained in the model, and it is a reasonable approximation of fractures in complex networks.Figure 4
Schematic diagram of a continuous fracture network.


### 3.1. Building Fracture Network Modeling

In order to establish the mathematical model of the fracture network, the coordinate system was established. The origin of the coordinates was taken as the shaft injection point, thex axis was the direction of the maximum horizontal principal stress, and the y axis was the direction of the minimum horizontal principal stress. Suppose that the half-length of elliptical network fracture on x axis and y axis is Lx and Ly, respectively, and the fracture spacing on x axis and y axis is dx and dy, respectively, then the number of fractures parallel to x axis and y axis is 2ny+1 and 2nx+1, where
(20)nm=intLmdmm=x,y.According to the elliptic equation, the half-length of thejth fracture parallel to the x axis is
(21)Lxj=Lx1−jdyLy2j=0,±1,⋯,±ny.The half-length of theith fracture parallel to the y axis is
(22)Lyi=Ly1−idxLx2i=0,±1,⋯,±nx.In the process of network fracture propagation, the fluid flows from the center of the ellipse to the edge of the ellipse and presents laminar flow in the fracture. The mathematical model describing network fracture propagation is established by considering fracture fluid filtration and fracture height variation.


#### 3.1.1. Volumetric Equilibrium Equation

According to the principle of mass conservation and ignoring fluid compressibility and initial filtration loss of fracturing fluid, the volume of fracturing injection is equal to the sum of the volume of network fractures and the filtration loss volume of fracturing fluid, that is(23)∫0tqτdτ=Vfxt+Vfyt−Vlxt+Vlyt,where Vfx and Vfy are the volume of fractures parallel to x axis and y axis, respectively. Vlx and Vly are the filtration volume of fracturing fluid for fractures parallel to x axis and y axis, respectively.According to the morphology and composition of the fracture network, the cross-section of fractures parallel to thex axis is approximately elliptical, and the volume of fractures is
(24)Vfx=∑j=−nyny∫−LxjLxjπ4wxjhxjdx,where wxj is the fracture width of the jth fracture parallel to the x axis. hxj is the fracture height of the jth fracture parallel to the x axis.Similarly, the volume of the fracture parallel to they axis is
(25)Vfy=∑i=−nxnx∫−LyiLyiπ4wyihyidy,where wyi is the fracture width of the ith fracture parallel to the y axis. hyi is the fracture height of the ith fracture parallel to the y axis.The volume of the fracture is equal to the sum of the volume of fractures parallel to thex axis and the volume of fractures parallel to the y axis, that is:
(26)Vf=∑j=−nyny∫−LxjLxjπ4wxjhxjdx+∑i=−nxnx∫−LyiLyiπ4wyihyidy.According to the morphology and composition of the fracture network, the fracturing fluid filtration volume of fractures parallel to thex axis can also be obtained by the following equation:
(27)Vlx=∑j=−nyny∫−LxjLxj∫0t2hxjCLt−τ0sdsdx.The fracturing fluid filtration volume of the fracture parallel to they axis is
(28)Vly=∑i=−nxnx∫−LyiLyi∫0t2hyiCLt−τ0sdsdy.


#### 3.1.2. Fracture Width Equation

In thejth fracture parallel to the x axis, considering the variation of crustal stress in the longitudinal direction, when the fracture penetrates the upper and lower layers (hxj≥Hp), the net pressure distribution in the fracture is
(29)Δpxjz=pxj−σhuaxj≤z≤lxjpxj−σhbxj≤z≤axjpxj−σhl−lxj≤z≤bxjj=0,±1,⋯,±ny,where
(30)lxj=hxj2=hxju+hxjl2,axj=Hp+hxjl−hxju2,bxj=−Hp−hxjl+hxju2.According to the elastic deformation theory, the net pressure in the fracture is decomposed into the sum of even stress distribution and odd stress distribution, and then, the fracture width profile is calculated by the England and Green formulas. After deduction, the fracture width is(31)wxj=41−v2lxjEpxj−σh−41−v2lxjπEσhu−σhcos−1axjlxj−axjlxjlnlxj+lxj2−axj2axj−41−v2lxjπEσhl−σhcos−1bxjlxj−bxjlxjlnlxj+lxj2−bxj2bxjj=0,±1,⋯,±ny.When the fracture does not penetrate the overlying and underlying strata (hxj<Hp) or the longitudinal variation of crustal stress is ignored, the fracture width of the jth fracture parallel to the x axis is
(32)wxj=21−v2hxjEpxj−σhj=0,±1,⋯,±ny.Similarly, when the fracture penetrates the upper and lower layers (hyi≥Hp) or the longitudinal variation of crustal stress is ignored, the fracture width of the ith fracture parallel to the y axis is
(33)wyi=41−v2lyiEpyi−σH−41−v2lyiπEσHu−σHcos−1ayilyi−ayilyilnlyi+lyi2−ayi2ayi−41−v2lyiπEσHl−σHcos−1byilyi−byilyilnlyi+lyi2−byi2byii=0,±1,⋯,±nx,where
(34)lyi=hyi2=hyiu+hyil2,ayi=Hp+hyil−hyiu2,byi=−Hp−hyil+hyiu2.When the fracture does not penetrate the upper and lower layers (hyi<Hp) or the longitudinal variation of crustal stress is ignored, the fracture width of the ith fracture parallel to the y axis is
(35)wyi=21−v2hyiEpyi−σHi=0,±1,⋯,±nx.


#### 3.1.3. Fracture Height Equation

According to the theory of linear elastic fracture mechanics, when the fracture penetrates the upper and lower layers (hxj≥Hp), the stress intensity factor at the upper and lower ends of the jth fracture parallel to the x axis is
(36)KIu=1πlxj∫−lxjlxjΔpxjlxj+zlxj−zdz,(37)KIl=−1πlxj∫lxj−lxjΔpxjlxj−zlxj+zdz.According to the net pressure distribution in the fracture and letKIu=KICu and KIl=KICl, the fracture height equation of the jth fracture parallel to the x axis is
(38)πlxjKICu=πlxjpxj−σhu+σhl2+σhl−σhlxj2−bxj2−lxjsin−1bxjlxj−σhu−σhlxj2−axj2−lxjsin−1axjlxj,πlxjKICl=πlxjpxj−σhu+σhl2−σhl−σhlxj2−bxj2+lxjsin−1bxjlxj+σhu−σhlxj2−axj2+lxjsin−1axjlxj.When the fracture does not penetrate the upper and lower layers (hxj<Hp) or the longitudinal variation of crustal stress is ignored, the fracture height equation of the jth fracture parallel to the x axis is
(39)KIC=pxj−σhπlxjj=0,±1,⋯,±ny.Similarly, when the fracture penetrates the upper and lower layers (hyi≥Hp), the fracture height equation of the ith fracture parallel to the y axis is
(40)πlyiKICu=πlyipyi−σHu+σHl2+σHl−σHlyi2−byi2−lyisin−1byilyi−σHu−σHlyi2−ayi2−lyisin−1ayilyi,πlyiKICl=πlyipyi−σHu+σHl2−σHl−σHlyi2−byi2+lyisin−1byilyi+σHu−σHlyi2−ayi2+lyisin−1ayilyi.When the fracture does not penetrate the upper and lower layers (hyi<Hp) or the longitudinal variation of crustal stress is ignored, the fracture height equation of the ith fracture parallel to the y axis is
(41)KIC=pyi−σHπlyii=0,±1,⋯,±nx.


#### 3.1.4. Fluid Flow Equation

According to the mechanics through porous media, the fluid flow equation in thejth fracture parallel to x axis is
(42)πλ∂ϕ∂t=1x∂∂xη1+ξxwxjkfxjμdx∂pxj∂xj=0,±1,⋯,±ny,where ξ is the axial ratio of the x axis to the y axis, ξ=Lx/Ly. η is the complete elliptic integrals of the second kind, η=∫0π/21−1−ξ2sin2sds. kfxj is the permeability of the jth fracture parallel to the x axis.Considering the flow in the fracture is in laminar, the permeability of thejth fracture parallel to the x axis is
(43)kfxj=π2hxj2pxj−σh48E2.Substitute into the fluid flow equation, it can be obtained:(44)∂ϕ∂t=π21−v2η1+ξ24λE3μxdx∂∂xxhxj3pxj−σh3∂pxj∂x.Considering the change of porosity with pressure, it is represented by rock compression coefficient, that is(45)CR=1ϕdϕdp.Integrate the above equation and the porosity can be obtained by:(46)ϕ=ϕ0eCRp−p0,where p0 is the initial pressure.Expand the above equation by series, omit the infinitesimal of higher order, and get(47)ϕ=ϕ0CRp−p0,so get
(48)∂ϕ∂t=ϕ0CR∂p∂t.Thus, in thejth fracture parallel to the x axis, the fluid flow equation is
(49)∂pxj∂t=Axx∂∂xBxpxj,x∂pxj∂xj=0,±1,⋯,±ny,where
(50)Axx=π21−v2η1+ξ24ϕ0CRλE3μxdx,Bxpxj,x=xhxj3pxj−σh3.The initial conditions and boundary conditions are(51)pxjt=0=σhpxjx=Lxj=σh,∂pxj∂xx=Lxj=0.Similarly, in theith fracture parallel to the y axis, the fluid flow equation is
(52)∂pyi∂t=Ayy∂∂yBypyi,y∂pyi∂yi=0,±1,⋯,±nx,where
(53)Ayy=π21−v2η1+ξ24ϕ0CRλE3μydy,Bypyi,y=yhyi3pyi−σH3.The initial conditions and boundary conditions are(54)pyit=0=σHpyiy=Lyi=σH,∂pyi∂yy=Lyi=0.


#### 3.1.5. Coupling Conditions

Each fracture of the network (whether parallel tox axis or y axis) is relatively independent and interrelated, constituting a complex network fracture system. Each fracture is described by its own fluid flow equation, fracture width equation, and fracture height equation and is controlled by volume balance equation. At the same time, there is coupling condition; that is, the pressure at the origin of coordinates is equal, expressed as
(55)px0x=0=py0y=0.


### 3.2. Numerical Solution of the Fracture Network Model

#### 3.2.1. Difference Method for Fluid Flow Equations

It is a nonlinear partial differential equation for the fluid flow equation of thejth fracture parallel to the x axis. We use the implicit scheme difference method to solve it, and the nonlinear problem is solved by the coefficient freezing method. Thus, the left end of the fluid flow equation is discretized as follows:
(56)∂pxj∂t≈pxjkn+1−pxjknΔt.The right end of the fluid flow equation is discretized as follows:(57)Ax∂∂xBx∂pxj∂x≈AxknBxk+1/2npxjk+1n+1−pxjkn+1/Δx−Bxk−1/2npxjkn+1−pxjk−1n+1/ΔxΔx,where Δt is the time step and Δx is the length step in the x direction.
(58)Bxk+1/2n≈Bxk+1n+Bxkn2,Bxk−1/2n≈Bxkn+Bxk−1n2.Thus, the difference equation of the fluid flow equation is(59)−ΔtAxknBxk+1/2nΔx2pxjk+1n+1+1+ΔtAxknBxk+1/2n+Bxk−1/2nΔx2pxjkn+1−ΔtAxknBxk−1/2nΔx2pxjk−1n+1=pxjkn.These are tridiagonal linear equations with dominant diagonal, which can be solved by the catch-up method combined with their initial conditions and boundary conditions.Similarly, for theith fracture parallel to the y axis, the difference equation of the fluid flow equation (3-108) can be expressed as follows:
(60)−ΔtAyknByk+1/2nΔy2pyik+1n+1+1+ΔtAyknByk+1/2n+Byk−1/2nΔy2pyikn+1−ΔtAyknByk−1/2nΔy2pyik−1n+1=pyikn,where Δy is the length step in the y direction.


#### 3.2.2. Forward Algorithm for the Numerical Solution

The main purpose of establishing and solving the fracture network model is known basic parameters (such as elastic modulusE, Poisson’s ratio ν and other rock parameters, fracturing fluid viscosity μ, injection rate q, injection time T, and other construction parameters); by solving the model, the fracture network parameters can be obtained (such as half-length Lx and Ly, fracture width, fracture height, and Vf).It is not hard to see that each fracture of the network is described by its own fluid flow equation, fracture height equation, and fracture width equation, which is a complete mathematical model. In other words, under the condition of known basic parameters, fracture spacing (dx and dy) and fracture half-length (Lx and Ly), the fracture length, fracture width, and fracture height and pressure in each fracture can be calculated by using this model. Therefore, by using the volume balance equation and coupling conditions (the pressure is equal at the intersection), the algorithm is established to obtain some parameters of the fracture spacing and the fracture half-length.The forward algorithm in this paper is to obtain the half-length (Lx and Ly) of fracture network, the length, width, height, pressure, and Vf of each fracture under the condition of known basic parameters and fracture spacing (dx and dy).The specific step of the forward algorithm is as follows:(1)
The coupling conditions are transformed to construct the following function:(61)Jd1Lx,Ly=px0x=0−py0y=0.As you can see that for anyLx and , Jd1 can be calculated by the function. Theoretical analysis and practical calculation show that if Jd1=0, then Lx =Ly. Therefore, for a fixed value Lx, according to the condition Jd1=0, it is easy to calculate the value of Ly by iterative method. At this time, the iteration formula is designed as follows:
(62)Ly=1±λLy,where λ is the iteration factor which is set to 0.1 during trial calculation.
(2)
According to the volume equilibrium equation, the following function is constructed:(63)Jd2Lx=∫0tqτdτ−Vfxt+Vfyt+Vlxt+Vlyt.It is not difficult to see that for any value ofLx, since the value of Ly is determined by the condition Jd1=0, Jd2 is just a function of Lx. Theoretical analysis and practical calculation show that if Jd2=0, then the value of Lx reaches volumetric balance. Therefore, according to the condition Jd2=0, it is easy to calculate the value of Lx by iterative method. At this time, the iteration formula is designed as follows:
(64)Lx=1±λLx.(3)
The forward algorithm advances step by step by time stepΔt. Because the implicit difference method is used to solve the fluid flow equation, the numerical calculation is unconditionally stable. However, the smaller the time step, the higher the accuracy of the calculation resultTo sum up, the forward algorithm designed is shown in Figure5.Figure 5
Flow chart of the forward algorithm for numerical solution.


#### 3.2.3. Inverse Algorithm for the Numerical Solution

The forward algorithm for the numerical solution of the network fracture model requires the fracture spacing (dx and dy) to be input as known parameters. Usually, the forward algorithm is used to check and invert the microseismic data of nearby wells, borrowing the fracture spacing of adjacent wells.In fact, parameter inversion is the inverse problem of parameter forward modeling, and the inversion algorithm is different with different parameters. The inversion algorithm in this paper is to obtain the fracture spacing (dx and dy), the fracture length, the fracture width, the fracture height, the pressure in the fracture, and Vf of each fracture through the algorithm under the condition that the basic parameters and the half-length (Lx and Ly) of the fracture network are known.The inversion algorithm of the fracture network numerical solution is similar to the forward algorithm, which is briefly described here. Take the coupling conditions of equal pressure at the origin of coordinates and construct the following function:(65)JL1dx,dy=px0x=0−py0y=0.The following iterative formula is adopted:(66)dy=1±λdy.At the same time, according to the volume equilibrium equation, the following function is constructed:(67)JL2dx=∫0tqτdτ−Vfxt+Vfyt+Vlxt+Vlyt.The following iterative formula is adopted:(68)dx=1±λdx.Therefore, the process of the inversion algorithm is shown in Figure6.Figure 6
Flow chart of the inverse algorithm for numerical solution.


## 3.1. Building Fracture Network Modeling

In order to establish the mathematical model of the fracture network, the coordinate system was established. The origin of the coordinates was taken as the shaft injection point, thex axis was the direction of the maximum horizontal principal stress, and the y axis was the direction of the minimum horizontal principal stress. Suppose that the half-length of elliptical network fracture on x axis and y axis is Lx and Ly, respectively, and the fracture spacing on x axis and y axis is dx and dy, respectively, then the number of fractures parallel to x axis and y axis is 2ny+1 and 2nx+1, where
(20)nm=intLmdmm=x,y.According to the elliptic equation, the half-length of thejth fracture parallel to the x axis is
(21)Lxj=Lx1−jdyLy2j=0,±1,⋯,±ny.The half-length of theith fracture parallel to the y axis is
(22)Lyi=Ly1−idxLx2i=0,±1,⋯,±nx.In the process of network fracture propagation, the fluid flows from the center of the ellipse to the edge of the ellipse and presents laminar flow in the fracture. The mathematical model describing network fracture propagation is established by considering fracture fluid filtration and fracture height variation.


### 3.1.1. Volumetric Equilibrium Equation

According to the principle of mass conservation and ignoring fluid compressibility and initial filtration loss of fracturing fluid, the volume of fracturing injection is equal to the sum of the volume of network fractures and the filtration loss volume of fracturing fluid, that is(23)∫0tqτdτ=Vfxt+Vfyt−Vlxt+Vlyt,where Vfx and Vfy are the volume of fractures parallel to x axis and y axis, respectively. Vlx and Vly are the filtration volume of fracturing fluid for fractures parallel to x axis and y axis, respectively.According to the morphology and composition of the fracture network, the cross-section of fractures parallel to thex axis is approximately elliptical, and the volume of fractures is
(24)Vfx=∑j=−nyny∫−LxjLxjπ4wxjhxjdx,where wxj is the fracture width of the jth fracture parallel to the x axis. hxj is the fracture height of the jth fracture parallel to the x axis.Similarly, the volume of the fracture parallel to they axis is
(25)Vfy=∑i=−nxnx∫−LyiLyiπ4wyihyidy,where wyi is the fracture width of the ith fracture parallel to the y axis. hyi is the fracture height of the ith fracture parallel to the y axis.The volume of the fracture is equal to the sum of the volume of fractures parallel to thex axis and the volume of fractures parallel to the y axis, that is:
(26)Vf=∑j=−nyny∫−LxjLxjπ4wxjhxjdx+∑i=−nxnx∫−LyiLyiπ4wyihyidy.According to the morphology and composition of the fracture network, the fracturing fluid filtration volume of fractures parallel to thex axis can also be obtained by the following equation:
(27)Vlx=∑j=−nyny∫−LxjLxj∫0t2hxjCLt−τ0sdsdx.The fracturing fluid filtration volume of the fracture parallel to they axis is
(28)Vly=∑i=−nxnx∫−LyiLyi∫0t2hyiCLt−τ0sdsdy.


### 3.1.2. Fracture Width Equation

In thejth fracture parallel to the x axis, considering the variation of crustal stress in the longitudinal direction, when the fracture penetrates the upper and lower layers (hxj≥Hp), the net pressure distribution in the fracture is
(29)Δpxjz=pxj−σhuaxj≤z≤lxjpxj−σhbxj≤z≤axjpxj−σhl−lxj≤z≤bxjj=0,±1,⋯,±ny,where
(30)lxj=hxj2=hxju+hxjl2,axj=Hp+hxjl−hxju2,bxj=−Hp−hxjl+hxju2.According to the elastic deformation theory, the net pressure in the fracture is decomposed into the sum of even stress distribution and odd stress distribution, and then, the fracture width profile is calculated by the England and Green formulas. After deduction, the fracture width is(31)wxj=41−v2lxjEpxj−σh−41−v2lxjπEσhu−σhcos−1axjlxj−axjlxjlnlxj+lxj2−axj2axj−41−v2lxjπEσhl−σhcos−1bxjlxj−bxjlxjlnlxj+lxj2−bxj2bxjj=0,±1,⋯,±ny.When the fracture does not penetrate the overlying and underlying strata (hxj<Hp) or the longitudinal variation of crustal stress is ignored, the fracture width of the jth fracture parallel to the x axis is
(32)wxj=21−v2hxjEpxj−σhj=0,±1,⋯,±ny.Similarly, when the fracture penetrates the upper and lower layers (hyi≥Hp) or the longitudinal variation of crustal stress is ignored, the fracture width of the ith fracture parallel to the y axis is
(33)wyi=41−v2lyiEpyi−σH−41−v2lyiπEσHu−σHcos−1ayilyi−ayilyilnlyi+lyi2−ayi2ayi−41−v2lyiπEσHl−σHcos−1byilyi−byilyilnlyi+lyi2−byi2byii=0,±1,⋯,±nx,where
(34)lyi=hyi2=hyiu+hyil2,ayi=Hp+hyil−hyiu2,byi=−Hp−hyil+hyiu2.When the fracture does not penetrate the upper and lower layers (hyi<Hp) or the longitudinal variation of crustal stress is ignored, the fracture width of the ith fracture parallel to the y axis is
(35)wyi=21−v2hyiEpyi−σHi=0,±1,⋯,±nx.


### 3.1.3. Fracture Height Equation

According to the theory of linear elastic fracture mechanics, when the fracture penetrates the upper and lower layers (hxj≥Hp), the stress intensity factor at the upper and lower ends of the jth fracture parallel to the x axis is
(36)KIu=1πlxj∫−lxjlxjΔpxjlxj+zlxj−zdz,(37)KIl=−1πlxj∫lxj−lxjΔpxjlxj−zlxj+zdz.According to the net pressure distribution in the fracture and letKIu=KICu and KIl=KICl, the fracture height equation of the jth fracture parallel to the x axis is
(38)πlxjKICu=πlxjpxj−σhu+σhl2+σhl−σhlxj2−bxj2−lxjsin−1bxjlxj−σhu−σhlxj2−axj2−lxjsin−1axjlxj,πlxjKICl=πlxjpxj−σhu+σhl2−σhl−σhlxj2−bxj2+lxjsin−1bxjlxj+σhu−σhlxj2−axj2+lxjsin−1axjlxj.When the fracture does not penetrate the upper and lower layers (hxj<Hp) or the longitudinal variation of crustal stress is ignored, the fracture height equation of the jth fracture parallel to the x axis is
(39)KIC=pxj−σhπlxjj=0,±1,⋯,±ny.Similarly, when the fracture penetrates the upper and lower layers (hyi≥Hp), the fracture height equation of the ith fracture parallel to the y axis is
(40)πlyiKICu=πlyipyi−σHu+σHl2+σHl−σHlyi2−byi2−lyisin−1byilyi−σHu−σHlyi2−ayi2−lyisin−1ayilyi,πlyiKICl=πlyipyi−σHu+σHl2−σHl−σHlyi2−byi2+lyisin−1byilyi+σHu−σHlyi2−ayi2+lyisin−1ayilyi.When the fracture does not penetrate the upper and lower layers (hyi<Hp) or the longitudinal variation of crustal stress is ignored, the fracture height equation of the ith fracture parallel to the y axis is
(41)KIC=pyi−σHπlyii=0,±1,⋯,±nx.


### 3.1.4. Fluid Flow Equation

According to the mechanics through porous media, the fluid flow equation in thejth fracture parallel to x axis is
(42)πλ∂ϕ∂t=1x∂∂xη1+ξxwxjkfxjμdx∂pxj∂xj=0,±1,⋯,±ny,where ξ is the axial ratio of the x axis to the y axis, ξ=Lx/Ly. η is the complete elliptic integrals of the second kind, η=∫0π/21−1−ξ2sin2sds. kfxj is the permeability of the jth fracture parallel to the x axis.Considering the flow in the fracture is in laminar, the permeability of thejth fracture parallel to the x axis is
(43)kfxj=π2hxj2pxj−σh48E2.Substitute into the fluid flow equation, it can be obtained:(44)∂ϕ∂t=π21−v2η1+ξ24λE3μxdx∂∂xxhxj3pxj−σh3∂pxj∂x.Considering the change of porosity with pressure, it is represented by rock compression coefficient, that is(45)CR=1ϕdϕdp.Integrate the above equation and the porosity can be obtained by:(46)ϕ=ϕ0eCRp−p0,where p0 is the initial pressure.Expand the above equation by series, omit the infinitesimal of higher order, and get(47)ϕ=ϕ0CRp−p0,so get
(48)∂ϕ∂t=ϕ0CR∂p∂t.Thus, in thejth fracture parallel to the x axis, the fluid flow equation is
(49)∂pxj∂t=Axx∂∂xBxpxj,x∂pxj∂xj=0,±1,⋯,±ny,where
(50)Axx=π21−v2η1+ξ24ϕ0CRλE3μxdx,Bxpxj,x=xhxj3pxj−σh3.The initial conditions and boundary conditions are(51)pxjt=0=σhpxjx=Lxj=σh,∂pxj∂xx=Lxj=0.Similarly, in theith fracture parallel to the y axis, the fluid flow equation is
(52)∂pyi∂t=Ayy∂∂yBypyi,y∂pyi∂yi=0,±1,⋯,±nx,where
(53)Ayy=π21−v2η1+ξ24ϕ0CRλE3μydy,Bypyi,y=yhyi3pyi−σH3.The initial conditions and boundary conditions are(54)pyit=0=σHpyiy=Lyi=σH,∂pyi∂yy=Lyi=0.


### 3.1.5. Coupling Conditions

Each fracture of the network (whether parallel tox axis or y axis) is relatively independent and interrelated, constituting a complex network fracture system. Each fracture is described by its own fluid flow equation, fracture width equation, and fracture height equation and is controlled by volume balance equation. At the same time, there is coupling condition; that is, the pressure at the origin of coordinates is equal, expressed as
(55)px0x=0=py0y=0.


## 3.1.1. Volumetric Equilibrium Equation

According to the principle of mass conservation and ignoring fluid compressibility and initial filtration loss of fracturing fluid, the volume of fracturing injection is equal to the sum of the volume of network fractures and the filtration loss volume of fracturing fluid, that is(23)∫0tqτdτ=Vfxt+Vfyt−Vlxt+Vlyt,where Vfx and Vfy are the volume of fractures parallel to x axis and y axis, respectively. Vlx and Vly are the filtration volume of fracturing fluid for fractures parallel to x axis and y axis, respectively.According to the morphology and composition of the fracture network, the cross-section of fractures parallel to thex axis is approximately elliptical, and the volume of fractures is
(24)Vfx=∑j=−nyny∫−LxjLxjπ4wxjhxjdx,where wxj is the fracture width of the jth fracture parallel to the x axis. hxj is the fracture height of the jth fracture parallel to the x axis.Similarly, the volume of the fracture parallel to they axis is
(25)Vfy=∑i=−nxnx∫−LyiLyiπ4wyihyidy,where wyi is the fracture width of the ith fracture parallel to the y axis. hyi is the fracture height of the ith fracture parallel to the y axis.The volume of the fracture is equal to the sum of the volume of fractures parallel to thex axis and the volume of fractures parallel to the y axis, that is:
(26)Vf=∑j=−nyny∫−LxjLxjπ4wxjhxjdx+∑i=−nxnx∫−LyiLyiπ4wyihyidy.According to the morphology and composition of the fracture network, the fracturing fluid filtration volume of fractures parallel to thex axis can also be obtained by the following equation:
(27)Vlx=∑j=−nyny∫−LxjLxj∫0t2hxjCLt−τ0sdsdx.The fracturing fluid filtration volume of the fracture parallel to they axis is
(28)Vly=∑i=−nxnx∫−LyiLyi∫0t2hyiCLt−τ0sdsdy.


## 3.1.2. Fracture Width Equation

In thejth fracture parallel to the x axis, considering the variation of crustal stress in the longitudinal direction, when the fracture penetrates the upper and lower layers (hxj≥Hp), the net pressure distribution in the fracture is
(29)Δpxjz=pxj−σhuaxj≤z≤lxjpxj−σhbxj≤z≤axjpxj−σhl−lxj≤z≤bxjj=0,±1,⋯,±ny,where
(30)lxj=hxj2=hxju+hxjl2,axj=Hp+hxjl−hxju2,bxj=−Hp−hxjl+hxju2.According to the elastic deformation theory, the net pressure in the fracture is decomposed into the sum of even stress distribution and odd stress distribution, and then, the fracture width profile is calculated by the England and Green formulas. After deduction, the fracture width is(31)wxj=41−v2lxjEpxj−σh−41−v2lxjπEσhu−σhcos−1axjlxj−axjlxjlnlxj+lxj2−axj2axj−41−v2lxjπEσhl−σhcos−1bxjlxj−bxjlxjlnlxj+lxj2−bxj2bxjj=0,±1,⋯,±ny.When the fracture does not penetrate the overlying and underlying strata (hxj<Hp) or the longitudinal variation of crustal stress is ignored, the fracture width of the jth fracture parallel to the x axis is
(32)wxj=21−v2hxjEpxj−σhj=0,±1,⋯,±ny.Similarly, when the fracture penetrates the upper and lower layers (hyi≥Hp) or the longitudinal variation of crustal stress is ignored, the fracture width of the ith fracture parallel to the y axis is
(33)wyi=41−v2lyiEpyi−σH−41−v2lyiπEσHu−σHcos−1ayilyi−ayilyilnlyi+lyi2−ayi2ayi−41−v2lyiπEσHl−σHcos−1byilyi−byilyilnlyi+lyi2−byi2byii=0,±1,⋯,±nx,where
(34)lyi=hyi2=hyiu+hyil2,ayi=Hp+hyil−hyiu2,byi=−Hp−hyil+hyiu2.When the fracture does not penetrate the upper and lower layers (hyi<Hp) or the longitudinal variation of crustal stress is ignored, the fracture width of the ith fracture parallel to the y axis is
(35)wyi=21−v2hyiEpyi−σHi=0,±1,⋯,±nx.


## 3.1.3. Fracture Height Equation

According to the theory of linear elastic fracture mechanics, when the fracture penetrates the upper and lower layers (hxj≥Hp), the stress intensity factor at the upper and lower ends of the jth fracture parallel to the x axis is
(36)KIu=1πlxj∫−lxjlxjΔpxjlxj+zlxj−zdz,(37)KIl=−1πlxj∫lxj−lxjΔpxjlxj−zlxj+zdz.According to the net pressure distribution in the fracture and letKIu=KICu and KIl=KICl, the fracture height equation of the jth fracture parallel to the x axis is
(38)πlxjKICu=πlxjpxj−σhu+σhl2+σhl−σhlxj2−bxj2−lxjsin−1bxjlxj−σhu−σhlxj2−axj2−lxjsin−1axjlxj,πlxjKICl=πlxjpxj−σhu+σhl2−σhl−σhlxj2−bxj2+lxjsin−1bxjlxj+σhu−σhlxj2−axj2+lxjsin−1axjlxj.When the fracture does not penetrate the upper and lower layers (hxj<Hp) or the longitudinal variation of crustal stress is ignored, the fracture height equation of the jth fracture parallel to the x axis is
(39)KIC=pxj−σhπlxjj=0,±1,⋯,±ny.Similarly, when the fracture penetrates the upper and lower layers (hyi≥Hp), the fracture height equation of the ith fracture parallel to the y axis is
(40)πlyiKICu=πlyipyi−σHu+σHl2+σHl−σHlyi2−byi2−lyisin−1byilyi−σHu−σHlyi2−ayi2−lyisin−1ayilyi,πlyiKICl=πlyipyi−σHu+σHl2−σHl−σHlyi2−byi2+lyisin−1byilyi+σHu−σHlyi2−ayi2+lyisin−1ayilyi.When the fracture does not penetrate the upper and lower layers (hyi<Hp) or the longitudinal variation of crustal stress is ignored, the fracture height equation of the ith fracture parallel to the y axis is
(41)KIC=pyi−σHπlyii=0,±1,⋯,±nx.


## 3.1.4. Fluid Flow Equation

According to the mechanics through porous media, the fluid flow equation in thejth fracture parallel to x axis is
(42)πλ∂ϕ∂t=1x∂∂xη1+ξxwxjkfxjμdx∂pxj∂xj=0,±1,⋯,±ny,where ξ is the axial ratio of the x axis to the y axis, ξ=Lx/Ly. η is the complete elliptic integrals of the second kind, η=∫0π/21−1−ξ2sin2sds. kfxj is the permeability of the jth fracture parallel to the x axis.Considering the flow in the fracture is in laminar, the permeability of thejth fracture parallel to the x axis is
(43)kfxj=π2hxj2pxj−σh48E2.Substitute into the fluid flow equation, it can be obtained:(44)∂ϕ∂t=π21−v2η1+ξ24λE3μxdx∂∂xxhxj3pxj−σh3∂pxj∂x.Considering the change of porosity with pressure, it is represented by rock compression coefficient, that is(45)CR=1ϕdϕdp.Integrate the above equation and the porosity can be obtained by:(46)ϕ=ϕ0eCRp−p0,where p0 is the initial pressure.Expand the above equation by series, omit the infinitesimal of higher order, and get(47)ϕ=ϕ0CRp−p0,so get
(48)∂ϕ∂t=ϕ0CR∂p∂t.Thus, in thejth fracture parallel to the x axis, the fluid flow equation is
(49)∂pxj∂t=Axx∂∂xBxpxj,x∂pxj∂xj=0,±1,⋯,±ny,where
(50)Axx=π21−v2η1+ξ24ϕ0CRλE3μxdx,Bxpxj,x=xhxj3pxj−σh3.The initial conditions and boundary conditions are(51)pxjt=0=σhpxjx=Lxj=σh,∂pxj∂xx=Lxj=0.Similarly, in theith fracture parallel to the y axis, the fluid flow equation is
(52)∂pyi∂t=Ayy∂∂yBypyi,y∂pyi∂yi=0,±1,⋯,±nx,where
(53)Ayy=π21−v2η1+ξ24ϕ0CRλE3μydy,Bypyi,y=yhyi3pyi−σH3.The initial conditions and boundary conditions are(54)pyit=0=σHpyiy=Lyi=σH,∂pyi∂yy=Lyi=0.


## 3.1.5. Coupling Conditions

Each fracture of the network (whether parallel tox axis or y axis) is relatively independent and interrelated, constituting a complex network fracture system. Each fracture is described by its own fluid flow equation, fracture width equation, and fracture height equation and is controlled by volume balance equation. At the same time, there is coupling condition; that is, the pressure at the origin of coordinates is equal, expressed as
(55)px0x=0=py0y=0.


## 3.2. Numerical Solution of the Fracture Network Model

### 3.2.1. Difference Method for Fluid Flow Equations

It is a nonlinear partial differential equation for the fluid flow equation of thejth fracture parallel to the x axis. We use the implicit scheme difference method to solve it, and the nonlinear problem is solved by the coefficient freezing method. Thus, the left end of the fluid flow equation is discretized as follows:
(56)∂pxj∂t≈pxjkn+1−pxjknΔt.The right end of the fluid flow equation is discretized as follows:(57)Ax∂∂xBx∂pxj∂x≈AxknBxk+1/2npxjk+1n+1−pxjkn+1/Δx−Bxk−1/2npxjkn+1−pxjk−1n+1/ΔxΔx,where Δt is the time step and Δx is the length step in the x direction.
(58)Bxk+1/2n≈Bxk+1n+Bxkn2,Bxk−1/2n≈Bxkn+Bxk−1n2.Thus, the difference equation of the fluid flow equation is(59)−ΔtAxknBxk+1/2nΔx2pxjk+1n+1+1+ΔtAxknBxk+1/2n+Bxk−1/2nΔx2pxjkn+1−ΔtAxknBxk−1/2nΔx2pxjk−1n+1=pxjkn.These are tridiagonal linear equations with dominant diagonal, which can be solved by the catch-up method combined with their initial conditions and boundary conditions.Similarly, for theith fracture parallel to the y axis, the difference equation of the fluid flow equation (3-108) can be expressed as follows:
(60)−ΔtAyknByk+1/2nΔy2pyik+1n+1+1+ΔtAyknByk+1/2n+Byk−1/2nΔy2pyikn+1−ΔtAyknByk−1/2nΔy2pyik−1n+1=pyikn,where Δy is the length step in the y direction.


### 3.2.2. Forward Algorithm for the Numerical Solution

The main purpose of establishing and solving the fracture network model is known basic parameters (such as elastic modulusE, Poisson’s ratio ν and other rock parameters, fracturing fluid viscosity μ, injection rate q, injection time T, and other construction parameters); by solving the model, the fracture network parameters can be obtained (such as half-length Lx and Ly, fracture width, fracture height, and Vf).It is not hard to see that each fracture of the network is described by its own fluid flow equation, fracture height equation, and fracture width equation, which is a complete mathematical model. In other words, under the condition of known basic parameters, fracture spacing (dx and dy) and fracture half-length (Lx and Ly), the fracture length, fracture width, and fracture height and pressure in each fracture can be calculated by using this model. Therefore, by using the volume balance equation and coupling conditions (the pressure is equal at the intersection), the algorithm is established to obtain some parameters of the fracture spacing and the fracture half-length.The forward algorithm in this paper is to obtain the half-length (Lx and Ly) of fracture network, the length, width, height, pressure, and Vf of each fracture under the condition of known basic parameters and fracture spacing (dx and dy).The specific step of the forward algorithm is as follows:(1)
The coupling conditions are transformed to construct the following function:(61)Jd1Lx,Ly=px0x=0−py0y=0.As you can see that for anyLx and , Jd1 can be calculated by the function. Theoretical analysis and practical calculation show that if Jd1=0, then Lx =Ly. Therefore, for a fixed value Lx, according to the condition Jd1=0, it is easy to calculate the value of Ly by iterative method. At this time, the iteration formula is designed as follows:
(62)Ly=1±λLy,where λ is the iteration factor which is set to 0.1 during trial calculation.
(2)
According to the volume equilibrium equation, the following function is constructed:(63)Jd2Lx=∫0tqτdτ−Vfxt+Vfyt+Vlxt+Vlyt.It is not difficult to see that for any value ofLx, since the value of Ly is determined by the condition Jd1=0, Jd2 is just a function of Lx. Theoretical analysis and practical calculation show that if Jd2=0, then the value of Lx reaches volumetric balance. Therefore, according to the condition Jd2=0, it is easy to calculate the value of Lx by iterative method. At this time, the iteration formula is designed as follows:
(64)Lx=1±λLx.(3)
The forward algorithm advances step by step by time stepΔt. Because the implicit difference method is used to solve the fluid flow equation, the numerical calculation is unconditionally stable. However, the smaller the time step, the higher the accuracy of the calculation resultTo sum up, the forward algorithm designed is shown in Figure5.Figure 5
Flow chart of the forward algorithm for numerical solution.


### 3.2.3. Inverse Algorithm for the Numerical Solution

The forward algorithm for the numerical solution of the network fracture model requires the fracture spacing (dx and dy) to be input as known parameters. Usually, the forward algorithm is used to check and invert the microseismic data of nearby wells, borrowing the fracture spacing of adjacent wells.In fact, parameter inversion is the inverse problem of parameter forward modeling, and the inversion algorithm is different with different parameters. The inversion algorithm in this paper is to obtain the fracture spacing (dx and dy), the fracture length, the fracture width, the fracture height, the pressure in the fracture, and Vf of each fracture through the algorithm under the condition that the basic parameters and the half-length (Lx and Ly) of the fracture network are known.The inversion algorithm of the fracture network numerical solution is similar to the forward algorithm, which is briefly described here. Take the coupling conditions of equal pressure at the origin of coordinates and construct the following function:(65)JL1dx,dy=px0x=0−py0y=0.The following iterative formula is adopted:(66)dy=1±λdy.At the same time, according to the volume equilibrium equation, the following function is constructed:(67)JL2dx=∫0tqτdτ−Vfxt+Vfyt+Vlxt+Vlyt.The following iterative formula is adopted:(68)dx=1±λdx.Therefore, the process of the inversion algorithm is shown in Figure6.Figure 6
Flow chart of the inverse algorithm for numerical solution.


## 3.2.1. Difference Method for Fluid Flow Equations

It is a nonlinear partial differential equation for the fluid flow equation of thejth fracture parallel to the x axis. We use the implicit scheme difference method to solve it, and the nonlinear problem is solved by the coefficient freezing method. Thus, the left end of the fluid flow equation is discretized as follows:
(56)∂pxj∂t≈pxjkn+1−pxjknΔt.The right end of the fluid flow equation is discretized as follows:(57)Ax∂∂xBx∂pxj∂x≈AxknBxk+1/2npxjk+1n+1−pxjkn+1/Δx−Bxk−1/2npxjkn+1−pxjk−1n+1/ΔxΔx,where Δt is the time step and Δx is the length step in the x direction.
(58)Bxk+1/2n≈Bxk+1n+Bxkn2,Bxk−1/2n≈Bxkn+Bxk−1n2.Thus, the difference equation of the fluid flow equation is(59)−ΔtAxknBxk+1/2nΔx2pxjk+1n+1+1+ΔtAxknBxk+1/2n+Bxk−1/2nΔx2pxjkn+1−ΔtAxknBxk−1/2nΔx2pxjk−1n+1=pxjkn.These are tridiagonal linear equations with dominant diagonal, which can be solved by the catch-up method combined with their initial conditions and boundary conditions.Similarly, for theith fracture parallel to the y axis, the difference equation of the fluid flow equation (3-108) can be expressed as follows:
(60)−ΔtAyknByk+1/2nΔy2pyik+1n+1+1+ΔtAyknByk+1/2n+Byk−1/2nΔy2pyikn+1−ΔtAyknByk−1/2nΔy2pyik−1n+1=pyikn,where Δy is the length step in the y direction.


## 3.2.2. Forward Algorithm for the Numerical Solution

The main purpose of establishing and solving the fracture network model is known basic parameters (such as elastic modulusE, Poisson’s ratio ν and other rock parameters, fracturing fluid viscosity μ, injection rate q, injection time T, and other construction parameters); by solving the model, the fracture network parameters can be obtained (such as half-length Lx and Ly, fracture width, fracture height, and Vf).It is not hard to see that each fracture of the network is described by its own fluid flow equation, fracture height equation, and fracture width equation, which is a complete mathematical model. In other words, under the condition of known basic parameters, fracture spacing (dx and dy) and fracture half-length (Lx and Ly), the fracture length, fracture width, and fracture height and pressure in each fracture can be calculated by using this model. Therefore, by using the volume balance equation and coupling conditions (the pressure is equal at the intersection), the algorithm is established to obtain some parameters of the fracture spacing and the fracture half-length.The forward algorithm in this paper is to obtain the half-length (Lx and Ly) of fracture network, the length, width, height, pressure, and Vf of each fracture under the condition of known basic parameters and fracture spacing (dx and dy).The specific step of the forward algorithm is as follows:(1)
The coupling conditions are transformed to construct the following function:(61)Jd1Lx,Ly=px0x=0−py0y=0.As you can see that for anyLx and , Jd1 can be calculated by the function. Theoretical analysis and practical calculation show that if Jd1=0, then Lx =Ly. Therefore, for a fixed value Lx, according to the condition Jd1=0, it is easy to calculate the value of Ly by iterative method. At this time, the iteration formula is designed as follows:
(62)Ly=1±λLy,where λ is the iteration factor which is set to 0.1 during trial calculation.
(2)
According to the volume equilibrium equation, the following function is constructed:(63)Jd2Lx=∫0tqτdτ−Vfxt+Vfyt+Vlxt+Vlyt.It is not difficult to see that for any value ofLx, since the value of Ly is determined by the condition Jd1=0, Jd2 is just a function of Lx. Theoretical analysis and practical calculation show that if Jd2=0, then the value of Lx reaches volumetric balance. Therefore, according to the condition Jd2=0, it is easy to calculate the value of Lx by iterative method. At this time, the iteration formula is designed as follows:
(64)Lx=1±λLx.(3)
The forward algorithm advances step by step by time stepΔt. Because the implicit difference method is used to solve the fluid flow equation, the numerical calculation is unconditionally stable. However, the smaller the time step, the higher the accuracy of the calculation resultTo sum up, the forward algorithm designed is shown in Figure5.Figure 5
Flow chart of the forward algorithm for numerical solution.


## 3.2.3. Inverse Algorithm for the Numerical Solution

The forward algorithm for the numerical solution of the network fracture model requires the fracture spacing (dx and dy) to be input as known parameters. Usually, the forward algorithm is used to check and invert the microseismic data of nearby wells, borrowing the fracture spacing of adjacent wells.In fact, parameter inversion is the inverse problem of parameter forward modeling, and the inversion algorithm is different with different parameters. The inversion algorithm in this paper is to obtain the fracture spacing (dx and dy), the fracture length, the fracture width, the fracture height, the pressure in the fracture, and Vf of each fracture through the algorithm under the condition that the basic parameters and the half-length (Lx and Ly) of the fracture network are known.The inversion algorithm of the fracture network numerical solution is similar to the forward algorithm, which is briefly described here. Take the coupling conditions of equal pressure at the origin of coordinates and construct the following function:(65)JL1dx,dy=px0x=0−py0y=0.The following iterative formula is adopted:(66)dy=1±λdy.At the same time, according to the volume equilibrium equation, the following function is constructed:(67)JL2dx=∫0tqτdτ−Vfxt+Vfyt+Vlxt+Vlyt.The following iterative formula is adopted:(68)dx=1±λdx.Therefore, the process of the inversion algorithm is shown in Figure6.Figure 6
Flow chart of the inverse algorithm for numerical solution.


## 4. Results and Discussion

### 4.1. Application

Take the well H03-11 as an example. The proposed continuous fracture network mechanics and model from this research have been applied to this well. The main parameters used during the simulation process for well H03-11 are listed in Tables2 and 3.Table 2
Reservoir and construction parameters.
ParameterValueDepth (m)1837Permeability (10−3 μm2)0.19Porosity (%)8.25Thickness (m)12Poisson’s ratio0.24Elastic modulus (MPa)32891Fracture spacing ofx axis (m)32.6Fracture spacing ofy axis (m)18.7Fracturing fluid viscosity (mPa.s)1.1Filtration coefficient (m/min0.5)0.0012Injection rate (m3/min)12Injection time of fracturing fluid (min)126Table 3
Rock mechanics parameters.
ParameterOverburden layerProduction layerUnderburden layerHorizontal minimum principal stress (MPa)29.5028.3431.63Horizontal maximum principal stress (MPa)33.8932.0635.72The formation of this well belongs to a naturally fractured developed reservoir, with the value of net pressure in the fracture larger thanσH−σh+St. According to the mechanical mechanism of fracture network formation proposed in this paper, well H03-11 can form a fracture network. Therefore, the continuous fracture network model established in this paper is appropriate for simulations of fracture propagation within well H03-11. Using the continuous fracture network model, the growth of fracture height was calculated. As shown in Figure 7, the heights of the upper and lower extensions are 19.53 m and 11.71 m, respectively. Therefore, the height of the fracture is 31.24 m.Figure 7
Fracture height profile extension.Lx and Ly are the half-length of the fracture in the x and y directions, respectively. H is the fracture height. The calculated results of the continuous fracture network are shown in Figure 8. 2Lx is 335.42 m, 2Ly is 105.72 m, and H is 31.24 m.Figure 8
Numerical simulation results of the continuous fracture network. The change of color in the figure represents the change of fracture height.


### 4.2. Verification and Discussion

The height of a fracture network changes with time as shown in Figure9. The total injection time of well H03-11 is 126 minutes. At the early time stage, the fracture height increases rapidly. The later-time fracture height grows slowly. In the production process, the height changes, so considering height variation, the simulation results are more accurate.Figure 9
Relationship between height and time.We discuss the changes ofLx, Ly, and Vf with different times; the results are shown in Figure 10. It can be seen from Figure 10 that the Lx and Ly increase rapidly at the early time stage and grow slowly at later-time. The growth of Vf is S-shaped; that is, the first growth is slow, the middle growth is fast, and the latter is slow.Figure 10
Change ofLx, Ly, and Vf with time.We also discussed the influence of fracturing fluid viscosity, filtration coefficient, and injection rate onLx, Ly, and Vf. The results are shown in Figure 11. It can be seen from Figure 11 that (1) the higher the viscosity and injection rate are, the better the effect of reservoir simulation is. The larger the filtration coefficient is, the worse the effect of reservoir simulation is. (2) Viscosity has little effect on Lx and Ly, but a greater effect on Vf, only because viscosity has a greater effect on fracture height. (3) The injection rate has a great influence on Vf, indicating that the injection rate should be increased in field case under the premise of ensuring safety.Figure 11
(a) Change ofLx, Ly, and Vf with fracturing fluid viscosity at 126 min. (b) Change of Lx, Ly, and Vf with filtration coefficient at 126 min. (c) Change of Lx, Ly, and Vf with injection rate at 126 min.
(a)(b)(c)Microseismic monitoring technology is the main method for accurate fracture network design. In order to verify the accuracy of the continuous fracture network model in this paper, the fracture was monitored by microseismic monitoring technology; the results are shown in Figure12. An industry-accepted technology, Meyer, was employed to predict the dimension of fractures with the same well H03-11 data. Calculated by the forward algorithm of continuous fracture network model, Meyer and the microseismic monitoring results are compared in Table 4.Figure 12
Results of microseismic monitoring of well H03-11.Table 4
Comparison of prediction results of well H03-11 by different technologies.
Applied technology2Lx2LyHResult (m)Error (%)Result (m)Error (%)Result (m)Error (%)Microseismic monitoring32411429Meyer362.2911.8296.7115.1624.6215.10Continuous fracture network335.423.52105.727.2631.247.72As shown in Table4, the prediction accuracy of the continuous fracture network is higher than that of Meyer. This is because the influence of fracture height and filtration coefficient is not taken into account in the Meyer. From the discussion in Figure 11, we know that these two factors have a great influence on the prediction results.In order to further prove the prediction accuracy of the continuous fracture network model in the field case, we applied this model to a horizontal well TP12-06. The model in this paper is only used to predict the size of fractures, and the directions of fractures are calculated by the pressure field model. The results of microseismic monitoring of TP12-06 are shown in Figure13(a). The comparison between the prediction results of microseismic monitoring and continuous fracture network model is shown in Figure 13(b) and Table 5. The average error of prediction of fracture length and fracture width is 3.64% and 5.70%, respectively.Figure 13
(a) Results of microseismic monitoring of TP12-06. (b) Comparison of prediction results of microseismic monitoring and continuous fracture network model.
(a)(b)Table 5
Comparison of prediction results of well TP12-06 by different technologies.
ClusterFracture lengthFracture widthAzimuth (°)Microseismic monitoring (m)Simulation result (m)Error (%)Microseismic monitoring (m)Simulation result (m)Error (%)1372384.473.35142134.675.16712304317.834.55151142.895.37333291302.583.98193181.495.96494302309.592.51201191.734.61365332343.723.53143136.644.45566314326.393.94210191.828.6656Average3.645.70


## 4.1. Application

Take the well H03-11 as an example. The proposed continuous fracture network mechanics and model from this research have been applied to this well. The main parameters used during the simulation process for well H03-11 are listed in Tables2 and 3.Table 2
Reservoir and construction parameters.
ParameterValueDepth (m)1837Permeability (10−3 μm2)0.19Porosity (%)8.25Thickness (m)12Poisson’s ratio0.24Elastic modulus (MPa)32891Fracture spacing ofx axis (m)32.6Fracture spacing ofy axis (m)18.7Fracturing fluid viscosity (mPa.s)1.1Filtration coefficient (m/min0.5)0.0012Injection rate (m3/min)12Injection time of fracturing fluid (min)126Table 3
Rock mechanics parameters.
ParameterOverburden layerProduction layerUnderburden layerHorizontal minimum principal stress (MPa)29.5028.3431.63Horizontal maximum principal stress (MPa)33.8932.0635.72The formation of this well belongs to a naturally fractured developed reservoir, with the value of net pressure in the fracture larger thanσH−σh+St. According to the mechanical mechanism of fracture network formation proposed in this paper, well H03-11 can form a fracture network. Therefore, the continuous fracture network model established in this paper is appropriate for simulations of fracture propagation within well H03-11. Using the continuous fracture network model, the growth of fracture height was calculated. As shown in Figure 7, the heights of the upper and lower extensions are 19.53 m and 11.71 m, respectively. Therefore, the height of the fracture is 31.24 m.Figure 7
Fracture height profile extension.Lx and Ly are the half-length of the fracture in the x and y directions, respectively. H is the fracture height. The calculated results of the continuous fracture network are shown in Figure 8. 2Lx is 335.42 m, 2Ly is 105.72 m, and H is 31.24 m.Figure 8
Numerical simulation results of the continuous fracture network. The change of color in the figure represents the change of fracture height.


## 4.2. Verification and Discussion

The height of a fracture network changes with time as shown in Figure9. The total injection time of well H03-11 is 126 minutes. At the early time stage, the fracture height increases rapidly. The later-time fracture height grows slowly. In the production process, the height changes, so considering height variation, the simulation results are more accurate.Figure 9
Relationship between height and time.We discuss the changes ofLx, Ly, and Vf with different times; the results are shown in Figure 10. It can be seen from Figure 10 that the Lx and Ly increase rapidly at the early time stage and grow slowly at later-time. The growth of Vf is S-shaped; that is, the first growth is slow, the middle growth is fast, and the latter is slow.Figure 10
Change ofLx, Ly, and Vf with time.We also discussed the influence of fracturing fluid viscosity, filtration coefficient, and injection rate onLx, Ly, and Vf. The results are shown in Figure 11. It can be seen from Figure 11 that (1) the higher the viscosity and injection rate are, the better the effect of reservoir simulation is. The larger the filtration coefficient is, the worse the effect of reservoir simulation is. (2) Viscosity has little effect on Lx and Ly, but a greater effect on Vf, only because viscosity has a greater effect on fracture height. (3) The injection rate has a great influence on Vf, indicating that the injection rate should be increased in field case under the premise of ensuring safety.Figure 11
(a) Change ofLx, Ly, and Vf with fracturing fluid viscosity at 126 min. (b) Change of Lx, Ly, and Vf with filtration coefficient at 126 min. (c) Change of Lx, Ly, and Vf with injection rate at 126 min.
(a)(b)(c)Microseismic monitoring technology is the main method for accurate fracture network design. In order to verify the accuracy of the continuous fracture network model in this paper, the fracture was monitored by microseismic monitoring technology; the results are shown in Figure12. An industry-accepted technology, Meyer, was employed to predict the dimension of fractures with the same well H03-11 data. Calculated by the forward algorithm of continuous fracture network model, Meyer and the microseismic monitoring results are compared in Table 4.Figure 12
Results of microseismic monitoring of well H03-11.Table 4
Comparison of prediction results of well H03-11 by different technologies.
Applied technology2Lx2LyHResult (m)Error (%)Result (m)Error (%)Result (m)Error (%)Microseismic monitoring32411429Meyer362.2911.8296.7115.1624.6215.10Continuous fracture network335.423.52105.727.2631.247.72As shown in Table4, the prediction accuracy of the continuous fracture network is higher than that of Meyer. This is because the influence of fracture height and filtration coefficient is not taken into account in the Meyer. From the discussion in Figure 11, we know that these two factors have a great influence on the prediction results.In order to further prove the prediction accuracy of the continuous fracture network model in the field case, we applied this model to a horizontal well TP12-06. The model in this paper is only used to predict the size of fractures, and the directions of fractures are calculated by the pressure field model. The results of microseismic monitoring of TP12-06 are shown in Figure13(a). The comparison between the prediction results of microseismic monitoring and continuous fracture network model is shown in Figure 13(b) and Table 5. The average error of prediction of fracture length and fracture width is 3.64% and 5.70%, respectively.Figure 13
(a) Results of microseismic monitoring of TP12-06. (b) Comparison of prediction results of microseismic monitoring and continuous fracture network model.
(a)(b)Table 5
Comparison of prediction results of well TP12-06 by different technologies.
ClusterFracture lengthFracture widthAzimuth (°)Microseismic monitoring (m)Simulation result (m)Error (%)Microseismic monitoring (m)Simulation result (m)Error (%)1372384.473.35142134.675.16712304317.834.55151142.895.37333291302.583.98193181.495.96494302309.592.51201191.734.61365332343.723.53143136.644.45566314326.393.94210191.828.6656Average3.645.70


## 5. Conclusion

This work contributes to a novel continuous fracture network model. The findings from this study allow the following conclusions to be drawn:(1)
By comprehensively considering the development situations of natural fractures, this study provides mechanical condition requirements related to fracture network formation. The mechanical condition requirements for a fracture network in a naturally fractured reservoir lie in net pressure in a constructed fracture exceeding the horizontal principal stress difference of the reservoir. For a fracture network in a reservoir without natural fracture development, the condition requirements lie in the fracture breaking in the rock body and the value of net pressure in the fracture being larger than the sum of horizontal principal stress difference and reservoir tensile strength(2)
A continuous fracture network model was established based on the mechanical mechanism for fracture network formation. A geometric fracture network model was constructed with the main fracture and branch fractures, mathematical equations were established to simulate geometric parameters of a shale gas continuous fracture network, and related calculation methods were derived(3)
In this study, a numerical simulation program of a continuous fracture network was developed, overall consideration was made for an increase in fracture height and filtration of fracturing liquid, and a simulation study was performed on a continuous fracture network. Actual microseismic monitoring data verifies the accuracy of this model. The simulation results are of great significance for fracturing design optimization of shale gas fracture networks


## 6. Recommendations

In future research, the following aspects can be further explored:(1)
The junction of fracture is assumed to be noninterference, but there is fracture interference in the actual situation, so this study needs to be added in the future(2)
It is assumed that fractures are perpendicular to each other. However, for an oblique fracture, how to establish a model in accordance with the real situation is worthy of further discussion



---

*Source: 4026200-2022-04-19.xml*