Abstract
Activating the natural fracture (NF) during hydraulic fracture (HF) propagation is a dominant factor for forming the complex fracture network in shale reservoirs. On basis of Renshaw and Pollard’s criterion for fracture intersection, we simulate the NF stability under dynamic disturbance of HF propagation in shale reservoirs with the cohesive zone model in ABAQUS and analyze the effects of various factors on the NF stability. The results show that activation of the NFs in shale hydraulic fracturing is controlled by NF length, NF density, intersection angle between the NF and the HF, horizontal principal stress difference, frac fluid viscosity, and injection rate. During HF propagation, the NF shorter than 0.2 m is less likely to be activated, and those longer than 0.5 m are more likely to be activated. As the NF density increases, not all NFs are activated. The NF at an approach angle larger than 10° is more likely to be activated. The NF is more likely to be activated under the low stress difference, and the NF cannot be activated when the stress difference is higher than 10 MPa. Both frac fluid viscosity and injection rate have effects on activation of NFs, and the NF is more likely to be activated by higher viscosity frac fluid and under higher injection rate. The results provide reference for fracture network stimulation in shale reservoirs.
1. Introduction
There is a great potential for shale gas development in China. China’s shale gas production exceeded in 2019. Compared with the sandstone reservoirs, shale gas reservoirs have extremely low matrix permeability and well-developed NFs. The stimulation technology in shale gas reservoirs is very different from that in sandstone reservoirs. The industrial development of shale gas reservoirs can only be realized through fracture network stimulation. Wide development of NFs in shale reservoirs is a prerequisite for generating the complex fracture network. In the original in situ stress field, the NFs are in equilibrium and cannot propagate spontaneously. Variation of the stress field around the NFs possibly causes rock failure and NF propagation. Access of frac fluid into the NFs causes NF failure and propagation as following conditions. (i) The NF is paralleled to the maximum horizontal principal stress, and its tensile failure occurs. (ii) The NF intersects the maximum horizontal principal stress with a certain angle, and both tensile and shear failure occur. (iii) The NF is not connected to the HF, and its shear failure occurs most likely under the action of the induced stress caused by the HF. Thus, whether a NF can be activated and connected with the HF and other NFs under disturbance of the HF determines complexity and distribution of the fracture network in the shale. The key to activate the NF lies in the mechanism of interaction between the NF and the HF, which is not fully understand currently.
A lot of efforts have been put in theoretical analysis, laboratory experiments, and numerical simulation on interaction between the NF and the HF. Theoretically, the WT criterion [1], the R-P criterion [2], and the GW criterion [3] were proposed to predict the behavior of NF failure. Warpinski and Teufel [1] proposed a criterion for interaction between the NF and the HF, and they determined NF shear and slippage, HF propagation along the NF, and NF tensile failure which causes serious fluid leakoff. Thiercelin and Makkhyu [4] proposed that the NF is reopened before the HF approaches, and they used the dislocation theory to establish a semianalytical model of reopening of NFs approaching the HF. They determined the instability center of NF reopening by obtaining the maximum tensile stress on the NF. Dahi and Olson [5] analyzed the mechanism of HF propagation with the finite element method, and they suggested that fracture initiation occurs from the NF tip when the HF approaches the NF. Renshaw and Pollard [2] established the criterion for intersection between the HF and the NF and proposed the condition of HF crossing NF. They proposed that HF interaction with the NF occurs in the form of crossing NF, opening NF, and arresting HF.
Lamont and Jessen [6] carried out physical simulation in the small rock to study the effect of NFs on HF propagation, and the HF crossed all NFs. Blanton [7] carried out experiments with the naturally fractured core sample in a large triaxial stress equipment and quantified the effects of the NF approach angle and the horizontal principal stress difference on the HF pattern. He found that NF opening, HF arrestion, and crossing NF occur when the HF approaches the NF. Beugelsdijk et al. [8] studied the effects of NFs on the HF propagation under different injection rate, frac fluid viscosity, and horizontal principal stress difference and found that the NF development degree has an obvious effect on the HF pattern. Chen et al. [9] carried out experiments to study the effect of micro- and macro-NF features on the HF propagation direction, and they established a new criterion for NF shearing, opening, and dis-stabilization. They illustrated the mechanism of interaction between the HF and the NF in the naturally fractured reservoirs. Due to the high cost and the difficulty in laboratory experiment on HF propagation, numerical simulation on the effect of NFs on the HF propagation behavior was carried out.
Generally, interaction between the HF and the NF is simulated with the finite element method. Gonzalez-Chavez et al. [10] studied intersection between HF and cohesive NF and analyzed the effects of approach angle, fluid injection rate, etc. on intersection of HF and NF. Dahi-Taleghani and Olson [11] established the numerical model for simultaneous fracture propagation. They found that occurrence and distribution of NFs are key factors, and the NF is possibly opened due to concentration of induced stress at the fracture tip before the HF approaches the NF. Maxwell and Cipolla [12] found that NF opening and slippage occur when the HF approaches the NF in microseismic monitoring. Warpinski et al. [13] proposed that the NF is likely to shear when the HF interferes with the NF. Daneshy [14] believed that small size NFs have little effect on HF propagation, and the medium and large size NFs have a large impact on HF propagation.
Some efforts also have been put in study on the NF stability under disturbance of the HF. Haifeng et al. [15] applied the principle of dynamic fracture propagation to study the mechanism of activating the NF and generating the complex fracture network in shale. Zhao et al. [16] applied the theory of elasticity and fracture mechanics to establish the model of stress field on the NF surface when the HF approaches the NF by considering the closure degree of the NF, and they obtained the analytical solution to the model. The HF affects the NF stability significantly when it approaches the NF. Considering the induced stress of the HF, they analyzed the sensitivity of horizontal principal stress difference, approach angle, approach distance, net pressure, NF closure degree, and HF length. The horizontal principal stress difference, approach angle, and net pressure have most significant effects, and the approach angle and approach distance affect the instability center of NF. Chuprakov et al. [17, 18] studied the physical and mechanical mechanisms of activating the NF when the HF approaches, intersects, and infiltrates the NF on basis of strain energy density theory and displacement discontinuity method. Zhong et al. [19] used the displacement discontinuity theory to analyze the characteristics of NF failure when the HF approaches the NF. On basis of the theory of linear elastic fracture mechanics, Hang et al. [20] used the PKN model to establish a model of the induced stress field when the HF approaches the NF by considering the fluid pressure distribution within the HF and the NF roughness. They studied the features of NF instability and failure as the HF approaches the NF. The NF stability is affected by the horizontal principal stress difference, pore pressure, frac fluid viscosity and fluid injection rate, approach distance, and NF spatial distribution and roughness. When the HF approaches the NF, both tensile and shear failure of NF occur possibly.
The previous study focuses on interaction between a HF and a single NF. The effect of the HF disturbance on the NF stability in the reservoirs with complex NF distribution has rarely been reported. When the HF interacts with the complex NF system, the local-induced stress field occurs around the NF and changes with time. The mechanism of activating the NF in the reservoir with the NF system is more complicated.
The paper is organized as follows. First, we create the model of complex NF with the cohesive unit in ABAQUS with Renshaw and Pollard’s criterion for fracture intersection and assigned the mechanical strength to the NFs. Second, we carry out numerical simulation on the NF stability during HF propagation. Third, we analyze the effects of various factors on the NF stability.
2. Physical Model
The physical model of hydraulic fracturing is established (Figure 1), with the horizontal maximum principal stress in the vertical direction, the horizontal minimum principal stress in the horizontal direction, the injection point in the middle of the bottom boundary, and the angle between the HF and the NF. It is assumed that the fluid within the fracture is an incompressible Newtonian fluid.

To increase the computation efficiency, achieve solution convergence, and improve computation accuracy in the HF region, the model is meshed by the unstructured quadrilateral grid with CPE4P bilinear quadrilateral elements for matrix and COH2D4P elements for fracture (Figure 2). The fracture initiation area and the surrounding area are meshed with relatively dense grids. Fracture propagation along a set path is realized by embedding cohesive units between cell grids in batches, and it is affected by the NF strength. The NF and HF distributions in the simulation model are set (Figure 3).


3. Mathematical Model
3.1. Criterion for Intersection between the NF and the HF
The NF is subject to the action of and in the far field. The angle between the fracture and is 0°-90°. The stress on the NF surface is decomposed into normal stress, shear stress, and fluid pressure within the fracture. Whether the NF is activated and how to activate the NF depends on interaction of these three forces (Figure 4). According to Renshaw and Pollard’s criterion for fracture intersection, the NF is regarded as a frictional interface, and the condition for no slippage along the NF is expressed as follows. where ; or ; is the roots of Eqs. (3) and (4), and it is obtained when , and is calculated in the following equations.

3.2. Constitutive Relation in the Cohesive Unit and Fluid Pressure Equation in the Fracture
Traditional linear elastic fracture mechanics is a main method for studying HF propagation in the solid material. When the stress intensity factor (SIF) at the fracture tip reaches the fracture toughness, the HF propagates forward. In linear elastic fracture mechanics, fracture tip stress has a singularity described by the SIF. Computation of the SIF at the fracture tip is time consuming in some configurations. To avoid these problems, HF propagation is characterized with the cohesive force model.
The cohesion model is an alternative method in analyzing HF propagation. Here, it is assumed that there is a process zone at the fracture tip, and fracture within the process zone is described by the traction-separation (T-S) constitutive model. This avoids the stress singularity at the fracture tip. In addition, the cohesion model is programmed as the cohesion element in the finite element method.
There are several types of T-S constitutive models, and the classic bilinear T-S model is applied here. Fracture initiation in the cohesive element is judged by the quadratic nominal stress rule expressed as follows (Figure 5).

(a) Tension direction

(b) Shear direction
The HF propagation pattern is simulated as the process that fluid enters the element and applies the fluid pressure to the fracture wall when the pore pressure cohesion element is damaged completely.
Flow of the in-compressible Newtonian fluid within the fracture is described by Poiseuille’s law as follows (Figure 6).

Fluid leakoff into the formation is described as follows:
4. Numerical Simulation and Results
4.1. Numerical Simulation Plan
Numerical simulation on the effects of NF length, NF density, NF angle, horizontal principal stress difference, injection rate, and frac fluid viscosity on the NF stability during HF propagation is carried out. Simulation is carried out with the NF length of 0.2 m, 0.5 m, 0.7 m, 1 m, 2 m, and 5 m; the NF at angle of 10°, 30°, 45°, 60°, 70°, and 80°; the NF density of 1/m, 2/m, and 3/m; the frac fluid viscosity of 3 mPa·s, 10 mPa·s, 25 mPa·s, and 50 mPa·s; and the injection rate of 0.5 m3/min, 0.75 m3/min 1 m3/min, 1.25 m3/min, and 2 m3/min, respectively. Finally, analysis of the effects of various factors on the NF stability is carried out. The formation parameters in simulation are listed in Table 1.
The NFs have the lower cementation degree and the lower tensile strength than the rock, and the HF propagation path is affected by the NF when the HF intersects with the NF. To explore the effects of NF length, NF density, NF angle, horizontal principal stress difference, injection rate, and frac fluid viscosity on HF propagation, the model is established in ABAQUS with parameters listed in Table 1.
4.2. Results
4.2.1. Effect of the NF Length
Numerical simulation is carried out with the NF angle of 45°, the horizontal principal stress difference of 8 MPa, the 3 mPa·s frac fluid, the injection rate of 1 m3/min, and the NF length of 0.2 m, 0.5 m, 0.7 m, 1 m, 2 m, and 5 m. The effects of the NF length on the NF stability is analyzed.
When the NF length is 0.2 m, the HF bypasses the NF, and the NF is stable. When the NF length increases to 0.5 m, the NF is activated, and the HF propagates along the NF. When the NF length increases to 5 m, and the HF propagates along one tip of the NF (Figure 7). The NFs with the length of 0.5 m and 0.7 m tend to be activated but have little effect on the HF orientation. The NFs longer than 1 m have significant effect on the HF propagation path. The longer NFs are more likely activated by HFs and accordingly have larger impact on the HF propagation path.

(a) NF length of 0.2 m

(b) NF length of 0.5 m

(c) NF length of 0.7 m

(d) NF length of 1 m

(e) NF length of 2 m

(f) NF length of 5 m
4.2.2. Effect of the NF Density
Numerical simulation is carried out with the NF angle of 45°, the NF length of 1 m, the horizontal principal stress difference of 8 MPa, the 3 mPa·s frac fluid, the injection rate of 1 m3/min, and the NF density of 1/m, 2/m, and 3/m (Figure 8).

(a) NF density of 1/m

(b) NF density of 2/m

(c) NF density of 3/m
When the NF density is 1/m, the NF is activated, and the HF is diverted and propagates to the NF. When the NF density is 2/m, the NF is activated, and the HF is diverted twice and propagates along no. 1 and no. 2 NFs. When the NF density is 3/m, only no. 2 NF is activated. The HF crosses no. 1 and no. 3 NFs, and then, it is diverted and propagates along no. 2 NF. As the NF density increases, not all NFs are connected with the HF. An increase in the NF changes the stress field significantly.
4.2.3. Effect of the NF Angle
Numerical simulation is carried out with the NF length of 1 m, the horizontal principal stress difference of 8 MPa, the 3 mPa·s frac fluid, the injection rate of 1 m3/min, and the NF angle of 10°, 30°, 45°, 60°, and 80°. The effect of the NF angle on the NF stability is analyzed (Figure 9).

(a) 10°

(b) 30°

(c) 45°

(d) 60°

(e) 80°
When the NF angle is 10°, the HF only connects with one tip of the NF, and the HF is diverted slightly. As the NF angle increases, the NF is activated, and the HF is arrested by the NF and propagates forward. The NFs at an approach angle of 45° and 60° have a significant effect on the HF propagation path, and the HF is diverted to the NFs and propagates forward in the 0° direction. When the HF encounters the NF at an approach angle of 80°, it propagates along the NF. The NFs at a larger angle are more likely to be activated and have a significant effect on the HF propagation path. The NFs at a large approach angle are conducive to forming the complex fracture network.
4.2.4. Effect of the Stress Difference
Numerical simulation is carried out with the NF angle of 45°, the NF length of 1 m, the 3 mPa·s frac fluid, the injection rate of 1 m3/min, and the horizontal principal stress difference of 8 MPa and 10 MPa. The effect of the stress difference on the NF stability is analyzed (Figure 10).

(a) Stress difference of 8 MPa

(b) Stress difference of 10 MPa
When the horizontal principal stress difference is 8 MPa, the NF is activated, and the HF is diverted and propagates along the NF. As the stress difference increases to 10 MPa, the NF remains stable. As a result, the NFs are more likely to be activated under the small stress difference. Compared with stress difference is 8 MPa, the stress difference is 10 MPa; HF is diverted a certain distance in the horizontal direction, while stress difference is 8 MPa; HF tends to propagates along the original direction.
4.2.5. Effects of Frac Fluid Viscosity
The frac fluid performance is a key factor in forming the complex fracture network. The high quality frac fluid provides excellent capacity of carrying proppants and creating the fracture. The higher frac fluid viscosity provides better carrying ability but brings greater friction and less fluid leakoff into the reservoir. Here, numerical simulation is carried out with the NF angle of 45°, the NF length of 1 m, the horizontal principal stress difference of 8 MPa, the injection rate of 1 m3/min, and the 3 mPa·s, 10 mPa·s, 25 mPa·s, and 50 mPa·s frac fluid. The effect of the frac fluid viscosity on the NF stability is analyzed (Figure 11).

(a) Frac fluid viscosity of 3 mPa·s

(b) Frac fluid viscosity 10 mPa·s

(c) Frac fluid viscosity 25 mPa·s

(d) Frac fluid viscosity 50 mPa·s
When fracturing the rock with the 3 mPa·s and 10 mPa·s fluid, the NF remains stable, and the HF does not connect with the NF. When the frac fluid viscosity increases to 25 mPa·s and 50 mPa·s, the NF is activated, and the HF connects with the NF and is diverted to the NF. When the frac fluid viscosity increases to 50 mPa·s, the HF propagates toward the NFs and is diverted. An appropriate increase of the fluid viscosity helps activate the NF.
4.2.6. Effect of the Injection Rate
An optimal injection rate helps achieve the optimal fracturing effect and save economic costs. Here, numerical simulation is carried out with the NF angle of 45°, the NF length of 1 m, the horizontal principal stress difference of 8 MPa, the 3 mPa·s frac fluid viscosity, the injection rate of 0.5 m3/min, 0.75 m3/min, 1 m3/min, 1.25 m3/min, and 2 m3/min. The effect of the injection rate on the NF stability is analyzed (Figure 12).

When the injection rate is 0.5 m3/min and 0.75 m3/min, the HF bypasses the NF, and the NF is stable. When the rate increases to 1 m3/min, the HF is diverted and propagates along the NFs, and the NFs are activated. As the injection rate increases to 2 m3/min, the HF propagates along the NF. The HF diversion degree under the rate of 2 m3/min is larger compared with that under the injection rate of 1 m3/min and 1.25 m3/min. The higher injection rate is favorable for NF activation and HF diversion.
4.2.7. Integrated Effects
Forming the complex fracture network in the shale reservoir is the result of interaction of the NF length, the NF density, the NF angle, and the injection rate, which all have an effect on the NF stability. Additional simulation is carried out with the 3 mPa·s fluid, the injection rate of 3 m3/min; the NF length and density of 0.2 m and 3/m, 1 m and 2/m, and 5 m and 1/m; and the NF angle of 10°, 45°, and 80°. The integrated effects of multiple factors on the NF stability are analyzed (Figure 13).

The NF with a length of 0.2 m, an angle of 45°, and a density of 3/m is activated under the injection rate of 3 m3/min. Compared with the low rate of 1 m3/min, where the NF remains stable, the high rate helps activate the shorter NF. The injection rate has a larger effect than the NF length. The NF at an angle of 45° and a length of 5 m is activated, and the HF only connects with one tip of the NF. The NF at an angle of 80° and with a length of 5 m is activated, and the HF is arrested by the NF and propagates forward. The angle has a greater effect than the NF length. The NF at an angle of 10° is not activated under all injection rate and NF density and length. The NF at an angle of 45° is activated. The NF angle has the greatest effect on the NF stability. The factors influencing the NF stability are ranked as the NF angle, the injection rate, the NF length, and the NF density.
5. Conclusions
(1)On basis of Renshaw and Pollard’s criterion for fracture intersection, we simulate the NF stability under the dynamic disturbance of HF propagation in shale with the cohesive zone model in ABAQUS and analyze the effects of various factors on the NF stability. Activation of NFs in shale is controlled by NF length, NF density, angle between the NF and the HF, horizontal principal stress difference, frac fluid viscosity, and injection rate(2)During HF propagation, the NFs longer than 0.5 m are more likely to be activated. The NF longer than 5 m is activated, but the HF propagates along one end of the NF. When the NF density is 3/m, only no. 2 NF is activated. As the NF density increases, not all NFs are activated. The NFs at an angle larger than 10° are more likely to be activated. Thus, the high angle of NFs is conducive to forming the complex fracture network. The NF is more likely to be activated under the low stress difference. When the horizontal principal stress difference is larger than 10 MPa, the NF cannot be activated. Both frac fluid viscosity and injection rate have effects on NF activation. The higher frac fluid viscosity and the higher injection rate are conductive to activating the NF(3)During HF propagation, the NF at an angle of 10° is not activated under all injection rate and NF density and length. The NF at an angle of 45° is activated. The factors influencing the NF stability are ranked as the NF angle, the injection rate, and the NF length
Nomenclature
: | Friction coefficient of the NF surface |
: | The cohesion force of the NF surface |
and : | The shear stress and normal stress on the NF surface |
: | The angle between the NF and the HF |
and : | Horizontal maximum and minimum principal stress in the far field |
: | The rock tensile strength, |
: | The roots of Eqs. (3) and (4), and it is obtained when the maximum principal stress , and is calculated in Eqs. (5)–(8) |
: | The type stress intensity factor at the fracture tip |
and : | The polar coordinates at the fracture tip |
, , and : | Components of the stress |
, , and : | The nominal stress in the normal, first, and second tangential directions of the cohesion element, Pa |
, , and : | The corresponding critical nominal stress |
: | The Macaulay bracket expressed as follows |
: | The fluid flow velocity within the fracture, m/s |
: | The fluid pressure within the fracture, Pa |
: | The fracture width, m |
: | The frac fluid viscosity, Pa·s |
and : | The fluid leakoff rate on the upper and lower surfaces of the cohesion element, m/s |
and : | The leakoff coefficient on the upper and lower surfaces, m/(Pa·s) |
: | The pressure in the middle node of the cohesion element, Pa |
and : | The pore pressure on the upper and lower surfaces of the cohesion element, Pa. |
Data Availability
The data used to support the findings of this study are included in the article.
Conflicts of Interest
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Acknowledgments
Xiangyi Yi was supported by the Fund of National Natural Science Foundation of China (Grant no. U20A20265).