# Phase Field Modeling of Multiple Fracture Growth in Natural Fractured Reservoirs

Authors: QiangSheng He; Chuang Liu
Journal: Geofluids
 (2023)

---

## Abstract

In recent years, hydraulic fracturing techniques have been widely used in extracting unconventional reservoir resources. During multicluster fracturing process, the stress shadowing effect can lead to a nonuniform distribution of fracturing fluid. In this paper, a two-dimensional multiple fracture propagation model is developed based on phase field method considering the distribution of fracturing fluid within each fracture during fracturing process. The distribution of fracturing fluid injected into each fracture is calculated through solving governing equations of perforation friction as well as wellbore friction. Numerical results demonstrate that the natural fractures in the reservoir have a significant influence on the distribution of injected fluid and the morphology of fracture network. In addition, higher injection rates and higher fluid viscosity stimulations are in favor of inducing uniform distribution of injected fluid into each cluster. The presented numerical model provides an effective tool to extend our understandings in generating fracture networks.


---

## Body

## 1. Introduction

There is a large amount of natural gas resources in ultralow permeability reservoirs such as shale [1]. Hydraulic fracturing technique combined with horizontal well is utilized to improve reservoir permeability through creating fracture network [2–4]. Previous results have verified that the natural fractures have a significant influence on the propagation morphology of hydraulic fractures [5–8]. Fracturing design strategies for conventional reservoirs are no longer applicable to naturally fractured reservoirs. Natural fractures are important storage spaces for oil and gas [9–11]. The presence of natural fractures also contributes to the production of complex fracture networks in stimulation treatments [12]. It has been shown that approximately 55% of natural hydrocarbon output is attributed to the activation of natural fractures [13]. In addition, the presence of natural fractures also has an impact on the distribution of fracturing fluid within multiple fractures. Therefore, studying the interaction between hydraulic fractures and natural fractures is significant for understanding complex fracture networks.In recent years, more and more researchers have studied the propagation behavior of hydraulic fractures when they intersect with natural fractures. [14, 15] carried out laboratory tests to discover the effect of stress differences and approach angles on hydraulic fracture and natural fracture propagation morphology. [16] studied the effect of reservoir discontinuities on hydraulic fracture extension based on mine recovery experiments and laboratory tests and proposed a W-T criterion. [17] considered the effect of different natural fracture angels and different wellbore spacing on the expansion path of hydraulic fractures. [18] used the extended finite element method (XFEM) to investigate the damage patterns and fracture propagation mechanisms at the intersection of hydraulic fracture and natural fracture. [19] used the finite discrete element method (FDEM) to investigate the basic laws of shear and tensile damage of natural fractures under the disturbing stresses of hydraulic fractures. However, these studies only consider the interaction between individual hydraulic fractures and natural fractures but do not consider the effect of stress interference between multiple fractures on fracture propagation. With regard to multicluster hydraulic fracturing, it is difficult to be achieved by means of laboratory tests. Numerical methods can effectively simulate the effect of multicluster fracture stress disturbance on hydraulic fracture propagation. [20] used a discrete fracture network (DFN) model to simulate multicluster hydraulic fractures in the presence of naturally fractured reservoirs. [21] used displacement discontinuity method (DDM) to build a three-dimensional model of multicluster hydraulic fracture propagation. Investigations were done into how the shape of natural fractures and the horizontal stress difference affected fracture morphology. A particle-based model was created using XSite numerical simulation software by [22] to study multicluster hydraulic fracture propagation within natural fractures. [23] proposed a fully coupled hydraulic fracturing model to study the influence of preexisting natural fractures on the direction, length, and width distribution of the fracture network. However, these numerical studies do not consider the effect of the mechanical interactions between multiple hydraulic fractures and natural fractures on the distribution of injected fluids into each fracture. Neglecting the nonuniform distribution of the fluids would result in an inaccurate prediction of the fracture network.During the simultaneous expansion of multicluster hydraulic fracturing, the flow distribution in each fracture is often uneven. An important reason for this uneven distribution of flow is the phenomenon of stress shadowing [20, 24]. Fluid distribution in multiple fractures is inversely proportional to flow resistance, which is a function of perforation friction, wellbore friction, and fluid pressure within the fractures. Stress shadowing effect imposes additional compressive stress on the internal fractures, increasing the flow resistance of the internal fractures and causing the internal fractures to distribute less flow and more fluid into the external fractures. To study the effects of stress shadowing and natural fractures on hydraulic fracture propagation during multifracture propagation, [25] uses a finite entry design technique to achieve simultaneous and uniform propagation of multiple fractures in the model by varying the diameter of the injection holes and the number of injection holes. [26] built a three-dimensional model to consider the effect of stress shadowing and multihorizontal well interwell stress interference on hydraulic fracture propagation in multiple fracture propagation. [27] used a complex fracture propagation model (XFRAC) to study the effect of group construction engineering parameters on the uniform propagation of fractures. [28] found that the stress shadowing effect caused by the interaction between multiple fractures increases the effective closure stress. Based on numerical simulations, [29] concluded that in the case of collinear cracks, there is a significant effect of crack spacing on the duration of crack propagation, and the smaller the spacing of fractures, the shorter the duration of crack propagation. In this paper, the phase field method (PFM) is used to simulate the simultaneous growth of multiple hydraulic fractures. The PFM is widely used to simulate crack growth due to its advantages [30–32]. It has been applied in modelling hydraulic fracture propagation in porous media [33, 34]. However, the effect of stress shadowing and natural fractures on the nonuniform fluid distribution during stimulation process based on PFM is yet to be estimated. In this study, we solve the fluid distribution in multiple clusters by utilizing Kirchhoff’s first and second laws that describe the balance of fluid pressure between wellbore and each fracture with perforation and wellbore friction.In this paper, we built a PFM model considering nonuniform distribution of injected fluid into each cluster to investigate the hydraulic fracture growth in naturally fractured reservoirs. The numerical model is verified by two analytical solutions. Finally, we investigated the effects of different orientations of natural fractures, horizontal stress differences, fracture energy, injection rates, fracturing fluid viscosity, and natural fracture density on multiple fracture propagation.


## 2. Mathematical Model

### 2.1. Phase Field Approximation of Cracks

[35] proposed that the sharp crack topology can be approximated as a diffusion crack topology, described by using the field variable ϕ, as shown in Figure 1. ϕ=1 means the domain is completely broken, and ϕ=0 means the domain is complete. The model considers a bounded domain Ω, with an outer boundary ∂Ω and an inner boundary Γ. The crack surface density function can be expressed as
(1)γϕ,∇ϕ=ϕ22l0+l02∂ϕ∂xi∂ϕ∂xi,where l0 is a length-scale parameter that limits the width of the transition domain between a completely fractured medium and an intact medium, and its value will affect the accuracy of the results, which needs to be determined by combining experience and sensitivity analysis for different problems. [35] pointed out that the value of l0 should be less than or equal to half of the mesh size in the fracture zone. For the hydraulic fracturing problems, we refer to relevant literature to determine the value of this parameter [34, 36, 37].Figure 1
Schematic diagram of cracks. (a) Sharp crack. (b) Phase field discrete crack.In fracture variational theory, for a pure single-phasic solid, the total potential energy can be represented by elastic strain energy and crack surface energy and external force work. However, for a porous solid, the influence of the fluid pressure of the solid also needs to be considered, and the total potential energy can be expressed as(2)Ψu,Γ=∫ΩψεεdΩ−∫Ωαp·∇·udΩ+∫ΓGcdS−∫Ωb·udΩ−∫∂Ωhif·udS,where ψεε is the elastic strain energy density, Gc is the critical energy release rate, p is the fluid pressure, α is the Biot coefficient, b is the body force, and f is the boundary traction imposed on the outer boundary ∂Ω.In a solid bodyΩ with discrete internal discontinuity Γ, the surface energy of a crack can be represented by a domain integral with the crack phase field:
(3)∫ΓGcdS=∫ΩGcϕ22l0+l02∂ϕ∂xi∂ϕ∂xidΩ.According to [38], it is assumed that the phase field affects only the stretched part of the elastic energy density. Therefore, the elastic energy density can be written as
(4)ψεε=1−k1−ϕ2+kψε+ε+ψε−ε,where 0<k≪1 is the stiffness coefficient, which prevents the stretched part of the elastic energy density from being 0 and avoids the singularity of the numerical value when the phase field ϕ tends to 1, and ψ+ε and ψ−ε are the tensile and compression parts of elastic energy, which can be calculated by
(5)ε=ε++ε−,ε±=∑a=1dεa±na⊗na,ψ±ε=λ2trε±2+μtrε±2,where ε+ and ε− are the tensile and compressive strain tensors, εa and na are the main strains and their directions, the operators •+ and •− are defined as •+=•+•/2 and •−=•−•/2, and λ and μ are the Lamé constants which can be calculated by
(6)λ=Ev1+v1−2v,μ=E21+v.According to Biot’s poroelasticity theory [39], the relationship between the total Cauchy stress σij and the effective stress σije considering the initial in situ stress σ0ij can be expressed as
(7)σij=σije+σ0ij−αPf.By introducing a phase field, the total Cauchy stress across the domain can be expressed as
(8)σij=σije+1−k1−ϕ2+kσ0ij−αPf,inΩ.Combining Equation (4) and Equation (5), the effective stress can be calculated by
(9)σije=∂ψεε∂ε=1−k1−ϕ2+kλtrε+I+2με++λtrε−I+2με−,where I is the identity tensor.The displacement field equilibrium equation considering body forces and boundary conditions can be written as(10)σijnj=fti,on∂Ωt.Finally, the elasticity matrix in the numerical implementation is represented as(11)De=D1111D1122D1133D1112D1123D1113D2211D2222D2233D2212D2223D2213D3311D3322D3333D3312D3323D3313D1211D1222D1233D1212D1223D1213D2311D2322D2333D2312D2323D2313D1311D1322D1333D1312D1323D1313,with Dijkl=D¯ijkl+D~ijkl.The componentD¯ijkl is written as
(12)D¯ijkl=λ1−k1−ϕ2+kHεtrε+Hε−trεδijδkl,where Hεx is a Heaviside function; at x>0, Hεx=1, and at x<0, Hεx=0; δij and δkl denote Kronecker deltas.TheD~ijkl component can be expressed as [40]
(13)D~ijkl=2μ1−k1−ϕ2+kPijkl++Pijkl−,with
(14)Pijkl+=∑a=13∑b=13Hεεaδabnainajnbknbl+∑a=13∑b≠a312εa+−εb+εa−εbnainbjnaknbl+nbknal,Pijkl−=∑a=13∑b=13Hε−εaδabnainajnbknbl+∑a=13∑b≠a312εa−−εb−εa−εbnainbjnaknbl+nbknal,where nai denotes the ith component of vector na.


### 2.2. Phase Field Evolution Equation

According to Equation (3), the energy equation of Equation (2) can be expressed as
(15)Ψu,Γ=∫ΩψεεdΩ−∫Ωαp·∇·udΩ+∫ΩGcϕ22l0+l02∂ϕ∂xi∂ϕ∂xidΩ−∫Ωb·udΩ−∫∂Ωhif·udS.When the first derivative of the energy equationΨu,Γ is equal to zero, Equation (15) takes the minimum value, which means that the slit Γx,t has cracked and extended within time t∈0,Tx∈Ω. Therefore, the governing equation of the phase field can be obtained as
(16)∂σijpor∂xj+bi=0,2l01−kψε+Gc+1ϕ−l02∇2ϕ=2l01−kψε+Gc,where σijpor is the Cauchy stress tensor. Fracture propagation in the phase field method is driven by the tensile part of the elastic energy [41].In order to prevent crack healing when the elastomer is under pressure or unloading, it is necessary to supplement the irreversible condition, and the irreversibility of the crack propagation evolution is ensured by introducing the historical state variableHx,t.
(17)Hx,t=maxs∈0,tψε+εx,s.Combined with Equation (17), the evolution equation in Equation (16) is rewritten as
(18)2l01−kHGc+1ϕ−l02∇2ϕ=2l01−kHGc.The phase field is guaranteed to increase monotonically asH increases. However, this equation does not consider the effect of the initial ground stress. In this research, a new historical reference field was used H as follows [34]:
(19)Hx,t=maxs∈0,tψε+εx,s+σ0:εx,s,inΩ×0,T.The initial hydraulic fracture in this study is considered to be a preexisting notches, and the initial historical field is used to create an initial phase field of [42]
(20)H0x=BGc2l01−2dx,ll0,dx,l≤l02,0,dx,l>l02,where B is a scalar parameter and dx,l is the shortest distance from any point of the model to the discrete fracture.


### 2.3. Fluid Flow in Porous Media

To calculate the flow equations in the domainΩ, the entire domain is divided into three parts [43]: Ωrt (unbroken domain), Ωft (fracture domain), and Ωtt (transition domain). Then, two phase field threshold values c1 and c2 are used to separate the three domains, if ϕ≥c2 denotes a fracture domain, ϕ≤c1 denotes an unbroken domain, and c1≤ϕ≤c2 denotes a transition domain. In this study, it is assumed that the transition domain can be obtained by interpolation of the fracture and unbroken domains.
(21)χrϕ=1,ϕ≤c1,c2−ϕc2−c1,c1<ϕ<c2,0,ϕ≥c2,χfϕ=0,ϕ≤c1,ϕ−c1c2−c1,c1<ϕ<c21,ϕ≥c2.,In this study, the flow field in the pore elastic medium is expressed in terms of Darcy’s law. The mass conservation equation for the entire domain is written as(22)ρS∂Pf∂t+∇·ρv=qm−ραχr∂εvol∂t,where ρ,S,qm,v,α denote the flow density, storage coefficient, source term, flow velocity, and Biot coefficient, respectively. εvol=∇·u is the volumetric strain.Biot’s coefficientα is expressed as
(23)α=αrχr+αfχf,where αr and αf are the Biot’s coefficients of the reservoir domain and fracture domain, respectively.The storage coefficientS is indicated as follows:
(24)S=εpc+α−εp1−αKVr,where εp,c,KVr are the porosity, fluid compressibility, and bulk modulus, respectively.Darcy’s law flow rate can be expressed as(25)v=−kμ∇Pf,where μ is the fluid viscosity and k is the effective permeability of the transition domain, which is expressed by
(26)k=krχr+kfχf.According to Equation (25), the fluid flow control equation (Equation (22)) for the fluid domain can be rewritten by using fluid pressure as
(27)ρS∂Pf∂t−∇·ρKμf∇Pf=qm−ραχr∂εvol∂t.


### 2.4. Governing Equations of Fluid Distribution

The propagation of multiple clusters of hydraulic fractures is driven by incompressible non-Newtonian fluids. The fluid flow in wellbores consists of two parts: fluid flowing along wellbores and fluid flowing within each fracture. Due to the effect of stress shadowing as well as natural fractures, the injected fluid is distributed nonuniformly within each fracture. For simultaneous fracturing, we assume that the hydrostatic pressure in the horizontal wellbore is constant, and the fluid flowing through the horizontal well into each fracture follows a similar pattern to that of current and voltage. Thus, Kirchhoff’s first and second laws are utilized to calculate the distribution of injected fluid dynamically [44, 45]. The fluid flow is determined by the wellbore friction, perforation friction, and fluid pressure within each fracture. In this paper, using the plane strain assumption to describe fracture propagation and fluid flow will indeed introduce some error into the results, while considering the large kinetic energy of injected fracturing fluid and the limited influence of gravity, it can be approximately simplified to the same flow state in all cross sections, which is the assumption adopted in the classical hydraulic fracturing KGD model [46]. This assumption is widely used in numerical models of hydraulic fracturing [21, 25, 47, 48]. Figure 2 shows a schematic representation of fluid flow in wellbores and fractures. Neglecting wellbore storage effects, the total injection rate is equal to the sum of the injection rates of all fractures:
(28)QT=∑i=1NQi.Figure 2
Illustration of fracturing fluid distribution and pressure drop.Kirchhoff’s second law describes the continuity of pressure along the wellbore, where the pressure in the wellbore heel is equal to the sum of the wellbore friction, perforation friction, and inlet pressure of the fracture [45]:
(29)P0=pw,i+ppf,i+∑Ij=1pcf,j,where pw,i is the inlet pressure of the fracture, ppf,i is the perforation friction, and pcf,i is the wellbore friction. ppf,i is proportional to the square of the flow rate, which is expressed by [49]
(30)ppf,i=0.807249ρnp,i2dp,i4C2Qi2,where Qi is the injection rate (m3/s); ρ is the fracturing fluid density; C denotes the perforation discharge coefficient which ranges from 0.56 to 0.89, it is set to be 0.56 in this work; np,i is the perforation number; and dp,i is the perforation diameter (m). The wellbore friction pressure drop can be calculated by [50]
(31)pcf,i=Ccf∑j=1ixj−xj−1Qwjn′,Qwj=QT−∑k=12j−1Qk,Qwj=QTatj=1,Ccf=23n′+2π−n′k′1+3n′n′n′D−3n′+1,where n′ is the fluid power law index, k′ is the fluid consistency index, D is the horizontal wellbore diameter, and xj is the distance from the heel point of the well to the jth fracture. For the calculation of each item in the pressure drop equation,ppf,i and pcf,i can be calculated by Equation (30) and Equation (31), respectively. The calculation of the fluid pressure pw,i in the fracture requires a large number of computer calculations in order to obtain the variation of fluid pressure in each fracture during the propagation of multiple fractures.


### 2.5. Segregated Scheme

The coupled equations listed above are solved using a segregated scheme that all the unknown variables are solved independently in an increment. The implicit generalized-α method [42] is used to integrate over time to ensure the unconditional stability of the calculation. The solving procedure consists of four separate steps as shown in Figure 3. The displacement tensor and pore pressure u,pij=0 are solved simultaneously in an increment. And the fluid distribution, the history strain, and the phase field variable are solved subsequently. The Newton-Raphson method is utilized to solve each segregated step from the initial guess of u,pij=0, Hij=0, ϕij=0, and Qij=0. At an iteration step j+1 within a time step i, the displacement tensor and pressure u,pij+1 are solved firstly based on the results of the previous step j (u,pij, Hij, ϕij, and Qij). Then, the Hij+1 is updated based on the calculated displacement tensor uij+1. Finally, a new phase filed ϕij+1 and injection distribution Qij+1 are achieved using the solutions of uij+1, pij+1, and Hij+1. The total relative error εr of the solutions can be calculated. The convergence condition requires the εr to be less than a user-defined critical valueε0(εr<ε0). The critical tolerance used in this work is ε0=1×10−4. As the condition is satisfied, a new time increment is implemented. Thus, all the solutions in the range of 0 to the total injection time can be obtained.Figure 3
The calculation process of the model.


## 2.1. Phase Field Approximation of Cracks

[35] proposed that the sharp crack topology can be approximated as a diffusion crack topology, described by using the field variable ϕ, as shown in Figure 1. ϕ=1 means the domain is completely broken, and ϕ=0 means the domain is complete. The model considers a bounded domain Ω, with an outer boundary ∂Ω and an inner boundary Γ. The crack surface density function can be expressed as
(1)γϕ,∇ϕ=ϕ22l0+l02∂ϕ∂xi∂ϕ∂xi,where l0 is a length-scale parameter that limits the width of the transition domain between a completely fractured medium and an intact medium, and its value will affect the accuracy of the results, which needs to be determined by combining experience and sensitivity analysis for different problems. [35] pointed out that the value of l0 should be less than or equal to half of the mesh size in the fracture zone. For the hydraulic fracturing problems, we refer to relevant literature to determine the value of this parameter [34, 36, 37].Figure 1
Schematic diagram of cracks. (a) Sharp crack. (b) Phase field discrete crack.In fracture variational theory, for a pure single-phasic solid, the total potential energy can be represented by elastic strain energy and crack surface energy and external force work. However, for a porous solid, the influence of the fluid pressure of the solid also needs to be considered, and the total potential energy can be expressed as(2)Ψu,Γ=∫ΩψεεdΩ−∫Ωαp·∇·udΩ+∫ΓGcdS−∫Ωb·udΩ−∫∂Ωhif·udS,where ψεε is the elastic strain energy density, Gc is the critical energy release rate, p is the fluid pressure, α is the Biot coefficient, b is the body force, and f is the boundary traction imposed on the outer boundary ∂Ω.In a solid bodyΩ with discrete internal discontinuity Γ, the surface energy of a crack can be represented by a domain integral with the crack phase field:
(3)∫ΓGcdS=∫ΩGcϕ22l0+l02∂ϕ∂xi∂ϕ∂xidΩ.According to [38], it is assumed that the phase field affects only the stretched part of the elastic energy density. Therefore, the elastic energy density can be written as
(4)ψεε=1−k1−ϕ2+kψε+ε+ψε−ε,where 0<k≪1 is the stiffness coefficient, which prevents the stretched part of the elastic energy density from being 0 and avoids the singularity of the numerical value when the phase field ϕ tends to 1, and ψ+ε and ψ−ε are the tensile and compression parts of elastic energy, which can be calculated by
(5)ε=ε++ε−,ε±=∑a=1dεa±na⊗na,ψ±ε=λ2trε±2+μtrε±2,where ε+ and ε− are the tensile and compressive strain tensors, εa and na are the main strains and their directions, the operators •+ and •− are defined as •+=•+•/2 and •−=•−•/2, and λ and μ are the Lamé constants which can be calculated by
(6)λ=Ev1+v1−2v,μ=E21+v.According to Biot’s poroelasticity theory [39], the relationship between the total Cauchy stress σij and the effective stress σije considering the initial in situ stress σ0ij can be expressed as
(7)σij=σije+σ0ij−αPf.By introducing a phase field, the total Cauchy stress across the domain can be expressed as
(8)σij=σije+1−k1−ϕ2+kσ0ij−αPf,inΩ.Combining Equation (4) and Equation (5), the effective stress can be calculated by
(9)σije=∂ψεε∂ε=1−k1−ϕ2+kλtrε+I+2με++λtrε−I+2με−,where I is the identity tensor.The displacement field equilibrium equation considering body forces and boundary conditions can be written as(10)σijnj=fti,on∂Ωt.Finally, the elasticity matrix in the numerical implementation is represented as(11)De=D1111D1122D1133D1112D1123D1113D2211D2222D2233D2212D2223D2213D3311D3322D3333D3312D3323D3313D1211D1222D1233D1212D1223D1213D2311D2322D2333D2312D2323D2313D1311D1322D1333D1312D1323D1313,with Dijkl=D¯ijkl+D~ijkl.The componentD¯ijkl is written as
(12)D¯ijkl=λ1−k1−ϕ2+kHεtrε+Hε−trεδijδkl,where Hεx is a Heaviside function; at x>0, Hεx=1, and at x<0, Hεx=0; δij and δkl denote Kronecker deltas.TheD~ijkl component can be expressed as [40]
(13)D~ijkl=2μ1−k1−ϕ2+kPijkl++Pijkl−,with
(14)Pijkl+=∑a=13∑b=13Hεεaδabnainajnbknbl+∑a=13∑b≠a312εa+−εb+εa−εbnainbjnaknbl+nbknal,Pijkl−=∑a=13∑b=13Hε−εaδabnainajnbknbl+∑a=13∑b≠a312εa−−εb−εa−εbnainbjnaknbl+nbknal,where nai denotes the ith component of vector na.


## 2.2. Phase Field Evolution Equation

According to Equation (3), the energy equation of Equation (2) can be expressed as
(15)Ψu,Γ=∫ΩψεεdΩ−∫Ωαp·∇·udΩ+∫ΩGcϕ22l0+l02∂ϕ∂xi∂ϕ∂xidΩ−∫Ωb·udΩ−∫∂Ωhif·udS.When the first derivative of the energy equationΨu,Γ is equal to zero, Equation (15) takes the minimum value, which means that the slit Γx,t has cracked and extended within time t∈0,Tx∈Ω. Therefore, the governing equation of the phase field can be obtained as
(16)∂σijpor∂xj+bi=0,2l01−kψε+Gc+1ϕ−l02∇2ϕ=2l01−kψε+Gc,where σijpor is the Cauchy stress tensor. Fracture propagation in the phase field method is driven by the tensile part of the elastic energy [41].In order to prevent crack healing when the elastomer is under pressure or unloading, it is necessary to supplement the irreversible condition, and the irreversibility of the crack propagation evolution is ensured by introducing the historical state variableHx,t.
(17)Hx,t=maxs∈0,tψε+εx,s.Combined with Equation (17), the evolution equation in Equation (16) is rewritten as
(18)2l01−kHGc+1ϕ−l02∇2ϕ=2l01−kHGc.The phase field is guaranteed to increase monotonically asH increases. However, this equation does not consider the effect of the initial ground stress. In this research, a new historical reference field was used H as follows [34]:
(19)Hx,t=maxs∈0,tψε+εx,s+σ0:εx,s,inΩ×0,T.The initial hydraulic fracture in this study is considered to be a preexisting notches, and the initial historical field is used to create an initial phase field of [42]
(20)H0x=BGc2l01−2dx,ll0,dx,l≤l02,0,dx,l>l02,where B is a scalar parameter and dx,l is the shortest distance from any point of the model to the discrete fracture.


## 2.3. Fluid Flow in Porous Media

To calculate the flow equations in the domainΩ, the entire domain is divided into three parts [43]: Ωrt (unbroken domain), Ωft (fracture domain), and Ωtt (transition domain). Then, two phase field threshold values c1 and c2 are used to separate the three domains, if ϕ≥c2 denotes a fracture domain, ϕ≤c1 denotes an unbroken domain, and c1≤ϕ≤c2 denotes a transition domain. In this study, it is assumed that the transition domain can be obtained by interpolation of the fracture and unbroken domains.
(21)χrϕ=1,ϕ≤c1,c2−ϕc2−c1,c1<ϕ<c2,0,ϕ≥c2,χfϕ=0,ϕ≤c1,ϕ−c1c2−c1,c1<ϕ<c21,ϕ≥c2.,In this study, the flow field in the pore elastic medium is expressed in terms of Darcy’s law. The mass conservation equation for the entire domain is written as(22)ρS∂Pf∂t+∇·ρv=qm−ραχr∂εvol∂t,where ρ,S,qm,v,α denote the flow density, storage coefficient, source term, flow velocity, and Biot coefficient, respectively. εvol=∇·u is the volumetric strain.Biot’s coefficientα is expressed as
(23)α=αrχr+αfχf,where αr and αf are the Biot’s coefficients of the reservoir domain and fracture domain, respectively.The storage coefficientS is indicated as follows:
(24)S=εpc+α−εp1−αKVr,where εp,c,KVr are the porosity, fluid compressibility, and bulk modulus, respectively.Darcy’s law flow rate can be expressed as(25)v=−kμ∇Pf,where μ is the fluid viscosity and k is the effective permeability of the transition domain, which is expressed by
(26)k=krχr+kfχf.According to Equation (25), the fluid flow control equation (Equation (22)) for the fluid domain can be rewritten by using fluid pressure as
(27)ρS∂Pf∂t−∇·ρKμf∇Pf=qm−ραχr∂εvol∂t.


## 2.4. Governing Equations of Fluid Distribution

The propagation of multiple clusters of hydraulic fractures is driven by incompressible non-Newtonian fluids. The fluid flow in wellbores consists of two parts: fluid flowing along wellbores and fluid flowing within each fracture. Due to the effect of stress shadowing as well as natural fractures, the injected fluid is distributed nonuniformly within each fracture. For simultaneous fracturing, we assume that the hydrostatic pressure in the horizontal wellbore is constant, and the fluid flowing through the horizontal well into each fracture follows a similar pattern to that of current and voltage. Thus, Kirchhoff’s first and second laws are utilized to calculate the distribution of injected fluid dynamically [44, 45]. The fluid flow is determined by the wellbore friction, perforation friction, and fluid pressure within each fracture. In this paper, using the plane strain assumption to describe fracture propagation and fluid flow will indeed introduce some error into the results, while considering the large kinetic energy of injected fracturing fluid and the limited influence of gravity, it can be approximately simplified to the same flow state in all cross sections, which is the assumption adopted in the classical hydraulic fracturing KGD model [46]. This assumption is widely used in numerical models of hydraulic fracturing [21, 25, 47, 48]. Figure 2 shows a schematic representation of fluid flow in wellbores and fractures. Neglecting wellbore storage effects, the total injection rate is equal to the sum of the injection rates of all fractures:
(28)QT=∑i=1NQi.Figure 2
Illustration of fracturing fluid distribution and pressure drop.Kirchhoff’s second law describes the continuity of pressure along the wellbore, where the pressure in the wellbore heel is equal to the sum of the wellbore friction, perforation friction, and inlet pressure of the fracture [45]:
(29)P0=pw,i+ppf,i+∑Ij=1pcf,j,where pw,i is the inlet pressure of the fracture, ppf,i is the perforation friction, and pcf,i is the wellbore friction. ppf,i is proportional to the square of the flow rate, which is expressed by [49]
(30)ppf,i=0.807249ρnp,i2dp,i4C2Qi2,where Qi is the injection rate (m3/s); ρ is the fracturing fluid density; C denotes the perforation discharge coefficient which ranges from 0.56 to 0.89, it is set to be 0.56 in this work; np,i is the perforation number; and dp,i is the perforation diameter (m). The wellbore friction pressure drop can be calculated by [50]
(31)pcf,i=Ccf∑j=1ixj−xj−1Qwjn′,Qwj=QT−∑k=12j−1Qk,Qwj=QTatj=1,Ccf=23n′+2π−n′k′1+3n′n′n′D−3n′+1,where n′ is the fluid power law index, k′ is the fluid consistency index, D is the horizontal wellbore diameter, and xj is the distance from the heel point of the well to the jth fracture. For the calculation of each item in the pressure drop equation,ppf,i and pcf,i can be calculated by Equation (30) and Equation (31), respectively. The calculation of the fluid pressure pw,i in the fracture requires a large number of computer calculations in order to obtain the variation of fluid pressure in each fracture during the propagation of multiple fractures.


## 2.5. Segregated Scheme

The coupled equations listed above are solved using a segregated scheme that all the unknown variables are solved independently in an increment. The implicit generalized-α method [42] is used to integrate over time to ensure the unconditional stability of the calculation. The solving procedure consists of four separate steps as shown in Figure 3. The displacement tensor and pore pressure u,pij=0 are solved simultaneously in an increment. And the fluid distribution, the history strain, and the phase field variable are solved subsequently. The Newton-Raphson method is utilized to solve each segregated step from the initial guess of u,pij=0, Hij=0, ϕij=0, and Qij=0. At an iteration step j+1 within a time step i, the displacement tensor and pressure u,pij+1 are solved firstly based on the results of the previous step j (u,pij, Hij, ϕij, and Qij). Then, the Hij+1 is updated based on the calculated displacement tensor uij+1. Finally, a new phase filed ϕij+1 and injection distribution Qij+1 are achieved using the solutions of uij+1, pij+1, and Hij+1. The total relative error εr of the solutions can be calculated. The convergence condition requires the εr to be less than a user-defined critical valueε0(εr<ε0). The critical tolerance used in this work is ε0=1×10−4. As the condition is satisfied, a new time increment is implemented. Thus, all the solutions in the range of 0 to the total injection time can be obtained.Figure 3
The calculation process of the model.


## 3. Model Verification

### 3.1. Evolution of Fracture Length

The first example is to validate the evolution of fracture length. A rectangular geometry with the size of60m×15m is used to model reservoir domain. An initial fracture with the length of 1.2 m is put in the middle of the domain. The injection rate is 10-3 m2/s. In accordance with the analytical solution, the initial in situ stress is set to be zero. The other parameters are given in Table 1. The numerical solution is compared with the analytical solution of the KGD model. The equation for the analytical solution of the fracture length is expressed as [51]
(32)l=2π2/3E′qtK′2/3,E′=E1−ν2,K′=42π1/2KIC,KIC=GcE1−v2.Table 1
Model validation parameters.
Parameter nameSymbolValueUnitDimension of the 2D domainL60mYoung’s modulusE20GPaPoisson’s ratioν0.2—Fracture energy release rateGc1000N/mFluid viscosityμ1×10−3Pa·sMatrix porosityφ0.05—Matrix permeabilityKr1×10−16m2Regularization parameterl00.3mInjection rateQ10−3m2/sAs shown in Figure4, the numerical solution is matched well with the analytical solution. The reason for the difference in the early fracture lengths is because the fracture propagation in the phase field model starts from 0.6 m, while the fracture propagation in the analytical solution starts from 0. Another reason is that the model medium used by the two methods is different. The PFM uses a poroelastic medium, while the analytical solution model uses an impermeable elastic solid.Figure 4
Comparison of numerical and analytical solutions.


### 3.2. Fracture Aperture for Varying Fluid Pressure

The second example is to validate the fracture aperture for varying fluid pressure. The input parameters are presented in Table2. There is an initial fracture of length 2lc in the plane at y=0. Under plane strain conditions, the displacement in the y direction resulting from the fluid pressure is expressed as [52]
(33)u+0,x=2plcEp1−x2lc21/2,where Ep=E/1−ν2.Table 2
Input parameters for validating fracture aperture.
Parameter nameSymbolValueUnitDimension of the 2D domainL40mYoung’s modulusE30GPaPoisson’s ratioν0.3—Fracture energy release rateGc2000N/mFluid viscosityμ1×10−3Pa·sMatrix permeabilityKr1×10−15m2Regularization parameterl03cmConstant flow source termqm100kg/(m3·s)Figure5 shows the displacement normal to fracture surfaces for varying fluid pressure. It is shown that most of the numerical results are consistent with the analytical solution. Larger errors are observed for the case of high fluid pressure. This is due to the effect of fluid leakage and effective stress, which is not considered in the analytical solution.Figure 5
Comparison of numerical and analytical results in displacement along the initial fracture.


## 3.1. Evolution of Fracture Length

The first example is to validate the evolution of fracture length. A rectangular geometry with the size of60m×15m is used to model reservoir domain. An initial fracture with the length of 1.2 m is put in the middle of the domain. The injection rate is 10-3 m2/s. In accordance with the analytical solution, the initial in situ stress is set to be zero. The other parameters are given in Table 1. The numerical solution is compared with the analytical solution of the KGD model. The equation for the analytical solution of the fracture length is expressed as [51]
(32)l=2π2/3E′qtK′2/3,E′=E1−ν2,K′=42π1/2KIC,KIC=GcE1−v2.Table 1
Model validation parameters.
Parameter nameSymbolValueUnitDimension of the 2D domainL60mYoung’s modulusE20GPaPoisson’s ratioν0.2—Fracture energy release rateGc1000N/mFluid viscosityμ1×10−3Pa·sMatrix porosityφ0.05—Matrix permeabilityKr1×10−16m2Regularization parameterl00.3mInjection rateQ10−3m2/sAs shown in Figure4, the numerical solution is matched well with the analytical solution. The reason for the difference in the early fracture lengths is because the fracture propagation in the phase field model starts from 0.6 m, while the fracture propagation in the analytical solution starts from 0. Another reason is that the model medium used by the two methods is different. The PFM uses a poroelastic medium, while the analytical solution model uses an impermeable elastic solid.Figure 4
Comparison of numerical and analytical solutions.


## 3.2. Fracture Aperture for Varying Fluid Pressure

The second example is to validate the fracture aperture for varying fluid pressure. The input parameters are presented in Table2. There is an initial fracture of length 2lc in the plane at y=0. Under plane strain conditions, the displacement in the y direction resulting from the fluid pressure is expressed as [52]
(33)u+0,x=2plcEp1−x2lc21/2,where Ep=E/1−ν2.Table 2
Input parameters for validating fracture aperture.
Parameter nameSymbolValueUnitDimension of the 2D domainL40mYoung’s modulusE30GPaPoisson’s ratioν0.3—Fracture energy release rateGc2000N/mFluid viscosityμ1×10−3Pa·sMatrix permeabilityKr1×10−15m2Regularization parameterl03cmConstant flow source termqm100kg/(m3·s)Figure5 shows the displacement normal to fracture surfaces for varying fluid pressure. It is shown that most of the numerical results are consistent with the analytical solution. Larger errors are observed for the case of high fluid pressure. This is due to the effect of fluid leakage and effective stress, which is not considered in the analytical solution.Figure 5
Comparison of numerical and analytical results in displacement along the initial fracture.


## 4. Numerical Results and Discussion

To study the effect of stress shadowing and natural fractures on fluid distribution and fracture geometry during multicluster hydraulic fracturing, a square domain with dimensions of60m×60m was created as shown in Figure 6. Six initial perforations of length 3 m were set on either side of the horizontal wellbore and injected at a total injection rate of 9.6 m2/min constant into the horizontal wellbore. In the setting of random natural fractures in the model, we consider bonded natural fractures, whose initial aperture is zero, but the fracture energy is not zero. The initial conditions of the model include in situ stress, pore pressure, and displacement field. Based on the plane strain model, we mainly focus on the horizontal stress and neglect the influence of vertical stress. Where the initial stress is σx=20MPa and σy=10MPa, the initial pore pressure and displacement field are zero. In this paper, the basic parameters of the numerical model adopt the parameters in [34, 53] research, as shown in Table 3.Figure 6
Schematic diagram of multicluster fracturing.Table 3
Numerical model parameters.
Parameter nameSymbolValueUnitDimension of the 2D domainL60mYoung’s modulusE60GPaPoisson’s ratioν0.3—Regularization parameterl00.5mThreshold constantsc10.4—Threshold constantsc21—Fluid densityρ1×103kg/m3Matrix porosityφ0.05—Biot’s modulus of matrixαr0.05—Fluid viscosityμ1×10−3Pa·sMatrix permeabilityKr1×10−15m2Injection rateQ9.6m2/minFluid compressibilityc1×10−81/PaPerforation numbersn15—Perforation diameterd6mmWellbore diameterD0.11mPerforation discharge coefficientC0.56—Power law indexn′0.7—Fluid consistency indexk′0.96—


### 4.1. Hydraulic Fracture Propagation in the Reservoir Containing Different Orientations of Natural Fractures

To study the effect of natural fracture orientations on fluid distribution and hydraulic fracture propagation, five cases with different orientations of natural fracturesθ=30°, θ=45°, θ=60°, θ=75°, and θ=90° were simulated. The input parameters of the model are listed in Table 3.Figure7 shows the fracture network morphology for varying orientations of natural fractures. It can be observed that as the θ is less than 60°, the hydraulic fracture is more likely to be arrested by the natural fractures and escape from an end of the natural fractures and go into the rock matrix along the direction of the maximum principal stress. As the θ greater than 60°, the hydraulic fractures mainly across the natural fractures. In addition, it is shown that the fracture lengths and fluid distribution are more uniform as the θ=60° compared with other cases (Figures 7 and 8(a)). The reason is that when θ is less than 60°, the hydraulic fractures propagate directly along the natural fractures in a horizontal direction. In the case of θ greater than 60°, the hydraulic fracture will propagate directly cross the natural fractures. In these cases, the natural fracture has no significant effect on the propagation of the hydraulic fractures. Therefore, the natural fractures do not have a great influence on the fluid distribution within the fractures, and the fluid distribution in each fracture is mainly determined by the stress shadowing. Figure 8(b) shows the hydraulic fracture propagation length variation in reservoirs containing different orientations of natural fractures. It can be found that when θ is equal to 60°, the multiple hydraulic fractures extend uniformly, but the fracture length is the shortest.Figure 7
Fracture network morphology at different natural fracture orientations. (a)θ=30°, (b) θ=45°, (c) θ=60°, (d) θ=75°, and (e) θ=90°.
(a)(b)(c)(d)(e)Figure 8
(a) Fluid distribution in fractures for varying orientations of natural fractures. (b) Hydraulic fracture length for varying orientations of natural fractures.
(a)(b)Figure9 shows the number of events that hydraulic fractures are crossing and arrested by natural fractures with different orientations. As the θ increases, the number of crossing events increases, while the arresting events decreases. This is because as the θ is low enough, the opening of natural fractures are controlled by the mode-II (shear failure), and the hydraulic fractures are likely to be arrested.Figure 9
Numbers of events that hydraulic fractures are crossing and arrested by natural fractures with natural fracture orientations.


### 4.2. Influence of Horizontal Stress Differences on Hydraulic Fracture Propagation

In this paper, horizontal stress differences (ΔσH) of 5 MPa, 10 MPa, and 15 MPa are considered to investigate the effect of stress difference on fluid distribution and fracture propagation. Figure 10 shows the fracture morphology with varying horizontal stress differences. The results show that at ΔσH=5MPa, the hydraulic fractures propagate deviating from the direction of the maximum horizontal principal stress and then intersecting with each other. In addition, each wing of the hydraulic fractures propagates asymmetrically on both sides, resulting in a complex fracture morphology. As the ΔσH increased to 15 MPa, the deviation of hydraulic fractures away from the direction of the horizontal direction is weakened compared with the other cases. And the hydraulic fracture tends to cross the natural fractures without a significant change in propagation direction.Figure 10
Fracture morphology in the formations with varying horizontal stress differences. (a)ΔσH=5MPa, (b) ΔσH=10MPa, and (c) ΔσH=15MPa.
(a)(b)(c)Figure11(a) shows the number of events that hydraulic fractures are crossing and arrested by natural fractures with different in situ stress differences. It is shown that the amount of crossing events increases with the in situ stress difference. Figure 11(b) shows the fracture propagation length for varying stress differences. At the same injection time, the total propagation length of the fracture increases with the increase of stress difference. The reason is that as the stress difference increases, the degree of fracture deflection gradually decreases and the fracture tends to extend along the direction of the maximum horizontal principal stress, resulting in an increase in the extension length of the fracture. Figure 12(a) demonstrates the fluid distribution in each fracture for varying horizontal stress differences. It is noticed that the stress differences do not have a significant effect on the fluid distribution within the fractures. Figure 12(b) shows the variations of fluid pressure within the fracture for varying horizontal stress differences. The fluid pressure within the fracture decreases with the horizontal stress difference. This is because under low horizontal stress difference conditions, hydraulic fractures are likely be arrested by natural fractures, whereas to extend along natural fractures would require higher fluid pressures.Figure 11
(a) Numbers of events that hydraulic fractures are crossing and arrested by natural fractures with different horizontal stress differences. (b) Fracture propagation length for varying stress differences.
(a)(b)Figure 12
(a) Fluid distribution in fractures and fluid pressure (b) for varying in situ stress differences.
(a)(b)


### 4.3. Effect of Fracture Energy of Natural Fractures on Hydraulic Fracture Propagation

The fracture energy ratio (η=Gc2/Gc1) defined by fracture energy of the rock matrix Gc1 and the natural fracture Gc2 has a significant effect on the morphology of hydraulic fractures. Three values of η, 0.08, 0.16, and 0.32, are considered in the simulations. Figure 13 shows the fracture network morphology for different η. The result shows that as the η equals to 0.08, implying a weak bonding strength of natural fractures, the hydraulic fractures are mainly arrested by the natural fractures. As the η increased to 0.32, the hydraulic fractures extended across the natural fractures, and the trajectories of the hydraulic fracture are more planar. It is implying that the propagation resistance of the interior fracture increases with the η.Figure 13
Fracture network morphology for varying fracture energy ratios. (a)η=0.08, (b) η=0.16, and (c) η=0.32.
(a)(b)(c)Figure14(a) shows the number of events that hydraulic fractures are crossing and arrested by natural fractures with different η. The results show that the number of arresting events decreases with the η, and the crossing events increases with the η. Figure 14(b) shows the fracture propagation length for varying fracture energy ratios. As the ratio of natural fracture to matrix fracture energy increases, the propagation length of the fracture increases. The reason is that with the increase of the fracture energy ratio, the hydraulic fracture is easier to cross the natural fracture and extend along the straight direction, so the length of the fracture is longer. Figure 15 demonstrates the fluid distribution in each fracture for different η. More fluid is allocated in the interior fracture for the case of low η. As the η increases, the influence of the natural fractures on the fluid distribution decreases; as a consequence, the flow distribution is dominated by the stress shadowing.Figure 14
(a) Numbers of events that hydraulic fractures are crossing and arrested by natural fractures with different fracture energy ratios. (b) Fracture propagation length for varying fracture energy ratios.
(a)(b)Figure 15
Fluid distribution in fractures for varying fracture energy.


### 4.4. Effect of Injection Rate on Fracture Propagation

The injection rate is an important controllable parameter in practical engineering. In this paper, three cases of injection ratesQ0, 8 m2/min, 9.6 m2/min, and 11.6m2/min, with the same injected volume were used. Figure 16 shows the fracture morphology for varying injection rates. The results show that at the case of Q0=8m2/min, the propagation of interior fracture is restricted because the allocated fluid in the interior fracture is less than the exterior fractures. As the Q0 increases, the allocated fluid in interior fracture increases, resulting in uniform lengths of the hydraulic fractures. Figure 17 shows the number of events that hydraulic fractures are crossing and arrested by natural fractures with different Q0. The amount of the events increases with the Q0.Figure 16
Fracture network morphology for varying injection rates. (a)Q0=8m2/min, (b) Q0=9.6m2/min, and (c) Q0=11.6m2/min.
(a)(b)(c)Figure 17
Numbers of events that hydraulic fractures are crossing and arrested by natural fractures with different injection rates.Figure18(a) shows the variation of the fluid distribution in the fractures for varying injection rates. The results show that at the case of low Q0, the fluid allocated to the interior fracture is less than exterior fractures. As the Q0 increases, the fluid distribution in each fracture is more uniform. The fluid distribution is controlled by the governing equations Equation (26) and Equation (27). The perforation friction is proportional to the square of the injection rate. Therefore, at the case of low Q0, the perforation friction is significantly less than the other cases, and the fluid distribution is dominated by the stress shadowing which restricts the propagation of interior fracture. Figure 18(b) shows the variations of the fluid pressure at the injection points for varying injection rates. The fluid pressure increases with the injection rate. Moreover, at the same injection volume, the injection time decreases with increasing injection rate.Figure 18
(a) Fluid distribution in fractures and fluid pressure (b) for varying injection rates.
(a)(b)


### 4.5. Effect of Fluid Viscosities on Fracture Propagation

Another important factor affecting the hydraulic fracture propagation morphology is the fracturing fluid viscosity. In order to obtain the effect of different fracturing fluid viscosities on fracture propagation morphology, three cases of fracturing fluid viscosities of 1 MPa·s, 5 MPa·s, and 10 MPa·s are simulated.Figure19 shows the fracture propagation morphology for varying fracturing fluid viscosities. When the fracturing fluid viscosity is equal to 1 MPa·s, the two wings of a fracture show an asymmetric geometry. As the fracturing fluid viscosity reached 10 MPa·s, the hydraulic fractures propagate uniformly. Figure 20(a) shows the variation of the flow distribution within each fracture for varying fracturing fluid viscosities. For the case of low fluid viscosity (μ=1MPa·s), the asymmetric propagation of the fractures results the fluid allocated to the upper fractures (Frac5 and Frac6) lesser than the middle fractures (Frac3 and Frac4). As the fluid viscosity increases, the distribution of the fluid in each fracture becomes more uniform.Figure 19
Fracture network morphology for varying fluid viscosities. (a)μ=1MPa·s, (b) μ=5MPa·s, and (c) μ=10MPa·s.
(a)(b)(c)Figure 20
(a) Fluid distribution in fractures and fluid pressure (b) for varying fluid viscosities.
(a)(b)Figure20(b) demonstrates the fluid pressure within the fracture for varying fluid viscosities. The fluid pressure increases with the fluid viscosity gradually. This is because the fluid leakoff volume into rock matrix decreases with the fluid viscosity, resulting in an increase of the fluid pressure. Figure 21 shows the number of events that hydraulic fractures are crossing and arrested by natural fractures with different fluid viscosities. The number of crossing events increases with the fluid viscosity.Figure 21
Numbers of events that hydraulic fractures are crossing and arrested by natural fractures for different fluid viscosities.


### 4.6. Sensitivity Effects of Natural Fracture Densities on Fluid Distribution as well as Fracture Network

Reservoir heterogeneity introduced by the natural fractures plays a significant role on the hydraulic fracture propagation [54]. In this section, we investigate the influence of natural fracture densities on multifracture propagation. The natural fracture density is defined as the number of natural fractures per area of the formation. Three cases of natural fracture densities (ρn), ρn=7.2×10−3number/m2, ρn=1.9×10−2number/m2, and ρn=3.9×10−2number/m2, were simulated.Figure22 shows the hydraulic fracture propagation morphology for varying natural fracture densities. The results show that when the ρn=7.2×10−3number/m2, the hydraulic fractures mainly propagate along the horizontal direction. As the natural fracture densities increase, the branching and coalescing are observed which facilitate inducing fracture network. Figure 23 demonstrates the variations of the fluid distribution for varying natural fracture densities. Under the condition of a lower ρn, the allocated fluid to the interior fracture is less than the exterior fractures. As the ρn increases, the effect of natural fractures on fluid distribution increases, which facilitates the propagation of the interior fracture.Figure 22
Fracture network morphology for varying natural fracture densities. (a)ρn=7.2×10−3number/m2, (b) ρn=1.9×10−2number/m2, and (c) ρn=3.9×10−2number/m2.
(a)(b)(c)Figure 23
Fluid distribution in fractures for varying natural fracture densities.


## 4.1. Hydraulic Fracture Propagation in the Reservoir Containing Different Orientations of Natural Fractures

To study the effect of natural fracture orientations on fluid distribution and hydraulic fracture propagation, five cases with different orientations of natural fracturesθ=30°, θ=45°, θ=60°, θ=75°, and θ=90° were simulated. The input parameters of the model are listed in Table 3.Figure7 shows the fracture network morphology for varying orientations of natural fractures. It can be observed that as the θ is less than 60°, the hydraulic fracture is more likely to be arrested by the natural fractures and escape from an end of the natural fractures and go into the rock matrix along the direction of the maximum principal stress. As the θ greater than 60°, the hydraulic fractures mainly across the natural fractures. In addition, it is shown that the fracture lengths and fluid distribution are more uniform as the θ=60° compared with other cases (Figures 7 and 8(a)). The reason is that when θ is less than 60°, the hydraulic fractures propagate directly along the natural fractures in a horizontal direction. In the case of θ greater than 60°, the hydraulic fracture will propagate directly cross the natural fractures. In these cases, the natural fracture has no significant effect on the propagation of the hydraulic fractures. Therefore, the natural fractures do not have a great influence on the fluid distribution within the fractures, and the fluid distribution in each fracture is mainly determined by the stress shadowing. Figure 8(b) shows the hydraulic fracture propagation length variation in reservoirs containing different orientations of natural fractures. It can be found that when θ is equal to 60°, the multiple hydraulic fractures extend uniformly, but the fracture length is the shortest.Figure 7
Fracture network morphology at different natural fracture orientations. (a)θ=30°, (b) θ=45°, (c) θ=60°, (d) θ=75°, and (e) θ=90°.
(a)(b)(c)(d)(e)Figure 8
(a) Fluid distribution in fractures for varying orientations of natural fractures. (b) Hydraulic fracture length for varying orientations of natural fractures.
(a)(b)Figure9 shows the number of events that hydraulic fractures are crossing and arrested by natural fractures with different orientations. As the θ increases, the number of crossing events increases, while the arresting events decreases. This is because as the θ is low enough, the opening of natural fractures are controlled by the mode-II (shear failure), and the hydraulic fractures are likely to be arrested.Figure 9
Numbers of events that hydraulic fractures are crossing and arrested by natural fractures with natural fracture orientations.


## 4.2. Influence of Horizontal Stress Differences on Hydraulic Fracture Propagation

In this paper, horizontal stress differences (ΔσH) of 5 MPa, 10 MPa, and 15 MPa are considered to investigate the effect of stress difference on fluid distribution and fracture propagation. Figure 10 shows the fracture morphology with varying horizontal stress differences. The results show that at ΔσH=5MPa, the hydraulic fractures propagate deviating from the direction of the maximum horizontal principal stress and then intersecting with each other. In addition, each wing of the hydraulic fractures propagates asymmetrically on both sides, resulting in a complex fracture morphology. As the ΔσH increased to 15 MPa, the deviation of hydraulic fractures away from the direction of the horizontal direction is weakened compared with the other cases. And the hydraulic fracture tends to cross the natural fractures without a significant change in propagation direction.Figure 10
Fracture morphology in the formations with varying horizontal stress differences. (a)ΔσH=5MPa, (b) ΔσH=10MPa, and (c) ΔσH=15MPa.
(a)(b)(c)Figure11(a) shows the number of events that hydraulic fractures are crossing and arrested by natural fractures with different in situ stress differences. It is shown that the amount of crossing events increases with the in situ stress difference. Figure 11(b) shows the fracture propagation length for varying stress differences. At the same injection time, the total propagation length of the fracture increases with the increase of stress difference. The reason is that as the stress difference increases, the degree of fracture deflection gradually decreases and the fracture tends to extend along the direction of the maximum horizontal principal stress, resulting in an increase in the extension length of the fracture. Figure 12(a) demonstrates the fluid distribution in each fracture for varying horizontal stress differences. It is noticed that the stress differences do not have a significant effect on the fluid distribution within the fractures. Figure 12(b) shows the variations of fluid pressure within the fracture for varying horizontal stress differences. The fluid pressure within the fracture decreases with the horizontal stress difference. This is because under low horizontal stress difference conditions, hydraulic fractures are likely be arrested by natural fractures, whereas to extend along natural fractures would require higher fluid pressures.Figure 11
(a) Numbers of events that hydraulic fractures are crossing and arrested by natural fractures with different horizontal stress differences. (b) Fracture propagation length for varying stress differences.
(a)(b)Figure 12
(a) Fluid distribution in fractures and fluid pressure (b) for varying in situ stress differences.
(a)(b)


## 4.3. Effect of Fracture Energy of Natural Fractures on Hydraulic Fracture Propagation

The fracture energy ratio (η=Gc2/Gc1) defined by fracture energy of the rock matrix Gc1 and the natural fracture Gc2 has a significant effect on the morphology of hydraulic fractures. Three values of η, 0.08, 0.16, and 0.32, are considered in the simulations. Figure 13 shows the fracture network morphology for different η. The result shows that as the η equals to 0.08, implying a weak bonding strength of natural fractures, the hydraulic fractures are mainly arrested by the natural fractures. As the η increased to 0.32, the hydraulic fractures extended across the natural fractures, and the trajectories of the hydraulic fracture are more planar. It is implying that the propagation resistance of the interior fracture increases with the η.Figure 13
Fracture network morphology for varying fracture energy ratios. (a)η=0.08, (b) η=0.16, and (c) η=0.32.
(a)(b)(c)Figure14(a) shows the number of events that hydraulic fractures are crossing and arrested by natural fractures with different η. The results show that the number of arresting events decreases with the η, and the crossing events increases with the η. Figure 14(b) shows the fracture propagation length for varying fracture energy ratios. As the ratio of natural fracture to matrix fracture energy increases, the propagation length of the fracture increases. The reason is that with the increase of the fracture energy ratio, the hydraulic fracture is easier to cross the natural fracture and extend along the straight direction, so the length of the fracture is longer. Figure 15 demonstrates the fluid distribution in each fracture for different η. More fluid is allocated in the interior fracture for the case of low η. As the η increases, the influence of the natural fractures on the fluid distribution decreases; as a consequence, the flow distribution is dominated by the stress shadowing.Figure 14
(a) Numbers of events that hydraulic fractures are crossing and arrested by natural fractures with different fracture energy ratios. (b) Fracture propagation length for varying fracture energy ratios.
(a)(b)Figure 15
Fluid distribution in fractures for varying fracture energy.


## 4.4. Effect of Injection Rate on Fracture Propagation

The injection rate is an important controllable parameter in practical engineering. In this paper, three cases of injection ratesQ0, 8 m2/min, 9.6 m2/min, and 11.6m2/min, with the same injected volume were used. Figure 16 shows the fracture morphology for varying injection rates. The results show that at the case of Q0=8m2/min, the propagation of interior fracture is restricted because the allocated fluid in the interior fracture is less than the exterior fractures. As the Q0 increases, the allocated fluid in interior fracture increases, resulting in uniform lengths of the hydraulic fractures. Figure 17 shows the number of events that hydraulic fractures are crossing and arrested by natural fractures with different Q0. The amount of the events increases with the Q0.Figure 16
Fracture network morphology for varying injection rates. (a)Q0=8m2/min, (b) Q0=9.6m2/min, and (c) Q0=11.6m2/min.
(a)(b)(c)Figure 17
Numbers of events that hydraulic fractures are crossing and arrested by natural fractures with different injection rates.Figure18(a) shows the variation of the fluid distribution in the fractures for varying injection rates. The results show that at the case of low Q0, the fluid allocated to the interior fracture is less than exterior fractures. As the Q0 increases, the fluid distribution in each fracture is more uniform. The fluid distribution is controlled by the governing equations Equation (26) and Equation (27). The perforation friction is proportional to the square of the injection rate. Therefore, at the case of low Q0, the perforation friction is significantly less than the other cases, and the fluid distribution is dominated by the stress shadowing which restricts the propagation of interior fracture. Figure 18(b) shows the variations of the fluid pressure at the injection points for varying injection rates. The fluid pressure increases with the injection rate. Moreover, at the same injection volume, the injection time decreases with increasing injection rate.Figure 18
(a) Fluid distribution in fractures and fluid pressure (b) for varying injection rates.
(a)(b)


## 4.5. Effect of Fluid Viscosities on Fracture Propagation

Another important factor affecting the hydraulic fracture propagation morphology is the fracturing fluid viscosity. In order to obtain the effect of different fracturing fluid viscosities on fracture propagation morphology, three cases of fracturing fluid viscosities of 1 MPa·s, 5 MPa·s, and 10 MPa·s are simulated.Figure19 shows the fracture propagation morphology for varying fracturing fluid viscosities. When the fracturing fluid viscosity is equal to 1 MPa·s, the two wings of a fracture show an asymmetric geometry. As the fracturing fluid viscosity reached 10 MPa·s, the hydraulic fractures propagate uniformly. Figure 20(a) shows the variation of the flow distribution within each fracture for varying fracturing fluid viscosities. For the case of low fluid viscosity (μ=1MPa·s), the asymmetric propagation of the fractures results the fluid allocated to the upper fractures (Frac5 and Frac6) lesser than the middle fractures (Frac3 and Frac4). As the fluid viscosity increases, the distribution of the fluid in each fracture becomes more uniform.Figure 19
Fracture network morphology for varying fluid viscosities. (a)μ=1MPa·s, (b) μ=5MPa·s, and (c) μ=10MPa·s.
(a)(b)(c)Figure 20
(a) Fluid distribution in fractures and fluid pressure (b) for varying fluid viscosities.
(a)(b)Figure20(b) demonstrates the fluid pressure within the fracture for varying fluid viscosities. The fluid pressure increases with the fluid viscosity gradually. This is because the fluid leakoff volume into rock matrix decreases with the fluid viscosity, resulting in an increase of the fluid pressure. Figure 21 shows the number of events that hydraulic fractures are crossing and arrested by natural fractures with different fluid viscosities. The number of crossing events increases with the fluid viscosity.Figure 21
Numbers of events that hydraulic fractures are crossing and arrested by natural fractures for different fluid viscosities.


## 4.6. Sensitivity Effects of Natural Fracture Densities on Fluid Distribution as well as Fracture Network

Reservoir heterogeneity introduced by the natural fractures plays a significant role on the hydraulic fracture propagation [54]. In this section, we investigate the influence of natural fracture densities on multifracture propagation. The natural fracture density is defined as the number of natural fractures per area of the formation. Three cases of natural fracture densities (ρn), ρn=7.2×10−3number/m2, ρn=1.9×10−2number/m2, and ρn=3.9×10−2number/m2, were simulated.Figure22 shows the hydraulic fracture propagation morphology for varying natural fracture densities. The results show that when the ρn=7.2×10−3number/m2, the hydraulic fractures mainly propagate along the horizontal direction. As the natural fracture densities increase, the branching and coalescing are observed which facilitate inducing fracture network. Figure 23 demonstrates the variations of the fluid distribution for varying natural fracture densities. Under the condition of a lower ρn, the allocated fluid to the interior fracture is less than the exterior fractures. As the ρn increases, the effect of natural fractures on fluid distribution increases, which facilitates the propagation of the interior fracture.Figure 22
Fracture network morphology for varying natural fracture densities. (a)ρn=7.2×10−3number/m2, (b) ρn=1.9×10−2number/m2, and (c) ρn=3.9×10−2number/m2.
(a)(b)(c)Figure 23
Fluid distribution in fractures for varying natural fracture densities.


## 5. Conclusions

In this paper, a simultaneous fracturing model is built based on phase field method. The injected fluid into each fracture is dynamically calculated considering perforation friction and wellbore friction. The effects of reservoir properties and controllable parameters on fluid distribution and hydraulic fracture propagation are presented. Hydraulic fracture trajectories mainly depend on the fluid distribution, the stress shadowing, and the interactions with natural fractures. Nonuniform distribution of injected fluid causes uneven growth of hydraulic fractures, especially for the case of low injection rate and fluid viscosity. The natural fractures not only alter the extending directions of the hydraulic fractures but also affect the fluid distribution in each fracture. High densities of natural fractures with low energy ratios could weaken the effect of stress shadowing on fracture propagation and promote the growth of interior fractures.The presented model provides an effective tool to simulate complex behaviors of multiple fracture propagation in naturally fractured reservoirs during stimulation in practical engineering. Although the above results only present three hydraulic fractures, the model can simulate arbitrary amount of hydraulic fractures and clusters without increasing modelling difficulties. In addition, nonlinear constitutive models such as plasticity and damage of rock matrix as well as reservoir heterogeneity can be integrated in the model seamlessly, which needs further work.



---

*Source: 4846474-2023-03-04.xml*