Failure Patterns and Mechanisms of Hydraulic Fracture Propagation Behavior in the Presence of Naturally Cemented Fractures

Daobing Wang,Fang Shi,Hao Qin,＊,Dongliang Sun and Bo Yu

1School of Mechanical Engineering,Beijing Key Laboratory of Pipeline Critical Technology and Equipment for Deep Water Oil&Gas Development,Beijing Institute of Petrochemical Technology,Beijing,102617,China

2Jiangsu Key Laboratory of Advanced Manufacturing Technology,Huaiyin Institute of Technology,Huaiyin,223003,China

ABSTRACT In this study,we use the extended finite element method(XFEM)with a consideration of junction enrichment functions to investigate the mechanics of hydraulic fractures related to naturally cemented fractures.In the proposed numerical model,the lubrication equation is adopted to describe the fluid flow within fractures.The fluid-solid coupling systems of the hydraulic fracturing problem are solved using the Newton-Raphson method.The energy release rate criterion is used to determine the cross/arrest behavior between a hydraulic fracture(HF)and a cemented natural fracture(NF).The failure patterns and mechanisms of crack propagation at the intersection of natural fractures are discussed.Simulation results show that after crossing an NF,the failure mode along the cemented NF path may change from the tensile regime to the shear or mixed-mode regime.When an advancing HF kinks back toward the matrix,the failure mode may gradually switch back to the tensile-dominated regime.Key factors,including the length of the upper/lower portion of the cemented NF,horizontal stress anisotropy,and the intersection angle of the crack propagation are investigated in detail.An uncemented or partially cemented NF will form a more complex fracture network than a cemented NF.This study provides insight into the formation mechanism of fracture networks in formations that contain cemented NF.

KEYWORDS Hydraulic fracturing;natural fractures;crack propagation;unconventional reservoirs;mechanical interaction;joints

1 Introduction

Hydraulic fracturing is a common technique for stimulating hydrocarbon production from tight reservoirs,improving waste disposal,and accelerating heat extraction from geothermal reservoirs[1–4].As exploitation technology improves,more energy from unconventional hydrocarbon resources such as shale gas,tight gas,and coal bed methane are expected to be unlocked[5–9].Natural fractures(NFs)can significantly enhance the effective permeability of a formation by connecting the primary hydraulic fracture(HF)to a network of natural fractures;hence,NFs can play a key role in boosting hydrocarbon production[10].Therefore,it is important to develop strategies for reactivating more NFs during hydraulic fracturing treatments[3,6,11–13].

NFs can be classified into two categories:Uncemented or partially cemented joints and faults,or fully cemented joints and faults[2,3,14].Cemented NFs often have a narrow thickness of less than 0.05 mm and are usually filled with calcite cement[1,15].The cement varies in composition and texture,and a potential HF extension path may exist alongside the weakly bonded interface[3,16].Using semicircular bend(SCB)experiments and finite element method(FEM)simulation,Wang et al.[3]found that the interaction behavior between an HF and cemented NFs is controlled by the rock-cement interfacial bond strength and the thickness of the NFs.Furthermore,numerical simulation results have shown that lower cement strength of NFs often leads to a more complex fracture network[14,15].

The interaction between HFs and NFs can have a considerable influence on fracture geometry,as revealed in experimental and field studies[10,17–19].By conducting a laboratory experiment on a Devonian shale and hydro-stone,Blanton[20]found that an HF is likely to cross preexisting fractures only under high differential stresses and at a high approach angle.According to observations made during mine-back experiments,the propagation of an HF can be arrested by NFs under moderate to low stress.Three possible types of intersection modes between HFs and NFs:NF slips under shear stress,HF arrested,and direct crossing or a crossing with an arrest[21].In addition,the shear strength of pre-existing fractures has an obvious impact on the direction of fracture propagation.Through a series of tri-axial fracturing experiments,Zhou et al.[17]observed that with an increase in the shear strength of NFs,the area between the crossed tendency line and dilated tendency line also increases.Based on experimental observations,Gu et al.[22]proposed an extended HF–NF crossing criterion at nonorthogonal angles,which could be used to judge whether HFs cross/divert into NFs under the condition of isotropic and homogeneous rock.The extended criterion accounted for the effect of the approach angle on crossing,and the artificial fracture was more likely to divert along the interface than to cross it when the approach angle was less than 90 degrees.

In terms of the numerical simulation of the interaction between HFs and NFs,different factors such as the approach angle and the cohesive and frictional properties of NFs have been analyzed to capture their mechanical behaviors in the geometry of HFs[1,14,23–25].The interaction between HFs and NFs may lead to arrest,crossing,or offset[7,22].After an HF–NF interaction,the fracture propagation mode changes from the tensile mode to the mixed mode with some shearing[26].Gu et al.[22]found that the fracture intersection is very sensitive to the approach angle,and when the angle of approach is about 90 degrees,an HF is more likely to cross the interface.Chuprakov et al.[18]developed a new OpenT model to investigate an HF contact with a pre-existing discontinuity,and their numerical results showed that injection parameters such as injection rate and fluid viscosity are major factors in the occurrence of crossing.This new model included a dependency on HF pumping characteristics and NF permeability,which had not been considered in other HF–NF interaction models.In addition,the debonding of an NF may occur when an HF is approaching the NF.Using XFEM,Dahi-Taleghani et al.[1,2]came to the conclusion that stress anisotropy may increase the possibility of opening parallel NFs but that it prevents debonded zones from coalescence with the HF.Thus,it may not enhance well performance in this case.Klimenko et al.[16]subsequently modified the above-mentioned XFEM model to consider HF propagation in both toughness-dominated and viscosity-dominated regimes.Recently,many scholars use XFEM model to simulate hydraulic fracturing problems due to its advantage of avoiding using a conforming mesh[27–30].

To the best of our knowledge,previous studies have concentrated mainly on HF–NF intersection behavior.Adopting the perspective of fracture mechanics,this paper focuses on the failure patterns and mechanisms of HF propagation behavior when an HF diverts into a cemented NF.These patterns and mechanisms are not fully understood at present.Using the XFEM technique with a consideration of junction enrichment functions,this study systematically investigates the mechanical properties of HF–NF intersection behavior,such as fracture toughness,maximum principal stress,and failure mode,according to energy release rate criteria.Crack propagation behaviors between cemented NFs and non-cemented or partially cemented NFs in hydraulic fracturing are also compared.This study offers new insight into the formation mechanism of fracture networks in cemented NF formations.

2 Problem Formulation

2.1 Governing Equations

As shown in Fig.1,we consider the 2D homogeneous,isotropic,and linear elastic properties of a formation that includes an HF with the interfaceΓHF and a cemented NF with the interfaceΓNF.The HF is filled with fracturing fluid injected at constant rateQ0.The HF is expected to intersect with the cemented NF after a period of time.The boundaryΓof the domainΩcomprises prescribed displacement boundaryΓuand prescribed traction boundaryΓt,i.e.,Γ=Γu∪Γt.The outward-pointing unit normal vectors of the fracture surfaces of the HF and NF are denoted asnΓHFandnΓNF,respectively.The effect of fluid leakoff on crack propagation is not considered in this study because of the ultra-low matrix permeability of the formation.We assume that fluid flow inside the HF is incompressible and that the propagation of the HF can be considered a quasi-static process in physics.In this model,no gap between the fracture tip and fluid front is assumed because of the low viscosity of fracturing fluids.

Figure 1:Schematics of a domain containing a hydraulic fracture and a cemented natural fracture

The stress equilibrium equation in the domain and associated boundary conditions can be written as[31–36]:

whereσdenotes the Cauchy stress tensor;bdenotes the body force;¯u denotes the prescribed displacement at the boundaryΓu;tdenotes the traction on the prescribed boundaryΓt;pdenotes fluid pressure in the HF;andtNFdenotes the contact traction vector on the NF surfaces.

We assume that rock behavior is linear elastic;therefore,the constitutive equation can be expressed as:

whereεis the strain tensor;andDdenotes the elastic matrix.In a plane strain problem,D=,whereEis the elastic modulus andνis Poisson’s ratio.

2.2 Crack Propagation Criterion

Maximum circumferential stress is used as the criterion to determine the fracture propagation direction during HF propagation[14,23,37].According to the criterion of maximum circumferential stress,the equivalent stress intensify factorKefor each tip is calculated at each time-step to determine whether the artificial fracture propagates at the corresponding tip[18,19].IfKeis greater than fracture toughnessKIC,the HF is propagating[3,38].According to linear elastic fracture mechanics(LEFM),the equivalent stress intensify factorKecan be written as:

whereKIandKIIdenote mode I and mode II fracture intensify factors,respectively.These factors are obtained by the interaction integral method[14,37];αdenotes the fracture propagation angle,which can be calculated in the local polar coordinate system at the crack tip as:

When an HF interacts with a cemented NF,the HF may be diverted toward the NF,as shown in Fig.2.In this study,the energy release rate criterion is used to determine the cross/arrest behavior between an HF and a cemented NF.In the local polar coordinate system,the energy release rate for any directionθwith respect to the fracture direction can be written as[1,37,39–41]

whereGθdenotes the energy release rate;E∗=E/(1−ν2)is Young’s modulus for plane strain;andθdenotes the polar angle with respect to the fracture tip.According to the criterion,the fracture propagates along the direction that leads to the maximum energy release rate.Therefore,the more likely path of two potential paths can be determined by comparing the ratio ofandat the intersection of a closed cemented NF with an HF.Ifis greater thanthe fracture will cross the cemented NF;otherwise,the HF will be arrested by the cemented NF[2].Here,denotes the rock fracture energy anddenotes the fracture energy of diagenetic cements or the fracture energy between diagenetic cements and the host rock,whichever is lower.

Figure 2:Schematics of a hydraulic fracture intersecting a cemented natural fracture

2.3 Fluid Flow within Hydraulic Fractures

Rather than solving the full Navier–Stokes equation for fluid flow inside fractures,we assume the flow inside a fracture to be flow between two parallel plates.This simplifies the governing equations to Poiseuille’s law to find fluid flowqand fluid pressurep,such that

wherewdenotes the fracture width,μdenotes the fluid viscosity,andsdenotes the coordinate along the crack.Under the assumption of incompressible flow conditions within the fracture and no fluid leakoff into the rock matrix,the mass conservation equation can be expressed as:

By substituting Eqs.(8)into(9),the lubrication equation for fracture flow can be obtained as

wheretandkdenote the injection time and fracture permeability,respectively.Fracture permeabilitykcan be mathematically expressed as follows:

The initial and boundary conditions for fluid flow are zero opening at the tip and the initial time.We further assume that the injection rate is kept constant and there is no flux at the fracture tips

Since the above conditions are Neumann’s boundary conditions,not Dirichlet’s boundary conditions,in order to solve the equations,we need to make sure that the fluid pressure satisfies the global mass conservation law to yield a unique solution,such that

whereΓfracdenotes the trajectory of HFs.

2.4 The Discretization Form of the Problem in XFEM

The extended finite element method(XFEM)was originally developed by Moes and Belytschko based on the partition of unity method(PUM),a key advantage of which is that the finite element mesh does not need to be updated to track the crack path in the problem of crack propagation.In XFEM,the solution space is enriched to differential equations with discontinuous functions.Accordingly,as shown in Fig.3,displacement u(x)in the domain can be approximated as[8,20,25,29,30].

whereSalldenotes the complete set of nodes in the mesh;Sfracdenotes the set of nodes whose supports are divided into two parts by the artificial fracture;Stipdenotes the set of nodes whose supports contain the fracture tip;Sjunctiondenotes the set of nodes whose supports contain two intersecting fractures;NuI(x)denotes the standard finite element shape functions of node I;uIdenotes the standard nodal displacement vector;and aI,and cIdenotes the enriched degree of freedom(DOF)vectors.H(x),Fl(x),andJ(x)denote the enrichment shape functions accounting for the displacement jump across the fracture surfaces,the singular displacement field near the fracture tips,and the displacement field around the intersection point of two fractures,respectively,such that:

whereψm(x)andψs(x)denote the signed distance functions of the main fracture and the secondary fracture,respectively.It can be seen thatJ(x)equals 1,−1,or 0 in different sub-domains created by the intersected fractures.

where Nwdenotes the matrix of shape function that transforms the nodal displacement into fracture opening;and U denotes the unknown nodal displacement vector.

Figure 3:Schematics of enriched nodes of two crossing cracks

By substituting the XFEM approximation of displacement and Eq.(17)into the variational form of the stress equilibrium equations,we obtain the corresponding discretization schemes as follows[23,42–44].

Accordingly,the fluid pressure in the hydro-fractures can be approximated as:

whereShfdenotes the set of nodes of the fluid pressure elements along the hydro-fracture;anddenotes the shape function of nodal fluid pressurepIat node I.

Similarly,we can substitute the above-mentioned FEM approximation into the variational forms of lubrication equation and use the forward Euler time discretization to address the time derivative in this equation.The corresponding discretization schemes of lubrication equation can be written as:

where Δtdenotes the time step;Pdenotes the unknown nodal pressure vector;and ΔU denotes the increment of vectorUduring time step Δt.

In Eq.(18),the global stiffness matrixKis defined as follows:

where Dcontdenotes the contact stiffness of fracture interfaces;and Bstdand Benr,respectively,denote the B matrix based on standard shape function and enrichment shape function,as shown in Eq.(13).

Accordingly,the coupling matrixQand the external force vector Fextare respectively defined as

In Eq.(20),flow matrixHand source termSare respectively written as:

2.5 The Fluid–Solid Coupling Strategy and the Related Algorithm

By combining Eqs.(18)with(20),we can obtain the fluid–solid coupling systems of the hydraulic fracturing problem as follows:

This fully coupled equation set is nonlinear since the cubic term of fracture permeabilitykexists in flow matrixHand the frictional interaction between fracture surfaces is also considered in global stiffness matrixK.Therefore,the Newton–Raphson iterative algorithm is utilized to solve the nonlinear problem at each time step.The corresponding residualRnand Jacobian matrixJnat iteration stepnare respectively expressed as:

Thus,the iterative scheme for the fully coupled equations at iteration stepncan be written as:

This iteration process converges only when nodal fracture opening vectorwand nodal fluid pressure vectorPsimultaneously satisfy the following convergence criterion[43,45,46].

where ‖·‖ denotes the L2-norm operator;andanddenote the specified tolerance and take the value of 10−3in this paper.

3 Results and Discussion

3.1 Model Verification

In order to verify our model,we compared the numerical results of 2D hydraulic fracturing based on the XFEM technique with the analytical solutions of the well-known KGD model(Geertsma &De Klerk).HFs can be categorized as toughness- or viscosity-dominated processes according to the dimensionless fracture toughness,which can be written as:

In the case ofKm>4,HFs are toughness-dominated,whereas in the case ofKm<0.5,they are viscosity-dominated[47–50].

In this verification model,the rectangular domain has a length of 100 m and a width of 180 m,and the injection point is located at the center of the edge-width.To reduce the amount of calculation that is necessary,the model is symmetric with respect to the edge-width.Using an injection rate of 0.001 m2/s,fracturing fluid with 1 mPa·s viscosity is injected into the rectangular area for 30 s.The input parameters of the model are listed in Tab.1.According to Eq.(31),Kmis equal to 0.313 for this model,and thus the hydro-fracture belongs to viscositydominated propagation.The 100 m×180 m rectangular domain is divided into 3,080 quadrilateral elements in total.The initial half-length of hydraulic fracture is equal to 1.25 m.Using the aforementioned XFEM technique,the fracture width along the hydro-fracture and fluid pressure changes over time at the injection point can be calculated,as shown in Fig.4.We observe that the XFEM results show very good agreement with the results of the analytical solutions of the KGD model[14,23,42,49,50].

Table 1:Input parameters of the 2D hydraulic fracturing problem

3.2 The Effect of the Length of the Upper/Lower Portion of the Cemented Natural Fracture on Crack Propagation

In the model shown in Fig.2,the constant length of the cemented NF is equal to 6 m under the state of isotropic stress,which is the sum of the length of the upper portion of the NF,Lupper,and the corresponding lower length,Llower.The direction of the maximum horizontal principal stress is along the vertical axis.The input parameters of the model are listed in Tab.2.The black line in Fig.4 shows the hydro-fracture propagation paths when an HF is intersecting with a closed cemented NF.We observe that in all cases,the advancing HF diverts along the lower side of the NF and then kinks back to propagate along the original fracturing direction near the lower end of the NF.Because the ratio ofGfracandGrockis less than unity,the hydro-fracture will grow along the cemented natural fracture according to the HF–NF criterion presented in Section 2.3.The advancing hydro-fracture diverts back to the original fracturing direction,which can be attributed to the combined action of the local crack tip stress field and the far-field stress field.According to Irwin’s relation,the conversion between the energy release rate and fracture toughness can be written as:

whereE∗=E/(1 −ν2)is the plane-strain elastic modulus;Eis the elastic modulus;andνis Poisson’s ratio.

As shown in Fig.5,there are two high tensile stress zones on the contour maps of maximum principle stress.One is near the intersection point of HF–NF;the other is in the neighborhood of crack tips.This shows that the failure mode is tensile-dominated near the two zones.However,the tensile stress value near the crack tips is greater than that near the intersection point,which is the reason that the diverted HF propagates along just one side of the NF.

Table 2:Input parameters of the model for HF–NF interaction

Figure 4:Results of the XFEM technique and of the analytical solutions:(a)Fluid pressure at the injection point(b)Fracture width

The distributions of fracture width and shear slippage along the propagation paths are shown in Fig.6.We observe that there is a saltation for the values of both the fracture opening and shear slippage near the HF–NF intersection point(fracture length=12.5 m).When the hydrofracture propagates from the intersection point to the lower end of the cemented NF,the value of the shear slippage is greater than that of the fracture opening,which indicates that the failure mode in this zone is shear-dominated.When the hydro-fracture kinks back toward the rock matrix from the end of the NF,the value of the shear slippage is close to zero,which shows that the failure mode is tensile-dominated.Therefore,the hydro-fracture propagation will shift from shear fracture to tensile fracture after the HF–NF intersection.We also observe that shear slippage in the zone of the lower portion of the NF decreases as the lower length of the NF,i.e.,Llower,increases,while the variation tendency of the fracture width in the same zone is the opposite.This indicates that an increase in the length of the lower portion of the NF promotes tensile failure and weakens shear failure.

Figure 5:Maximum principal stress at different lengths of the lower and upper parts of a natural fracture:(a) Llower=3 m, Lupper=3 m;(b) Llower=2 m, Lupper=4 m;(c) Llower=5 m, Lupper=1 m;(d) Llower=6 m, Lupper=0 m

The curves of stress intensify factorsKIandKIIof the lower crack tip at each fracturing step are plotted in Fig.7.When the fracturing step moves from 6 to 8(corresponding to the lower portion of the NF),the value ofKIIpeaks while the value ofKIreaches its minimum.After fracturing step 8,KIIis reduced to zero,andKIis much greater thanKII.The variation in stress intensify factors is consistent with that of fracturing width and shear slippage.

The net pressure distribution along the fracture length at different levels of NF length is shown in Fig.8.We observe that the net pressure along the hydro-fracture increases as the length of the lower portion of the NF decreases,except in the case ofLlower=6 m.The reason is that opening the lower portion of NF requires a higher shear stress for the shorter length of the lower portion of the NF in Fig.4.While the length of the lower portion of NF is equal to 6 m,the hydro-fracture just reaches the lower end of the NF at fracturing step 10.Therefore,it requires more energy to make the advancing hydro-fracture kink back toward the rock matrix.

Figure 6:Distribution of fracture width and shear slippage along the hydro-fracture length at different levels of natural fracture length:(a)Fracture width(b)Shear slippage

Figure 7:Distributions of stress intensify factors of the lower crack tip for each fracturing step at different levels of natural fracture length:(a) KI;(b) KII

Figure 8:Net pressure distribution along the fracture length at different levels of natural fracture length

3.3 The Effect of Horizontal Stress Anisotropy on Crack Propagation

In this section,we use the parameters listed in Tab.1 to investigate the effect of horizontal stress anisotropy on the HF–NF interaction.Thex-axis along the HF coincides with the perpendicular bisector of the cemented NF with a length of 6 m.The direction of maximum principal stress is along the verticaly-axis direction,and its crack propagation paths at different levels of horizontal differential principal stress are as shown in Fig.9.These paths correspond to 0,1,3,and 5 MPa.We observe that horizontal stress anisotropy has an obvious impact on the growth of the advancing hydro-fracture.In all cases,the HF will first be arrested by the NF before propagating only along the lower portion of the NF.The HF diverts into the rock matrix when it subsequently arrives at the lower tip of the NF.Under the combined action of the far–field stress and the local crack tip stress field,the stronger the horizontal stress anisotropy is,the more likely the diverted HF is to deviate from the direction of the original HF,i.e.,the horizontalx-axis direction.

Figure 9:Maximum principal stress at different levels of horizontal differential stress:(a) σH=5 MPa, σh=5 MPa;(b) σH=5 MPa, σh=4 MPa;(c) σH=5 MPa, σh=2 MPa;(d) σH=5 MPa,σh=0 MPa

From the contour maps of maximum principal stress shown in Fig.9,we are able to observe that tensile stress zones exist at the intersection point of HF–NF and at the lower tip portion of the NF,which can explain why the advancing HF is arrested by the NF before propagating in a single side direction.With the increase in horizontal stress anisotropy,the tensile stress zone expands below the lower portion of the NF,as shown by the green/yellow color region of Fig.9c.Therefore,the diverted HF grows away from the horizontal direction more easily.

The distribution of fracture width and shear slippage along the hydro-fracture length at different levels of horizontal differential stress is shown in Fig.10.We observe the same result as in Fig.6,i.e.,there is abrupt change in both the fracture opening and shear slippage in the vicinity of the HF–NF intersection point.When the HF enters the lower portion of the NF,the fracture width is approximately equal to 1 mm,while the shear slippage ranges from 2 to 3 mm.Obviously,the failure mode is shear-dominated in this segment of the lower portion of the NF.When the HF is diverted into the rock matrix,the fracture opening is greater than the value of the shear slippage.This indicates that the failure mode of the diverted HF is tensiledominated.Therefore,we conclude that the failure mode shifts from the shear-dominated regime to the tensile-dominated regime after the HF–NF interaction,which is consistent with the result in Fig.6.In addition,under the condition of increased horizontal stress anisotropy,the values of the fracture opening and shear slippage decrease when the fracture length is less than 12.5 m(corresponding to the primary hydro-fracture).However,when the advancing HF is arrested by the NF,the fracture opening is augmented and the shear slippage decreases with the increase in horizontal stress anisotropy.In this study,we demonstrate that strong stress anisotropy weakens the tendency of shear failure while enhancing the tendency of tensile failure in the segment of the lower portion of the NF.

Figure 10:Distribution of fracture width and shear slippage along the hydro-fracture length at different levels of horizontal differential stress:(a)Fracture width;(b)Shear slippage

The curves of stress intensify factorsKIandKIIof the lower crack tip at each fracturing step are shown in Fig.11.From the sixth to the eighth fracturing step,the HF propagates along the lower portion of the NF.We can see that the value ofKIIpeaks at the eighth fracturing step,whileKIreaches its minimum value during the same steps.Afterwards,KIIsharply decreases andKIsharply increases.This indicates that the HF failure mode shifts from the shear regime to the tensile regime,which is consistent with the analytical results of the fracture opening and shear slippage in Fig.10.

Figure 11:Distributions of stress intensify factors of the lower crack tip for each fracturing step at different levels of natural fracture length:(a) KI;(b) KII

The net pressure distribution along the fracture length at different levels of horizontal differential stress is shown in Fig.12.We observe that the net pressure along the hydro-fracture decreases as horizontal stress anisotropy increases,except in the case of Δσ=3 MPa.However,the net pressure increases when the diverted HF kinks back toward the rock matrix.This is because higher tensile stress is needed to split the rock matrix.Therefore,it requires more energy to make the advancing hydro-fracture kink back toward the rock matrix.

Figure 12:Net pressure distribution along the fracture length at different levels of horizontal differential stress

3.4 The Effect of the Approach Angle on Crack Propagation

In this section,the effect of the approach angle on the propagation of an HF is numerically investigated based on input parameters similar to those listed in Tab.1.In the isotropic stress state,we observe that the approach angle has a strong impact on the extension paths of the HF,as shown in Fig.13.For all four cases,i.e.,approach angle=90 degrees,60 degrees,45 degrees,and 30 degrees,the HF is arrested by the cemented NF and can then propagate on only one side along the lower/upper portion of the NF before finally kinking back toward the rock matrix.For the latter three cases,the HF propagates along the upper side of the NF,while it extends along the lower side in the case ofβ=90 degrees.We also observe that with the decrease in the approach angle,the diverted HF is more likely to kink back toward the rock matrix.

The contour maps of maximum principal stress at different approach angles are shown in Fig.13.We observe that there is a tensile stress zone in the vicinity of one side of the crack tips,while the relatively weak tensile zone disappears with a decrease in the approach angle.This indicates that it is most difficult to make the primary HF divert/kink back toward the rock matrix in the case ofβ=90 degrees.

The distribution of fracture width and shear slippage along the hydro-fracture length at different approach angles is shown in Fig.14.We observe that there is an abrupt change in both the fracture opening and shear slippage along the fracture length.The value of this sudden change peaks in the case ofβ=90 degrees.In the cases of decreased approach angles,the values are comparatively small.We verify that the HF is easy to divert in the context of a small approach angle.Forβ=90 degrees,the fracture opening at the interval with a 12.5 to 15 m fracture length is close to zero,and the corresponding shear slippage is about 4 mm.Forβ=60 degrees,45 degrees,and 30 degrees,the fracture opening is more than 2 mm,and the corresponding shear slippage is about 2 mm.This demonstrates the failure mode transition from the shear regime to the tensile regime with the decrease in the approach angle.

The curves of stress intensify factorsKIandKIIof a single-side crack tip at each fracturing step are shown in Fig.15.We can see that the value ofKIIpeaks around the sixth to eighth fracturing step,while the value ofKIreaches its minimum during the same steps.Afterwards,KIIsharply decreases whileKImaintains its relatively high value.This indicates that the HF failure mode shifts from the shear regime to the tensile regime,which agrees well with the result in Fig.14.

The net pressure distribution along the fracture length at different levels of horizontal differential stress is shown in Fig.16.For the approach angleβ=90 degrees,the net pressure at the interval with a 12.5 to 15 m fracture length undergoes a sudden change,while the corresponding sudden change for the non-orthogonal approach angle becomes weak.For the non-orthogonal approach angle,the net pressure experiences a sudden change only when the HF kinks back toward the rock matrix.

3.5 Discussion

The main difference between frictional natural fracture and cemented natural fracture is that natural fracture contains calcite cement[14,15].There is not calcite cement in frictional natural fractures,but the two fracture surfaces are in contact.In order to better understand the mechanical interaction between HFs and cemented NFs,we compared the crack propagation behavior of uncemented NFs and cemented NFs,as shown in Fig.17.In the case of an interaction between an HF and an uncemented NF,the advancing HF will propagate along the two sides of the NF.This is quite different from the case of a cemented NF.The tensile stress zone can be only seen in the vicinity of the crack tips after approaching the frictional NF.This indicates that the failure mode is tensile-dominated throughout the interaction process.A comparison of the fracture opening and shear slippage between uncemented NFs and cemented NFs is shown in Fig.17b.For convenience,the crack propagation path is depicted along the upper side of the NF in both cases.It is obvious that their fracture openings are very close,while the shear slippage of the cemented NF is much greater than that of the uncemented NF.This is very consistent with the results of failure mode.Wang et al.’s simulated numerical results indicated that uncemented NFs attract the propagating crack,causing a higher stress intensity factor compared to that of the cemented NFs.The findings in this study support those of Wang et al.[3].

Figure 13:Maximum principal stress at different approach angles:(a) β=90 degrees;(b) β=60 degrees;(c) β=45 degrees;(d) β=30 degrees

We numerically simulate the crack propagation paths when an HF is approaching two NFs,as shown in Fig.18.For frictional NFs,the HF will extend along the upper and lower directions,thereby forming a complex fracture network.For cemented NFs,the HF can only propagate along one side direction.Therefore,the fracture network of cemented NFs is relatively simple compared to that of uncemented NFs.

Figure 14:Distribution of fracture width and shear slippage along the hydro-fracture length at different approach angles:(a)Fracture width;(b)Shear slippage

Figure 15:Distributions of stress intensify factors of the lower crack tip for each fracturing step at different levels of natural fracture length:(a) KI;(b) KII

Figure 16:Net pressure distribution along the fracture length at different approach angles

Figure 17:Comparison of crack propagation behavior between uncemented NFs and cemented NFs:(a)The maximum principal stress for the interaction between an HF and uncemented NF(b)The fracture opening and shear slippage between an uncemented NF and a cemented NF

Figure 18:Comparison of crack propagation paths between uncemented NFs and cemented NFs:(a)One HF and two uncemented NFs;(b)One HF and two cemented NFs

4 Conclusion

Using the XFEM technique,this paper systematically analyzed the mechanical mechanism of HF behavior in cemented formations.Mechanical factors include maximum principal stress,stress intensify factors,and shear slippage.Key factors in crack propagation,such as the length of the upper/lower portion of the cemented NF,horizontal stress anisotropy,and the approach angle,were investigated in detail.Our preliminary conclusions are as follows:

(1)When an HF encounters a cemented NF,the growing HF is more likely to propagate along only one side of the NF.The corresponding mechanical mechanism lies in that the tensile stress zone is mainly located on one side of the crack tip of the NF,while the value of the stress intensify factor on the other side is approximately equal to zero.After crossing the intersection point,the failure mode along the cemented NF shifters from the tensiledominated regime to the shear-dominated regime.Meanwhile,shear slippage is greater than the fracture opening along the NF path,and the value of the fracturing width is very small.When the advancing HF kinks back toward the rock matrix,the failure mode shifts back to the tensile-dominated regime.

(2)Key factors such as the length of the upper/lower portion of the cemented NF,horizontal stress anisotropy,and the approach angle have a non-negligible effect on crack propagation paths.Under the given conditions,an increase in the length of the lower portion of the cemented NF or horizontal stress anisotropy will promote tensile failure and weaken shear failure,and the net pressure along the HF increases as the length of the lower portion of the NF decreases.The approach angle has a strong impact on the extension paths of the HF.With a decrease in the approach angle,the failure mode shifts from the shear regime to the tensile regime.

(3)When an HF encounters an uncemented NF,the advancing HF will propagate along the two sides of the NF,which is quite different from the case of the cemented NF.The failure mode is mainly tensile-dominated throughout the interaction process.The shear slippage of the cemented NF is much greater than that of the uncemented NF.This is very consistent with the results of failure mode.Uncemented NFs will form more complex fracture networks than cemented NFs.The results in this paper provide new insight into the mechanisms of fracture network generation.Future research should consider the primary and secondary relations of various factors using a combination of experiments and numerical calculations.

Acknowledgement:The authors wish to express their appreciation to the reviewers for their helpful suggestions which greatly improved the presentation of this paper.

Funding Statement:This work was financially supported by the National Science Foundation of China(Grant Nos.51804033 and 51936001),Natural Science Foundation of Jiangsu Province(Grant No.BK20170457),Program of Great Wall Scholar(Grant No.CIT&TCD20180313)and Jointly Projects of Beijing Natural Science Foundation,Beijing Municipal Education Commission(Grant No.KZ201810017023),and Beijing Youth Talent Support Program(CIT &TCD201804037).

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.