# The Propagation Behavior of Hydraulic Fracture in Rock Mass with Cemented Joints

Authors: Yongxiang Zheng; Jianjun Liu; Yun Lei
Journal: Geofluids
 (2019)

---

## Abstract

The formation of the fracture network in shale hydraulic fracturing is the key to the successful development of shale gas. In order to analyze the mechanism of hydraulic fracturing fracture propagation in cemented fractured formations, a numerical simulation about fracture behavior in cemented joints was conducted based firstly on the block discrete element. And the critical pressure of three fracture propagation modes under the intersection of hydraulic fracturing fracture and closed natural fracture is derived, and the parameter analysis is carried out by univariate analysis and the response surface method (RSM). The results show that at a low intersecting angle, hydraulic fractures will turn and move forward at the same time, forming intersecting fractures. At medium angles, the cracks only turn. At high angles, the crack will expand directly forward without turning. In conclusion, low-angle intersecting fractures are more likely to form complex fracture networks, followed by medium-angle intersecting fractures, and high-angle intersecting fractures have more difficulty in forming fracture networks. The research results have important theoretical guiding significance for the hydraulic fracturing design.


---

## Body

## 1. Introduction

The porosity and permeability of shale are extremely low. To achieve the effective development of shale gas, fractures are required to provide flow channels for gas. Therefore, hydraulic fracturing is applied to increase fractures in shale. Natural weakness planes contribute to the formation of a complex hydraulic fracture network. Fortunately, many shale outcrops, cores, and image logs show that there are weakness planes including joints, faults, bed-parallel fractures, early compacted fractures, and fractures associated with concretions in shale reservoir [1–3]. There is also some evidence of natural fractures as shown in Figure 1. Those outcrops from north China provide proofs of cemented fractures and natural fracture networks in rocks (Figures 1(a) and 1(b)), while the cores from the Changning-Weiyuan Shale provide proofs of open or cemented natural fractures in deep formation (Figures 1(c) and 1(d)). Those weakness planes including open fractures and cemented cracks are heterogeneous. They can reduce rock strength and may divert hydraulic fracture propagation. Therefore, the interaction between the weakness planes and hydraulic fractures is significant for the formation of a complex fracture network.Figure 1
Weakness planes in rocks.
(a)
Cemented fractures of outcrops
(b)
Natural fractures network of outcrops
(c)
Natural fractures of cores (depth 3635 m)
(d)
Cemented fractures of cores (depth 1531 m)The experimental study primarily focuses on investigating the behaviors of a hydraulic fracture intersecting with a natural fracture. Hydraulic fracturing tests at a laboratory scale have been carried out at different approximation angles and differential stresses, which were run for a condition where the hydraulic fracture opened existing fractures to determine if more constructive crossing interaction could be obtained. For hydraulic fractures, when the angle of approach angle was 30 degrees, existing fractures tend to open and prevent the induced fracture from crossing. At an angle of 60 and 90, hydraulic fractures tend to cross existing fractures when the differential stress was high enough. Based on the results, an elastic solution has been used as a basis for a hydraulic/natural fracture interaction criterion for fracture opening [4]. But the shear slip of the discontinuous surface was ignored in Blanton’s criterion which only considers the condition that the fractures opened when the fluid pressure overcame the normal stress on the discontinuous surfaces. Therefore, Warpinski and Teufel gave the critical pressure of the shear slip on the fracture surfaces [5]. Then, a series of experimental investigations designed to assess the conditions required for the crossing were presented. And a stress solution of the linear elastic fracture mechanics near a fracture tip was used to calculate the criterion of the hydraulic/natural fracture interaction [6]. Besides, other related experiments about the natural/hydraulic fracture interaction were performed with a certain purpose. A criterion has been developed for fractures crossing frictional interfaces with cohesion at nonorthogonal intersection angles. Laboratory experiments were conducted to validate the criterion [7]. Bahorich et al. performed 9 hydraulic fracturing experiments in gypsum cement blocks. And the cemented natural fractures were represented by plaster discontinuities [8]. And a series of semicircular bend tests on Marcellus Shale core samples containing calcite-filled natural fractures (veins) were performed to investigate the influence of weak planes on hydraulic fracture propagation [9].Additional numerical analyses have also been performed to develop propagation criteria for determining propagation direction when a hydraulic fracture encounters natural fractures. For example, the discrete element method [1, 10, 11], extended finite element method [12–14], cohesive zoned method [1, 15], phase field method [16], and displacement discontinuity method [17] were implemented to perform the numerical simulation which analyzes the crossing or arresting behavior of hydraulic and natural fractures. The discrete element method is one of the effective methods for hydraulic fracturing. Five typical coal models were established to simulate hydraulic fracturing in the coal seam based on the two-dimensional particle flow code (PFC2D) by Wang et al. [10]. Hydraulic fracture propagation and proppant transport in a two-layer formation with stress drop were studied by Zhang and Dontsov [18].Previous studies mainly focus on the influence of damaged natural fractures on hydraulic fractures. The jointed surface of the weakness planes is not cracked due to cementation, but its strength is weaker. When the fracturing fluid reaches the weakness plane, the initiation and propagation of induced fractures are likely along the weakness plane. Therefore, the interaction between the weakness plane and the hydraulic fracture has an important influence on fracture propagation.The purpose of this paper is to comprehensively consider the expansion path and morphology of hydraulic fractures intersecting with natural closed fractures under the shear and tensile action of fracture surfaces. The expansion patterns of hydraulic fractures after they meet with intersecting closed fractures are of great significance for the formation of complex fracture networks after fracturing. In this paper, a numerical simulation about fracture behavior in cemented joints was conducted based firstly on block discrete element. Then, the failure modes of hydraulic fractures in fractured formation are discussed, and the failure conditions under different failure modes are given. Based on the different failure conditions of the intersection, the fracture morphology with two boundaries after fracture failure at different intersecting angles is discussed. Besides, the influence of parameters on fracture morphology was analyzed by univariate analysis and response surface method (RSM).


## 2. Numerical Model

### 2.1. Block Discrete Element Method [18, 19]

The block discrete element method is a three-dimensional numerical modeling for advanced geotechnical analysis. The block discrete element method simulates the response of discontinuous media (such as jointed rock or masonry bricks) that is subject to either static or dynamic loading. The discontinuous material is represented as an assemblage of discrete blocks. Models may contain a mix of rigid or deformable blocks. Deformable blocks are defined by a continuum mesh of finite-difference zones, with each zone behaving according to a prescribed linear or nonlinear stress-strain law. The relative motion of the discontinuities is also governed by linear or nonlinear force-displacement relations for movement in both the normal and shear directions. Joint models and properties can be assigned separately to individual, or sets of, discontinuities.The block discrete element method is based on the distinct element method (DEM) for discontinuum modeling. It differs from particle-based methods in its ability to represent a zero-initial porosity condition, as well as interlocked irregular block shapes that provide resistance to block rotation (moments) after contact breakage [20]. The block discrete element method consists of the continuum and discrete model simulations. Continuum blocks need not be represented by particles. Therefore, the calculation amount of this method is less than that of particle-based methods. This method can deal with larger models with lower equipment requirements than particle-based methods.


#### 2.1.1. The Constitutive Model of Joints

The thickness of the joint is far less than the scale of its plane, so its deformation characteristics are described by the stress-displacement relationship. The constitutive relation of the structural plane mainly studies the relationship between the stress of the structural plane and its normal deformation and tangential deformation. There are two stresses on the structure surface: normal stressσ and shear stress τ, and two corresponding displacements: normal displacement δn and tangent displacement δs. The stress-displacement relation matrix is expressed as
(1)στ=KnKnsKsnKsδnδs,where Kn=∂σ/∂δn, which is the normal stiffness coefficient and represents the effect of normal displacement on normal stress; Ks=∂τ/∂δs, which is the shear stiffness coefficient and represents the effect of shear displacement on shear stress; Kns=∂σ/∂δs, which is the dilation stiffness coefficient and represents the effect of shear displacement on normal stress; and Ksn=∂τ/∂δn, which represents the effect of normal displacement on shear stress.The Goodman element is adopted in this calculation; i.e.,Ksn and Kns are set as 0 and Kn and Ks are respectively used to describe the normal deformation and tangential deformation of the structural surface.The Mohr-Coulomb equation for shear strength of the structure plane is(2)τ=c+fσ,where f=tanφ and c and φ are the cohesive force and friction angle of the structural plane, respectively.


#### 2.1.2. Contact Friction Joint Model

The contact friction joint model and the coulomb slip model are used in this calculation. It is assumed that the normal stress incrementFn (positive pressure) and the shear vector increment Fs between blocks are proportional to the normal displacement un and the tangential displacement increment us, respectively, in the elastic stage. The contact stiffness of the normal and tangential springs is kn and ks, respectively. Assuming that there is no tension in the joint and Coulomb’s law is satisfied, then the relation can be expressed as
(3)Fn=knun,whileun≤0,Fn=0,whileun>0,Fs=ksus,whileFcs≤fFcn+cL,Fs=signu̇sfFcn+cL,whileFcs>fFcn+cL,where Fn and Fs are the normal and tangential components of contact force Fc, respectively; kn and ks are the normal and tangential stiffness of joints, un and us are the normal and tangential relative displacements of joints, respectively; F and c are the friction coefficient and cohesive force of joint material, respectively; and L is the length of the contact surface.


#### 2.1.3. Fluid Flow in Crack [10, 21]

It is assumed that the fracture surface is impermeable and fluid flows only within the fracture surface. The flow of fluid in the crack satisfies the cube law. Using the modified cube law,(4)q=−uh3ρg12μ∇h=−Kh∇h,where q is the fluid flow rate per crack width, uh is the equivalent hydraulic aperture, ρ is the density of the fluid, μ is the fluid dynamic viscosity, ∇h is the hydraulic gradient, and Kh is the hydraulic conductivity.In the hydromechanics coupling of fractures, the influence of mechanical deformation on fracture permeability is mainly manifested as the change of fracture aperture. In the elastic stage, the crack aperture is expressed as an equation related to the effective stress.(5)uhe=uh0+f1Δσn′kn.In the plastic stage, considering the expansion caused by a joint slip, the crack opening can be expressed as(6)uhp=uh0+f1Δσn′kn+f1uhs,where uh0 is initial fracture aperture, Δσn′ is effective stress, and f1 is a coefficient about the influence of stiffness and is set as 1 in this study.


#### 2.1.4. Motion Equation of Blocks

The deformation of blocks is nonnegligible. The block is divided into tetrahedral units. The vertices of a tetrahedral element are called grid difference points. The motion equation is established at each node as follows:(7)u̇i=∫sσijnjds+Fim+gi,where s is the outside surface, m is the mass on the grid point, gi is the acceleration of gravity, and Fi is the resultant force of an external force applied to a node.The node force∑Fi is 0 in the balance state. Otherwise, there is an acceleration of the node according Newton’s second law.
(8)u̇it+Δt/2=u̇it−Δt/2+∑FitΔtm.For every time step, strain and rotation are related with displacement. Their forms are as follows:(9)ε̇ij=12u̇i,j+u̇j,i,θij=12u̇i,j−u̇j,i.The incremental method is applied in calculation. It is not limited to small strain problems. Therefore, the constitution model of deformation blocks takes incremental form. The equation is as follows:(10)Δσije=λΔεvδij+2μΔεij,where λ and μ are Lamé coefficients, and δij is the Kronecker function.


### 2.2. Numerical Model Description

Based on the block discrete element theory, the propagation behavior of hydraulic fractures in rock mass with cemented joints is simulated by the following model. The model including two cutters and one vertical joint that can simulate the behavior of fracture crossing and fracture containment is shown in Figure2. As shown in Figure 2(b), the intersection angles between vertical and two cutters are α1 and α2. The size of the model is 20 m in 3 dimensions. A 3D model of rock with joints and elements is shown in Figure 2(c). The distance between the two cutter centers is 4 m. And the parameters are shown in Table 1.Figure 2
Model for numerical simulation: (a) 3D model of rock with joints; (b) a side view of joints; (c) 3D model with elements.
(a)
(b)
(c)Table 1
The values of parameters.
Parameter
Value
Parameter
Value
Elasticity modulus
20 GPa
Fluid viscosity
0.0015 cp
Poisson’s ratio
0.25
Fluid density
1000 kg/m3
Density of rock
2600 kg/m3
Depth
3000 m
Joint friction
20°
Joint tension
0 MPa
Joint cohesion
0 MPa
Initial vertical stress
σz=hpg=76.518MPaIn this simulation, theα1 and α2 are set with the same angle. 30°, 45°, and 90° are chosen to investigate the effect of the intersection angle on crack propagation. And the vertical joint is divided into 3 parts by two cutters. The fluid-driven crack is always contained in the middle part of the vertical joint. The simulation results are shown in the next section.


## 2.1. Block Discrete Element Method [18, 19]

The block discrete element method is a three-dimensional numerical modeling for advanced geotechnical analysis. The block discrete element method simulates the response of discontinuous media (such as jointed rock or masonry bricks) that is subject to either static or dynamic loading. The discontinuous material is represented as an assemblage of discrete blocks. Models may contain a mix of rigid or deformable blocks. Deformable blocks are defined by a continuum mesh of finite-difference zones, with each zone behaving according to a prescribed linear or nonlinear stress-strain law. The relative motion of the discontinuities is also governed by linear or nonlinear force-displacement relations for movement in both the normal and shear directions. Joint models and properties can be assigned separately to individual, or sets of, discontinuities.The block discrete element method is based on the distinct element method (DEM) for discontinuum modeling. It differs from particle-based methods in its ability to represent a zero-initial porosity condition, as well as interlocked irregular block shapes that provide resistance to block rotation (moments) after contact breakage [20]. The block discrete element method consists of the continuum and discrete model simulations. Continuum blocks need not be represented by particles. Therefore, the calculation amount of this method is less than that of particle-based methods. This method can deal with larger models with lower equipment requirements than particle-based methods.


### 2.1.1. The Constitutive Model of Joints

The thickness of the joint is far less than the scale of its plane, so its deformation characteristics are described by the stress-displacement relationship. The constitutive relation of the structural plane mainly studies the relationship between the stress of the structural plane and its normal deformation and tangential deformation. There are two stresses on the structure surface: normal stressσ and shear stress τ, and two corresponding displacements: normal displacement δn and tangent displacement δs. The stress-displacement relation matrix is expressed as
(1)στ=KnKnsKsnKsδnδs,where Kn=∂σ/∂δn, which is the normal stiffness coefficient and represents the effect of normal displacement on normal stress; Ks=∂τ/∂δs, which is the shear stiffness coefficient and represents the effect of shear displacement on shear stress; Kns=∂σ/∂δs, which is the dilation stiffness coefficient and represents the effect of shear displacement on normal stress; and Ksn=∂τ/∂δn, which represents the effect of normal displacement on shear stress.The Goodman element is adopted in this calculation; i.e.,Ksn and Kns are set as 0 and Kn and Ks are respectively used to describe the normal deformation and tangential deformation of the structural surface.The Mohr-Coulomb equation for shear strength of the structure plane is(2)τ=c+fσ,where f=tanφ and c and φ are the cohesive force and friction angle of the structural plane, respectively.


### 2.1.2. Contact Friction Joint Model

The contact friction joint model and the coulomb slip model are used in this calculation. It is assumed that the normal stress incrementFn (positive pressure) and the shear vector increment Fs between blocks are proportional to the normal displacement un and the tangential displacement increment us, respectively, in the elastic stage. The contact stiffness of the normal and tangential springs is kn and ks, respectively. Assuming that there is no tension in the joint and Coulomb’s law is satisfied, then the relation can be expressed as
(3)Fn=knun,whileun≤0,Fn=0,whileun>0,Fs=ksus,whileFcs≤fFcn+cL,Fs=signu̇sfFcn+cL,whileFcs>fFcn+cL,where Fn and Fs are the normal and tangential components of contact force Fc, respectively; kn and ks are the normal and tangential stiffness of joints, un and us are the normal and tangential relative displacements of joints, respectively; F and c are the friction coefficient and cohesive force of joint material, respectively; and L is the length of the contact surface.


### 2.1.3. Fluid Flow in Crack [10, 21]

It is assumed that the fracture surface is impermeable and fluid flows only within the fracture surface. The flow of fluid in the crack satisfies the cube law. Using the modified cube law,(4)q=−uh3ρg12μ∇h=−Kh∇h,where q is the fluid flow rate per crack width, uh is the equivalent hydraulic aperture, ρ is the density of the fluid, μ is the fluid dynamic viscosity, ∇h is the hydraulic gradient, and Kh is the hydraulic conductivity.In the hydromechanics coupling of fractures, the influence of mechanical deformation on fracture permeability is mainly manifested as the change of fracture aperture. In the elastic stage, the crack aperture is expressed as an equation related to the effective stress.(5)uhe=uh0+f1Δσn′kn.In the plastic stage, considering the expansion caused by a joint slip, the crack opening can be expressed as(6)uhp=uh0+f1Δσn′kn+f1uhs,where uh0 is initial fracture aperture, Δσn′ is effective stress, and f1 is a coefficient about the influence of stiffness and is set as 1 in this study.


### 2.1.4. Motion Equation of Blocks

The deformation of blocks is nonnegligible. The block is divided into tetrahedral units. The vertices of a tetrahedral element are called grid difference points. The motion equation is established at each node as follows:(7)u̇i=∫sσijnjds+Fim+gi,where s is the outside surface, m is the mass on the grid point, gi is the acceleration of gravity, and Fi is the resultant force of an external force applied to a node.The node force∑Fi is 0 in the balance state. Otherwise, there is an acceleration of the node according Newton’s second law.
(8)u̇it+Δt/2=u̇it−Δt/2+∑FitΔtm.For every time step, strain and rotation are related with displacement. Their forms are as follows:(9)ε̇ij=12u̇i,j+u̇j,i,θij=12u̇i,j−u̇j,i.The incremental method is applied in calculation. It is not limited to small strain problems. Therefore, the constitution model of deformation blocks takes incremental form. The equation is as follows:(10)Δσije=λΔεvδij+2μΔεij,where λ and μ are Lamé coefficients, and δij is the Kronecker function.


## 2.1.1. The Constitutive Model of Joints

The thickness of the joint is far less than the scale of its plane, so its deformation characteristics are described by the stress-displacement relationship. The constitutive relation of the structural plane mainly studies the relationship between the stress of the structural plane and its normal deformation and tangential deformation. There are two stresses on the structure surface: normal stressσ and shear stress τ, and two corresponding displacements: normal displacement δn and tangent displacement δs. The stress-displacement relation matrix is expressed as
(1)στ=KnKnsKsnKsδnδs,where Kn=∂σ/∂δn, which is the normal stiffness coefficient and represents the effect of normal displacement on normal stress; Ks=∂τ/∂δs, which is the shear stiffness coefficient and represents the effect of shear displacement on shear stress; Kns=∂σ/∂δs, which is the dilation stiffness coefficient and represents the effect of shear displacement on normal stress; and Ksn=∂τ/∂δn, which represents the effect of normal displacement on shear stress.The Goodman element is adopted in this calculation; i.e.,Ksn and Kns are set as 0 and Kn and Ks are respectively used to describe the normal deformation and tangential deformation of the structural surface.The Mohr-Coulomb equation for shear strength of the structure plane is(2)τ=c+fσ,where f=tanφ and c and φ are the cohesive force and friction angle of the structural plane, respectively.


## 2.1.2. Contact Friction Joint Model

The contact friction joint model and the coulomb slip model are used in this calculation. It is assumed that the normal stress incrementFn (positive pressure) and the shear vector increment Fs between blocks are proportional to the normal displacement un and the tangential displacement increment us, respectively, in the elastic stage. The contact stiffness of the normal and tangential springs is kn and ks, respectively. Assuming that there is no tension in the joint and Coulomb’s law is satisfied, then the relation can be expressed as
(3)Fn=knun,whileun≤0,Fn=0,whileun>0,Fs=ksus,whileFcs≤fFcn+cL,Fs=signu̇sfFcn+cL,whileFcs>fFcn+cL,where Fn and Fs are the normal and tangential components of contact force Fc, respectively; kn and ks are the normal and tangential stiffness of joints, un and us are the normal and tangential relative displacements of joints, respectively; F and c are the friction coefficient and cohesive force of joint material, respectively; and L is the length of the contact surface.


## 2.1.3. Fluid Flow in Crack [10, 21]

It is assumed that the fracture surface is impermeable and fluid flows only within the fracture surface. The flow of fluid in the crack satisfies the cube law. Using the modified cube law,(4)q=−uh3ρg12μ∇h=−Kh∇h,where q is the fluid flow rate per crack width, uh is the equivalent hydraulic aperture, ρ is the density of the fluid, μ is the fluid dynamic viscosity, ∇h is the hydraulic gradient, and Kh is the hydraulic conductivity.In the hydromechanics coupling of fractures, the influence of mechanical deformation on fracture permeability is mainly manifested as the change of fracture aperture. In the elastic stage, the crack aperture is expressed as an equation related to the effective stress.(5)uhe=uh0+f1Δσn′kn.In the plastic stage, considering the expansion caused by a joint slip, the crack opening can be expressed as(6)uhp=uh0+f1Δσn′kn+f1uhs,where uh0 is initial fracture aperture, Δσn′ is effective stress, and f1 is a coefficient about the influence of stiffness and is set as 1 in this study.


## 2.1.4. Motion Equation of Blocks

The deformation of blocks is nonnegligible. The block is divided into tetrahedral units. The vertices of a tetrahedral element are called grid difference points. The motion equation is established at each node as follows:(7)u̇i=∫sσijnjds+Fim+gi,where s is the outside surface, m is the mass on the grid point, gi is the acceleration of gravity, and Fi is the resultant force of an external force applied to a node.The node force∑Fi is 0 in the balance state. Otherwise, there is an acceleration of the node according Newton’s second law.
(8)u̇it+Δt/2=u̇it−Δt/2+∑FitΔtm.For every time step, strain and rotation are related with displacement. Their forms are as follows:(9)ε̇ij=12u̇i,j+u̇j,i,θij=12u̇i,j−u̇j,i.The incremental method is applied in calculation. It is not limited to small strain problems. Therefore, the constitution model of deformation blocks takes incremental form. The equation is as follows:(10)Δσije=λΔεvδij+2μΔεij,where λ and μ are Lamé coefficients, and δij is the Kronecker function.


## 2.2. Numerical Model Description

Based on the block discrete element theory, the propagation behavior of hydraulic fractures in rock mass with cemented joints is simulated by the following model. The model including two cutters and one vertical joint that can simulate the behavior of fracture crossing and fracture containment is shown in Figure2. As shown in Figure 2(b), the intersection angles between vertical and two cutters are α1 and α2. The size of the model is 20 m in 3 dimensions. A 3D model of rock with joints and elements is shown in Figure 2(c). The distance between the two cutter centers is 4 m. And the parameters are shown in Table 1.Figure 2
Model for numerical simulation: (a) 3D model of rock with joints; (b) a side view of joints; (c) 3D model with elements.
(a)
(b)
(c)Table 1
The values of parameters.
Parameter
Value
Parameter
Value
Elasticity modulus
20 GPa
Fluid viscosity
0.0015 cp
Poisson’s ratio
0.25
Fluid density
1000 kg/m3
Density of rock
2600 kg/m3
Depth
3000 m
Joint friction
20°
Joint tension
0 MPa
Joint cohesion
0 MPa
Initial vertical stress
σz=hpg=76.518MPaIn this simulation, theα1 and α2 are set with the same angle. 30°, 45°, and 90° are chosen to investigate the effect of the intersection angle on crack propagation. And the vertical joint is divided into 3 parts by two cutters. The fluid-driven crack is always contained in the middle part of the vertical joint. The simulation results are shown in the next section.


## 3. Results and Discussion

### 3.1. Fracture Propagation in Rocks with Cemented Joints

The time-dependent injection pressure was recorded. The result is shown in Figure3. As shown in the figure, the pressure increases rapidly at the initial period. The crack was initialized when it reached maximum pressure at 56.4 MPa. Then, there is a vibration period after crack open. At last, the pressure tended to become smooth and steady. This is consistent with the trend of actual monitoring data on site, which means that the simulation results are suitable for reference research.Figure 3
The history of injection pressure with time.One of the main purposes of numerical simulation is to investigate the behavior of a fluid-driven crack in the rock with cemented joints. In other words, the shape of open fractures after fracking is one of the most important goals. The shape is described by the aperture of fractures in this study. The results are plotted in Figure4. Figures 4(a)–4(c) is the spatial distribution of fractures with different angles between the vertical joint and two cutters. As shown in the pictures, with different intersection angles, the final shape is varied. According to a predecessor’s studies, it is easy to understand that the crack extends straight forward with a high intersection angle, while the crack is always captured by the meeting natural fracture but would not stretch forward when the intersection angle is large enough. However, at the situation that two angles are 30°, it penetrates again. Besides, it is also captured by natural fractures. Why can the fluid-driven crack penetrate again? This is determined by the stress field around the fractures. In Section 3.2, the fracture propagation criterion based on elastic mechanics is discussed.Figure 4
The aperture of joints: the spatial distribution of fractures with different angles between the vertical joint and two cutters: (a) 30-30, (b) 45-45, and (c) 90-90; the aperture of vertical joints with different angles between the vertical joint and two cutters: (d) 30-30, (e) 45-45, and (f) 90-90.
(a)
(b)
(c)
(d)
(e)
(f)Another noteworthy feature different from the 2D model simulation is that fracture height is restrained by two cutters. In the 2D model, the fluid-driven crack must extend forward in their own form, penetrating or captured, when approaching two joints. However, in the 3D model, it is not obligatory for the fluid-driven crack extending forward as in the 2D. Figure4(d)–4(f) is the aperture of vertical joints with different angles at 30°, 60°, and 90°, respectively. As shown in the figures, no matter what the angles, the vertical propagation is restricted by the two cutters. And it extends left and right. The fracture containment problem is another problem still to be solved.


### 3.2. Critical Pressure for the Failure of the Weakness Plane

There are many weak surface structures, such as bedding and fractures, in shale due to the action of geological tectonic stress. It is more likely to be damaged in fracturing for the reason that the strength of the weakness plane is low. The fracture will be preferentially expanded along the weakness plane. Therefore, the intersection of the hydraulic fracture and the weakness plane is the key problem of fracture propagation in the formation with the weakness plane.The intersection model of the hydraulic fracture and the natural fracture is shown in Figure5. Suppose there are two intersecting closed natural fractures, NF1 and NF2, at an approach angle of α. The assumptions in this paper are as follows. The hydraulic fracture extends along the fracture plane NF2 parallel to the maximum principal stress. The hydraulic fracture tip reaches the intersection of the two fractures but has not yet passed the intersection point. NF1 and NF2 have certain tensile and shear strength before cracking, and the tensile and the shear strength disappear after cracking. To simplify the calculation, this paper selects the same method as the W&T criterion and assumes that the fracture is blunted at the interface, so that the singularity of the hydraulic fracture tip can be ignored. More detailed information can be found in [5].Figure 5
The intersection model of the hydraulic fracture and the natural fracture.The failure modes of NF1 and NF2 under fluid action are different with different angles of fracture and principal stress. In general, the failure form of NF1 may be a tensile failure or a shear failure. After tensile failure, cracks will open directly. However, cracks do not necessarily open after shear failure. It also requires that the fluid pressure can overcome the normal stress on the fracture surface and prop up the fracture. The failure mode of NF2 without shear stress on the fracture surface can only be tensile failure. The failure forms are summarized in Table2. Different failure forms should meet the following conditions.Table 2
The criteria for different failure forms of fractures.
Fracture
Failure forms
Criteria
NF1
Shear failure and opening
P=P∣P>Pshear∪P>Popen
Tensile failure
P=P∣P>Pten1
NF2
Tensile failure
P=P∣P>Pten2P is the fluid pressure at the hydraulic fracture tip. Pshear is the critical pressure to shear failure of NF1. Popen is the opening pressure of NF1 without bond strength. Pten1 is the critical tensile failure pressure of NF1 considering the bond strength. Pten2 is the critical tensile failure pressure of NF2 considering the bond strength.Once the critical conditions are determined, the failure forms of fracture propagation can be determined after obtaining the failure conditions of different forms of fractures. Therefore, the critical pressure of different fracture failure forms will be calculated in the following part of this section.According to Warpinski and Teufel, the critical fluid pressure in the crack during shear failure of NF1 can be calculated.(11)pshear=2τ0−σ1−σ3sin2α+kfcos2β−kf2kf+σ3,where σn is the normal stress on the fracture surface, τ is the shear stress on the fracture surface, σ1 is the maximum principal stress, σ3 is the minimum principal stress, α is the approach angle, τ0 is the inherent shear strength of the interface, and kf is the coefficient of friction.The strength of the weakness plane disappears after cracking. When the fluid pressure inside the crack is greater than the normal stress on the fracture surface, the fracture surface opens. Therefore, the fluid pressure in the crack after shear failure is(12)popen=σn=σ1sin2α+σ3cos2α=σ1−σ3sin2α+σ3.Tensile failure occurs when the fluid pressure in NF1 is greater than the sum of normal stress and tensile strength of the fracture surface. At this point, the critical pressure of fluid in the joint is(13)pten1=σn+σt=σ1−σ3sin2α+σ3+σt,where σt is the tensile strength of the cemented fracture.Tensile failure occurs when the fluid pressure of fracture NF2 is greater than the sum of normal stress on the surface and tensile strength of the fracture. At this time, the critical fluid pressure of tensile failure for fracture NF2 is(14)pten2=σ3+σt.The pressure conditions for crack opening can be determined by the above formulas. When the natural fracture NF1 opens, the hydraulic fracture is captured by fracture NF1. When fracture NF2 opens, the hydraulic fracture passes through fracture NF1. The critical pressure of fracture failure at different angles of NF1 and NF2 is calculated to determine the fracture expansion form at the intersection point.


### 3.3. The Behavior of Hydraulic Fracture When Approaching Intersecting Cemented Joints

The critical pressure of different failure forms of cracks can be calculated by the formulas (11)–(14) and Table 2. The parameters selected in this article are shown in Table 3. And the results are shown in Figure 6.Table 3
The values of parameters.
Parameter
Value
Parameter
Value
Inherent shear strength,τ0
2 MPa
Coefficient of friction,kf
0.2
Maximum principle stress,σ1
20 MPa
Tensile strength of cemented joint,σt
2 MPa
Minimum principle stress,σ3
12 MPaFigure 6
The critical pressure under different failure conditions.Obviously, for NF1, the critical pressure of shear failure decreases first and then increases with the increase of the intersecting angle. The critical pressure of tensile failure and opening after shear failure increase with the increase of the intersecting angle. However, the critical pressure of the NF2 tensile failure is independent of the intersection angle.At a small angle (0° and 10.67°), as the pressure in the crack rises, fluid pressure first reachedPopen. However, NF1 is not broken by the shear at this time, so it cannot be opened. The pressure then reached Pten1 and Pten2. At this point, the two values are close; both NF1 and NF2 undergo tensile failure. Hydraulic fractures extend to both NF1 and NF2, forming cross fractures. At the angle range of 10.67° to 15°, the fracture fluid pressure in the first achieves Popen, and no fractures are failures now. Then, the pressure continued to rise to Pshear and shear failure occurred in NF1. At this point, NF1 shear failure and cracking occurred and hydraulic fracture only extended in NF1. When the fracture intersection angle is 15° to 37.7°, the fluid pressure first reaches Pshear, and then NF1 shear failure occurs. At this point, the pressure has not reached the opening pressure Popen. As the pressure rises to Popen, NF1, which has been shear damaged, cracks. When the intersection angle is 37.7° to 68°, the shear failure of NF1 occurs first, but it does not open, because the fluid pressure at this time is less than Popen. With the increase of pressure, the fluid pressure reaches Pten2 and the tensile failure of NF2 occurs. Since NF1 is not open, hydraulic cracks will only propagate along NF2. When the intersection angle is between 68 and 90 degrees, the fluid pressure first reaches Pten2 and the hydraulic fracture will only extend in NF2. It is noteworthy that the values of Pten2 and Popen are close to each other in the range of 37.7 degrees, at which time NF1 and NF2 will open simultaneously.To sum up, with the change of the intersecting angle, the expansion morphology of intersecting cracks can be divided into three categories. At low angles (0-10.67), NF1 and NF2 open at the same time, forming a cross-state. At the middle angle (10.67-37.7), shear cracking occurred in NF1 and no cracking occurred in NF2. The fracture morphology shows turning. At high angles (37.7-90), NF2 cracks and NF1 does not open. The fracture shape is a single fracture extending forward. In summary, the intersection angle can be divided into three parts according to the fracture propagation pattern, namely, low angle, medium angle, and high angle. The boundary between the low angle and the middle angle is defined as theL-M boundary (10.67 in Figure 6). And the boundary between the medium angle and the high angle is defined as the M-H boundary (37.7 in Figure 6). It should be noted that Popen is close to Pten2 at the M-H boundary, so the NF1 shear failure and the NF2 tensile failure may occur at the same time nearby. In addition, the failure mechanism of the crack is different even if the crack propagates in the same form. The concrete failure form should be judged according to its failure mechanism.From Figure6, it can be seen that in the low-angle region or near the M-H boundary, the crack number changes from 1 to 3. In the middle-angle area, the crack number is changed from the original 1 to 2. In both cases, the number of hydraulic cracks crossing the intersection point is higher than that before crossing. Whether the number of cracks increases or not is of great significance to the formation of a fracture network after fracturing. Therefore, it is more advantageous to form a complex fracture network at low and medium angles. When the angle is high, the cracks extend forward directly and the number of cracks does not increase, so it is not conducive to the formation of the fracture network. In summary, it is most advantageous to form a fracture network in the low intersection angle and angle boundary range, followed by a medium angle, and difficult to form an effective fracture network at a high angle. If the range of favorable angles can be widened, the effect of fracturing can be improved. Therefore, this paper carries out parameter analysis to clarify the influence of different parameters on the L-M boundary and M-H boundary.


## 3.1. Fracture Propagation in Rocks with Cemented Joints

The time-dependent injection pressure was recorded. The result is shown in Figure3. As shown in the figure, the pressure increases rapidly at the initial period. The crack was initialized when it reached maximum pressure at 56.4 MPa. Then, there is a vibration period after crack open. At last, the pressure tended to become smooth and steady. This is consistent with the trend of actual monitoring data on site, which means that the simulation results are suitable for reference research.Figure 3
The history of injection pressure with time.One of the main purposes of numerical simulation is to investigate the behavior of a fluid-driven crack in the rock with cemented joints. In other words, the shape of open fractures after fracking is one of the most important goals. The shape is described by the aperture of fractures in this study. The results are plotted in Figure4. Figures 4(a)–4(c) is the spatial distribution of fractures with different angles between the vertical joint and two cutters. As shown in the pictures, with different intersection angles, the final shape is varied. According to a predecessor’s studies, it is easy to understand that the crack extends straight forward with a high intersection angle, while the crack is always captured by the meeting natural fracture but would not stretch forward when the intersection angle is large enough. However, at the situation that two angles are 30°, it penetrates again. Besides, it is also captured by natural fractures. Why can the fluid-driven crack penetrate again? This is determined by the stress field around the fractures. In Section 3.2, the fracture propagation criterion based on elastic mechanics is discussed.Figure 4
The aperture of joints: the spatial distribution of fractures with different angles between the vertical joint and two cutters: (a) 30-30, (b) 45-45, and (c) 90-90; the aperture of vertical joints with different angles between the vertical joint and two cutters: (d) 30-30, (e) 45-45, and (f) 90-90.
(a)
(b)
(c)
(d)
(e)
(f)Another noteworthy feature different from the 2D model simulation is that fracture height is restrained by two cutters. In the 2D model, the fluid-driven crack must extend forward in their own form, penetrating or captured, when approaching two joints. However, in the 3D model, it is not obligatory for the fluid-driven crack extending forward as in the 2D. Figure4(d)–4(f) is the aperture of vertical joints with different angles at 30°, 60°, and 90°, respectively. As shown in the figures, no matter what the angles, the vertical propagation is restricted by the two cutters. And it extends left and right. The fracture containment problem is another problem still to be solved.


## 3.2. Critical Pressure for the Failure of the Weakness Plane

There are many weak surface structures, such as bedding and fractures, in shale due to the action of geological tectonic stress. It is more likely to be damaged in fracturing for the reason that the strength of the weakness plane is low. The fracture will be preferentially expanded along the weakness plane. Therefore, the intersection of the hydraulic fracture and the weakness plane is the key problem of fracture propagation in the formation with the weakness plane.The intersection model of the hydraulic fracture and the natural fracture is shown in Figure5. Suppose there are two intersecting closed natural fractures, NF1 and NF2, at an approach angle of α. The assumptions in this paper are as follows. The hydraulic fracture extends along the fracture plane NF2 parallel to the maximum principal stress. The hydraulic fracture tip reaches the intersection of the two fractures but has not yet passed the intersection point. NF1 and NF2 have certain tensile and shear strength before cracking, and the tensile and the shear strength disappear after cracking. To simplify the calculation, this paper selects the same method as the W&T criterion and assumes that the fracture is blunted at the interface, so that the singularity of the hydraulic fracture tip can be ignored. More detailed information can be found in [5].Figure 5
The intersection model of the hydraulic fracture and the natural fracture.The failure modes of NF1 and NF2 under fluid action are different with different angles of fracture and principal stress. In general, the failure form of NF1 may be a tensile failure or a shear failure. After tensile failure, cracks will open directly. However, cracks do not necessarily open after shear failure. It also requires that the fluid pressure can overcome the normal stress on the fracture surface and prop up the fracture. The failure mode of NF2 without shear stress on the fracture surface can only be tensile failure. The failure forms are summarized in Table2. Different failure forms should meet the following conditions.Table 2
The criteria for different failure forms of fractures.
Fracture
Failure forms
Criteria
NF1
Shear failure and opening
P=P∣P>Pshear∪P>Popen
Tensile failure
P=P∣P>Pten1
NF2
Tensile failure
P=P∣P>Pten2P is the fluid pressure at the hydraulic fracture tip. Pshear is the critical pressure to shear failure of NF1. Popen is the opening pressure of NF1 without bond strength. Pten1 is the critical tensile failure pressure of NF1 considering the bond strength. Pten2 is the critical tensile failure pressure of NF2 considering the bond strength.Once the critical conditions are determined, the failure forms of fracture propagation can be determined after obtaining the failure conditions of different forms of fractures. Therefore, the critical pressure of different fracture failure forms will be calculated in the following part of this section.According to Warpinski and Teufel, the critical fluid pressure in the crack during shear failure of NF1 can be calculated.(11)pshear=2τ0−σ1−σ3sin2α+kfcos2β−kf2kf+σ3,where σn is the normal stress on the fracture surface, τ is the shear stress on the fracture surface, σ1 is the maximum principal stress, σ3 is the minimum principal stress, α is the approach angle, τ0 is the inherent shear strength of the interface, and kf is the coefficient of friction.The strength of the weakness plane disappears after cracking. When the fluid pressure inside the crack is greater than the normal stress on the fracture surface, the fracture surface opens. Therefore, the fluid pressure in the crack after shear failure is(12)popen=σn=σ1sin2α+σ3cos2α=σ1−σ3sin2α+σ3.Tensile failure occurs when the fluid pressure in NF1 is greater than the sum of normal stress and tensile strength of the fracture surface. At this point, the critical pressure of fluid in the joint is(13)pten1=σn+σt=σ1−σ3sin2α+σ3+σt,where σt is the tensile strength of the cemented fracture.Tensile failure occurs when the fluid pressure of fracture NF2 is greater than the sum of normal stress on the surface and tensile strength of the fracture. At this time, the critical fluid pressure of tensile failure for fracture NF2 is(14)pten2=σ3+σt.The pressure conditions for crack opening can be determined by the above formulas. When the natural fracture NF1 opens, the hydraulic fracture is captured by fracture NF1. When fracture NF2 opens, the hydraulic fracture passes through fracture NF1. The critical pressure of fracture failure at different angles of NF1 and NF2 is calculated to determine the fracture expansion form at the intersection point.


## 3.3. The Behavior of Hydraulic Fracture When Approaching Intersecting Cemented Joints

The critical pressure of different failure forms of cracks can be calculated by the formulas (11)–(14) and Table 2. The parameters selected in this article are shown in Table 3. And the results are shown in Figure 6.Table 3
The values of parameters.
Parameter
Value
Parameter
Value
Inherent shear strength,τ0
2 MPa
Coefficient of friction,kf
0.2
Maximum principle stress,σ1
20 MPa
Tensile strength of cemented joint,σt
2 MPa
Minimum principle stress,σ3
12 MPaFigure 6
The critical pressure under different failure conditions.Obviously, for NF1, the critical pressure of shear failure decreases first and then increases with the increase of the intersecting angle. The critical pressure of tensile failure and opening after shear failure increase with the increase of the intersecting angle. However, the critical pressure of the NF2 tensile failure is independent of the intersection angle.At a small angle (0° and 10.67°), as the pressure in the crack rises, fluid pressure first reachedPopen. However, NF1 is not broken by the shear at this time, so it cannot be opened. The pressure then reached Pten1 and Pten2. At this point, the two values are close; both NF1 and NF2 undergo tensile failure. Hydraulic fractures extend to both NF1 and NF2, forming cross fractures. At the angle range of 10.67° to 15°, the fracture fluid pressure in the first achieves Popen, and no fractures are failures now. Then, the pressure continued to rise to Pshear and shear failure occurred in NF1. At this point, NF1 shear failure and cracking occurred and hydraulic fracture only extended in NF1. When the fracture intersection angle is 15° to 37.7°, the fluid pressure first reaches Pshear, and then NF1 shear failure occurs. At this point, the pressure has not reached the opening pressure Popen. As the pressure rises to Popen, NF1, which has been shear damaged, cracks. When the intersection angle is 37.7° to 68°, the shear failure of NF1 occurs first, but it does not open, because the fluid pressure at this time is less than Popen. With the increase of pressure, the fluid pressure reaches Pten2 and the tensile failure of NF2 occurs. Since NF1 is not open, hydraulic cracks will only propagate along NF2. When the intersection angle is between 68 and 90 degrees, the fluid pressure first reaches Pten2 and the hydraulic fracture will only extend in NF2. It is noteworthy that the values of Pten2 and Popen are close to each other in the range of 37.7 degrees, at which time NF1 and NF2 will open simultaneously.To sum up, with the change of the intersecting angle, the expansion morphology of intersecting cracks can be divided into three categories. At low angles (0-10.67), NF1 and NF2 open at the same time, forming a cross-state. At the middle angle (10.67-37.7), shear cracking occurred in NF1 and no cracking occurred in NF2. The fracture morphology shows turning. At high angles (37.7-90), NF2 cracks and NF1 does not open. The fracture shape is a single fracture extending forward. In summary, the intersection angle can be divided into three parts according to the fracture propagation pattern, namely, low angle, medium angle, and high angle. The boundary between the low angle and the middle angle is defined as theL-M boundary (10.67 in Figure 6). And the boundary between the medium angle and the high angle is defined as the M-H boundary (37.7 in Figure 6). It should be noted that Popen is close to Pten2 at the M-H boundary, so the NF1 shear failure and the NF2 tensile failure may occur at the same time nearby. In addition, the failure mechanism of the crack is different even if the crack propagates in the same form. The concrete failure form should be judged according to its failure mechanism.From Figure6, it can be seen that in the low-angle region or near the M-H boundary, the crack number changes from 1 to 3. In the middle-angle area, the crack number is changed from the original 1 to 2. In both cases, the number of hydraulic cracks crossing the intersection point is higher than that before crossing. Whether the number of cracks increases or not is of great significance to the formation of a fracture network after fracturing. Therefore, it is more advantageous to form a complex fracture network at low and medium angles. When the angle is high, the cracks extend forward directly and the number of cracks does not increase, so it is not conducive to the formation of the fracture network. In summary, it is most advantageous to form a fracture network in the low intersection angle and angle boundary range, followed by a medium angle, and difficult to form an effective fracture network at a high angle. If the range of favorable angles can be widened, the effect of fracturing can be improved. Therefore, this paper carries out parameter analysis to clarify the influence of different parameters on the L-M boundary and M-H boundary.


## 4. Parameter Analysis

### 4.1. Effect of Minimum Principal Stress

From (11)–(14), we can see that σ3 has a separate term in each formula. When the difference of in situ stress is constant, if only the minimum principal stress is changed, the critical pressure will change the same value and their difference will not change. Therefore, the change of the minimum horizontal principal stress will not cause the change of the angle range. Therefore, in the case of constant in situ stress difference, the influence of the minimum principal stress on the angle range division of fracture propagation can be neglected. However, the increase of the minimum principal stress will lead to an increase of the critical pressure, so the increase of the minimum principal stress will increase the difficulty of crack opening.


### 4.2. The Effect of Inherent Shear Strength of the Interface

The inherent shear strength influences the shear behavior of cracks. The higher the inherent shear strength is, the greater the shear strength of cracks is. According to the propagation criterion of intersecting cemented cracks, the inherent shear strength only affectsPshear and the two are positively correlated.The two angle boundaries vary withτ0 as shown in Figure 7. It can be seen that with the increase of τ0, the L-M boundary value increases gradually, while the M-H boundary keeps constant. Figure 7(a) shows that when τ0 is small, the L-M boundary is close to 0, which is due to the shear failure of NF1 caused by the decrease of Pshear, leading to the left movement of the intersection of Pshear and Pten2. When τ0 is too large, the shear strength of NF1 is too high and the shear failure of NF1 will not occur. As can be seen from Figure 7(b), when τ0=4, Pshear and Pten2 do not intersect and the middle-angle range disappears. In short, as τ0 increases, the L-M boundary moves to the right and the low-angle range expands. The fracture network after fracturing is the most complex in the low-angle range. Therefore, the increase of τ0 can improve the complexity of the fracture network after fracturing.Figure 7
The effect of the inherent shear strength of the interface.
(a)
(b)


### 4.3. The Effect of the Friction Coefficient

The influence of the friction coefficient on each boundary is shown in Figure8(a). As can be seen from the figure, the friction coefficient only affects the L-M boundary. With the increase of kf, the L-M boundary moves to the left and the range of the middle angle becomes wider. The changes of critical pressures under different kf conditions are shown in Figure 8(b). As shown in the figure, as kf increased, the Pshear curve increased in the middle part while it decreased in the two ends. This causes the left intersection point of curve Pshear and Pten2 to move to the left, resulting in the decrease of the L-M value.Figure 8
The effect of the friction coefficient.
(a)
(b)


### 4.4. The Effect of the Differential Stress

As mentioned in Section4.1, the change of minimum principal stress has no effect on the boundary ranges. It is the differential stresses that affect the boundary range. The critical pressures (Figures 9(a)–9(c)) under different differential stress conditions (0 MPa, 2 MPa, and 5 MPa) were calculated, and the change curves of the L-M boundary and the M-H boundary with differential stress under 4-10 MPa were given (Figure 9(d)).Figure 9
The effect of the differential stress.
(a)
(b)
(c)
(d)Figure9(a) is the curve of each critical pressure when the differential stress is 0. It can be seen from the figure that under the condition of uniform stress, the intersection angle has no influence on the critical pressure. Moreover, Pten1 and Pten2 were equal and less than Pshear. Under this condition, both NF1 and NF2 will undergo tensile failure at the same time, forming cross-cracks. This indicates that under the condition of uniform stress, it will be advantageous to form complex fracture networks. Figure 9(b) shows the critical pressure curve when the stress difference is 2 MPa. It can be seen from the figure that under the condition of nonuniform stress field, the critical conditions of NF1 will change with the angle of intersection. The Pshear was lower than that under uniform stress. This indicates that the increase of differential stress reduces the difficulty of the shear failure of NF1. Under the condition of lower differential stress, as Pshear was greater than Pten1 and Pten2, NF1 had not yet undergone shear failure. Therefore, in the case of low angle, NF1 and NF2 suffer from tensile failure at the same time, while in the case of high angle, only NF2 suffers from tensile failure. Figure 9(c) shows the critical pressure curves when the differential stress is 5 MPa. According to the figure, with the increase of the different stress, the critical financial pressure curve will become a classical three-interval model. Figure 9(d) shows the boundary change with the change of differential stress. As can be seen from the figure, both the L-M boundary and the M-H boundary decrease with the increase of differential stress. Since the angle range larger than the M-H boundary cannot form an effective fracture network, it indicates that it is more difficult to form a fracture network under the condition of high differential stress. In addition, when the boundary is less than the L-M boundary, the number of final cracks formed is the largest. The reduction of the L-M boundary will lead to the reduction of the complexity of the joint network. In conclusion, under the condition of uniform stress, a complex joint network will be formed. Low differential stress conditions are conducive to the formation of complex fracture networks, while high stress conditions are not conducive to the formation of fracture networks. The complexity of the fracture network formed under high differential stress is reduced.


### 4.5. The Effect of the Tensile Strength of Fractures

The tensile strength of cracks affects their tensile cracking. The greater the tensile strength, the higher the critical pressure for tensile cracking. The influence of the fracture tensile strength on each critical pressure is shown in Figure10. Figure 10(a) shows the change of the boundary with the tensile strength. It can be seen from the figure that with the increase of the tensile strength, the L-M boundary decreases slightly, the M-H boundary increases, and the range of the middle angle expands. This indicates that the increase of the fracture tensile strength is conducive to the formation of the fracture network. According to Figure 10(b), as the fracture strength increased, the values of Pten1 and Pten2 increased, while Popen and Pshear remained constant. As the difficulty of fracture tensile failure increases, shear failure is more likely to occur under the same pressure. Shear failure is the main cause of fracture steering. Therefore, the appropriate increase of fracture tensile strength is conducive to fracture diversion. But increased tensile strength also makes fracturing more difficult.Figure 10
The effect of the tensile strength of fractures.
(a)
(b)


### 4.6. Multivariate Regression Analysis Based On the Response Surface Method

As mentioned in Section3.3, the premise of forming a complex fracture network is the increase of the number of fractures. When the angle is less than the L-M boundary, the number of fractures changes from 1 to 3, and in the middle-angle range, the number of fractures changes from 1 to 2, both of which are conducive to the formation of the complex fracture network. So it makes sense to study the boundary. The above research can analyze the influence of a single factor on the two kinds of the boundary but cannot comprehensively evaluate the comprehensive influence of multiple factors. The response surface method can be used to study the dependence of the response variable R on independent variables. If the number of independent variables is small (no more than 3), then the response surface method is one of the best methods. The response surface method can not only study the influence of individual variables on response variables but also study the influence of interaction between independent variables on response variables. The general form of the linear regression equation fitted by the response surface method is
(15)R=∑i=1nkixi+∑i=1n∑j=1nkijxixj+k0,where n is the number of variables, xi is the variable, and ki is the coefficient to be solved.The selection and combination of independent variables are called the experimental design. The central composite design (CCD) is a common method for experimental design. The design tables of two experiments are shown in Tables4 and 5. Specific variable selection will be based on different response variables.Table 4
Design of experiment and response value for the L-M boundary.
No.
Factor 1
Factor 2
Factor 3
Response 1
A
B
C
R1
1
8
2.4
0.3
14.62
2
8
2
0.5
7.79
3
8
2
0.3
10.91
4
6.4
2
0.3
14.13
5
8
1.6
0.3
7.55
6
8.95137
2.23784
0.418921
9.86
7
8
2
0.3
10.91
8
8
2
0.3
10.91
9
8
2
0.1
13.74
10
7.04863
2.23784
0.181079
17.16
11
8
2
0.3
10.91
12
7.04863
1.76216
0.181079
12.22
13
9.6
2
0.3
8.92
14
7.04863
2.23784
0.418921
13.06
15
8.95137
2.23784
0.181079
12.97
16
8.95137
1.76216
0.418921
6.26
17
8.95137
1.76216
0.181079
9.41
18
7.04863
1.76216
0.418921
8.10
19
8
2
0.3
10.91
20
8
2
0.3
10.91Table 5
Design of experiment and response value for the M-H boundary.
No.
Factor 1
Factor 2
Response 1
A
B
R1
1
6.86863
2.28284
35
2
8
2
30
3
8
2
30
4
9.6
2
27
5
8
2
30
6
8
1.6
27
7
8
2
30
8
8
2.4
33
9
9.13137
2.28284
30
10
9.13137
1.71716
26
11
6.86863
1.71716
30
12
6.4
2
34
13
8
2
30


#### 4.6.1. L-M Boundary

According to the previous analysis, for the L-M boundary, the influence of the tensile strength of joints is negligible. Therefore, in the analysis of the L-M boundary, three independent variables,A, B, and C, were selected—A: principal stress difference; B: inherent shear strength, τ0; and C: coefficient of friction, kf. The selection of variables and the experimental design are shown in Tables 4 and 5. Then, the RSM for the L-M boundary was conducted. The regression equation is shown in equation (16). And the analysis of variance is listed in Table 6.
(16)R1=10.73−3.04A+16.69B−30.68C−1.51AB+2.16AC+0.26BC+0.24A2+1.06B2−3.74C2.Table 6
Analysis of variance table.
Source
Sum of squares
Mean square
F-value
p value
Prob>F
Model
139.16
15.46
10080.89
<0.0001
Significant
A
31.69
31.69
20658.06
<0.0001
B
61.37
61.37
40010.96
<0.0001
C
43.90
43.90
28624.17
<0.0001
AB
0.94
0.94
611.84
<0.0001
AC
0.48
0.48
313.07
<0.0001
BC
4.5e-4
4.5e-4
0.29
0.5999
A2
0.67
0.67
437.30
<0.0001
B2
0.052
0.052
34.02
0.0002
C2
0.040
0.040
26.35
0.0004The modelF-value of 10080.89 implies the model is significant. There is only a 0.01% chance that an F-value this large could occur due to noise. Values of “Prob>F” less than 0.0500 indicate that model terms are significant. And the significance order is B>C>A.The factor effect graph shows the linear effect of changing the level of a single factor. It is constructed by predicting the responses for the low (-1) and high (+1) levels of a factor. As shown in Figure11, the L-M boundary is negatively correlated with principal stress difference and coefficient of friction and positively correlated with inherent shear strength. The result is consistent with the univariate analysis above.Figure 11
Factor effect graph of variances for the L-M boundary.
(a)
(b)
(c)


#### 4.6.2. M-H Boundary

From the univariate analysis, it is known that the M-H boundary is not affected byτ0 and kf. Therefore, in the analysis of the M-H boundary, only two independent variables A and B are used for analysis. They are A: principal stress difference, and B: tensile strength of cemented joint, σt. The selection of variables and the experimental design are shown in Tables 4 and 5. Then, the RSM for the L-M boundary was conducted. The regression equation is shown in equation (17). And the analysis of variance is listed in Table 7.
(17)R2=30.13−4.23A+16.91B−0.70AB+0.22A2−0.73B2.Table 7
Analysis of variance table.
Source
Sum of squares
Mean square
F-value
p value
Prob>F
Model
91.47
18.29
35479.45
<0.0001
Significant
A
45.98
45.98
89173.59
<0.0001
B
44.67
44.67
86638.72
<0.0001
AB
0.20
0.20
392.73
<0.0001
A2
0.55
0.55
1067.20
<0.0001
B2
0.024
0.024
46.57
0.0002The modelF-value of 35479.45 implies the model is significant. There is only a 0.01% chance that an F-value this large could occur due to noise. Values of “Prob>F” less than 0.0500 indicate model terms are significant. In this case A, B, AB, A2, and B2 are significant model terms. And the significance order is A>B.The factor effect graph of variances for the M-H boundary is shown in Figure12. The M-H boundary is negatively correlated with principal stress difference and positively correlated with the tensile strength of the cemented joint. The result is consistent with the univariate analysis above.Figure 12
Factor effect graph of variances for the M-H boundary.


## 4.1. Effect of Minimum Principal Stress

From (11)–(14), we can see that σ3 has a separate term in each formula. When the difference of in situ stress is constant, if only the minimum principal stress is changed, the critical pressure will change the same value and their difference will not change. Therefore, the change of the minimum horizontal principal stress will not cause the change of the angle range. Therefore, in the case of constant in situ stress difference, the influence of the minimum principal stress on the angle range division of fracture propagation can be neglected. However, the increase of the minimum principal stress will lead to an increase of the critical pressure, so the increase of the minimum principal stress will increase the difficulty of crack opening.


## 4.2. The Effect of Inherent Shear Strength of the Interface

The inherent shear strength influences the shear behavior of cracks. The higher the inherent shear strength is, the greater the shear strength of cracks is. According to the propagation criterion of intersecting cemented cracks, the inherent shear strength only affectsPshear and the two are positively correlated.The two angle boundaries vary withτ0 as shown in Figure 7. It can be seen that with the increase of τ0, the L-M boundary value increases gradually, while the M-H boundary keeps constant. Figure 7(a) shows that when τ0 is small, the L-M boundary is close to 0, which is due to the shear failure of NF1 caused by the decrease of Pshear, leading to the left movement of the intersection of Pshear and Pten2. When τ0 is too large, the shear strength of NF1 is too high and the shear failure of NF1 will not occur. As can be seen from Figure 7(b), when τ0=4, Pshear and Pten2 do not intersect and the middle-angle range disappears. In short, as τ0 increases, the L-M boundary moves to the right and the low-angle range expands. The fracture network after fracturing is the most complex in the low-angle range. Therefore, the increase of τ0 can improve the complexity of the fracture network after fracturing.Figure 7
The effect of the inherent shear strength of the interface.
(a)
(b)


## 4.3. The Effect of the Friction Coefficient

The influence of the friction coefficient on each boundary is shown in Figure8(a). As can be seen from the figure, the friction coefficient only affects the L-M boundary. With the increase of kf, the L-M boundary moves to the left and the range of the middle angle becomes wider. The changes of critical pressures under different kf conditions are shown in Figure 8(b). As shown in the figure, as kf increased, the Pshear curve increased in the middle part while it decreased in the two ends. This causes the left intersection point of curve Pshear and Pten2 to move to the left, resulting in the decrease of the L-M value.Figure 8
The effect of the friction coefficient.
(a)
(b)


## 4.4. The Effect of the Differential Stress

As mentioned in Section4.1, the change of minimum principal stress has no effect on the boundary ranges. It is the differential stresses that affect the boundary range. The critical pressures (Figures 9(a)–9(c)) under different differential stress conditions (0 MPa, 2 MPa, and 5 MPa) were calculated, and the change curves of the L-M boundary and the M-H boundary with differential stress under 4-10 MPa were given (Figure 9(d)).Figure 9
The effect of the differential stress.
(a)
(b)
(c)
(d)Figure9(a) is the curve of each critical pressure when the differential stress is 0. It can be seen from the figure that under the condition of uniform stress, the intersection angle has no influence on the critical pressure. Moreover, Pten1 and Pten2 were equal and less than Pshear. Under this condition, both NF1 and NF2 will undergo tensile failure at the same time, forming cross-cracks. This indicates that under the condition of uniform stress, it will be advantageous to form complex fracture networks. Figure 9(b) shows the critical pressure curve when the stress difference is 2 MPa. It can be seen from the figure that under the condition of nonuniform stress field, the critical conditions of NF1 will change with the angle of intersection. The Pshear was lower than that under uniform stress. This indicates that the increase of differential stress reduces the difficulty of the shear failure of NF1. Under the condition of lower differential stress, as Pshear was greater than Pten1 and Pten2, NF1 had not yet undergone shear failure. Therefore, in the case of low angle, NF1 and NF2 suffer from tensile failure at the same time, while in the case of high angle, only NF2 suffers from tensile failure. Figure 9(c) shows the critical pressure curves when the differential stress is 5 MPa. According to the figure, with the increase of the different stress, the critical financial pressure curve will become a classical three-interval model. Figure 9(d) shows the boundary change with the change of differential stress. As can be seen from the figure, both the L-M boundary and the M-H boundary decrease with the increase of differential stress. Since the angle range larger than the M-H boundary cannot form an effective fracture network, it indicates that it is more difficult to form a fracture network under the condition of high differential stress. In addition, when the boundary is less than the L-M boundary, the number of final cracks formed is the largest. The reduction of the L-M boundary will lead to the reduction of the complexity of the joint network. In conclusion, under the condition of uniform stress, a complex joint network will be formed. Low differential stress conditions are conducive to the formation of complex fracture networks, while high stress conditions are not conducive to the formation of fracture networks. The complexity of the fracture network formed under high differential stress is reduced.


## 4.5. The Effect of the Tensile Strength of Fractures

The tensile strength of cracks affects their tensile cracking. The greater the tensile strength, the higher the critical pressure for tensile cracking. The influence of the fracture tensile strength on each critical pressure is shown in Figure10. Figure 10(a) shows the change of the boundary with the tensile strength. It can be seen from the figure that with the increase of the tensile strength, the L-M boundary decreases slightly, the M-H boundary increases, and the range of the middle angle expands. This indicates that the increase of the fracture tensile strength is conducive to the formation of the fracture network. According to Figure 10(b), as the fracture strength increased, the values of Pten1 and Pten2 increased, while Popen and Pshear remained constant. As the difficulty of fracture tensile failure increases, shear failure is more likely to occur under the same pressure. Shear failure is the main cause of fracture steering. Therefore, the appropriate increase of fracture tensile strength is conducive to fracture diversion. But increased tensile strength also makes fracturing more difficult.Figure 10
The effect of the tensile strength of fractures.
(a)
(b)


## 4.6. Multivariate Regression Analysis Based On the Response Surface Method

As mentioned in Section3.3, the premise of forming a complex fracture network is the increase of the number of fractures. When the angle is less than the L-M boundary, the number of fractures changes from 1 to 3, and in the middle-angle range, the number of fractures changes from 1 to 2, both of which are conducive to the formation of the complex fracture network. So it makes sense to study the boundary. The above research can analyze the influence of a single factor on the two kinds of the boundary but cannot comprehensively evaluate the comprehensive influence of multiple factors. The response surface method can be used to study the dependence of the response variable R on independent variables. If the number of independent variables is small (no more than 3), then the response surface method is one of the best methods. The response surface method can not only study the influence of individual variables on response variables but also study the influence of interaction between independent variables on response variables. The general form of the linear regression equation fitted by the response surface method is
(15)R=∑i=1nkixi+∑i=1n∑j=1nkijxixj+k0,where n is the number of variables, xi is the variable, and ki is the coefficient to be solved.The selection and combination of independent variables are called the experimental design. The central composite design (CCD) is a common method for experimental design. The design tables of two experiments are shown in Tables4 and 5. Specific variable selection will be based on different response variables.Table 4
Design of experiment and response value for the L-M boundary.
No.
Factor 1
Factor 2
Factor 3
Response 1
A
B
C
R1
1
8
2.4
0.3
14.62
2
8
2
0.5
7.79
3
8
2
0.3
10.91
4
6.4
2
0.3
14.13
5
8
1.6
0.3
7.55
6
8.95137
2.23784
0.418921
9.86
7
8
2
0.3
10.91
8
8
2
0.3
10.91
9
8
2
0.1
13.74
10
7.04863
2.23784
0.181079
17.16
11
8
2
0.3
10.91
12
7.04863
1.76216
0.181079
12.22
13
9.6
2
0.3
8.92
14
7.04863
2.23784
0.418921
13.06
15
8.95137
2.23784
0.181079
12.97
16
8.95137
1.76216
0.418921
6.26
17
8.95137
1.76216
0.181079
9.41
18
7.04863
1.76216
0.418921
8.10
19
8
2
0.3
10.91
20
8
2
0.3
10.91Table 5
Design of experiment and response value for the M-H boundary.
No.
Factor 1
Factor 2
Response 1
A
B
R1
1
6.86863
2.28284
35
2
8
2
30
3
8
2
30
4
9.6
2
27
5
8
2
30
6
8
1.6
27
7
8
2
30
8
8
2.4
33
9
9.13137
2.28284
30
10
9.13137
1.71716
26
11
6.86863
1.71716
30
12
6.4
2
34
13
8
2
30


### 4.6.1. L-M Boundary

According to the previous analysis, for the L-M boundary, the influence of the tensile strength of joints is negligible. Therefore, in the analysis of the L-M boundary, three independent variables,A, B, and C, were selected—A: principal stress difference; B: inherent shear strength, τ0; and C: coefficient of friction, kf. The selection of variables and the experimental design are shown in Tables 4 and 5. Then, the RSM for the L-M boundary was conducted. The regression equation is shown in equation (16). And the analysis of variance is listed in Table 6.
(16)R1=10.73−3.04A+16.69B−30.68C−1.51AB+2.16AC+0.26BC+0.24A2+1.06B2−3.74C2.Table 6
Analysis of variance table.
Source
Sum of squares
Mean square
F-value
p value
Prob>F
Model
139.16
15.46
10080.89
<0.0001
Significant
A
31.69
31.69
20658.06
<0.0001
B
61.37
61.37
40010.96
<0.0001
C
43.90
43.90
28624.17
<0.0001
AB
0.94
0.94
611.84
<0.0001
AC
0.48
0.48
313.07
<0.0001
BC
4.5e-4
4.5e-4
0.29
0.5999
A2
0.67
0.67
437.30
<0.0001
B2
0.052
0.052
34.02
0.0002
C2
0.040
0.040
26.35
0.0004The modelF-value of 10080.89 implies the model is significant. There is only a 0.01% chance that an F-value this large could occur due to noise. Values of “Prob>F” less than 0.0500 indicate that model terms are significant. And the significance order is B>C>A.The factor effect graph shows the linear effect of changing the level of a single factor. It is constructed by predicting the responses for the low (-1) and high (+1) levels of a factor. As shown in Figure11, the L-M boundary is negatively correlated with principal stress difference and coefficient of friction and positively correlated with inherent shear strength. The result is consistent with the univariate analysis above.Figure 11
Factor effect graph of variances for the L-M boundary.
(a)
(b)
(c)


### 4.6.2. M-H Boundary

From the univariate analysis, it is known that the M-H boundary is not affected byτ0 and kf. Therefore, in the analysis of the M-H boundary, only two independent variables A and B are used for analysis. They are A: principal stress difference, and B: tensile strength of cemented joint, σt. The selection of variables and the experimental design are shown in Tables 4 and 5. Then, the RSM for the L-M boundary was conducted. The regression equation is shown in equation (17). And the analysis of variance is listed in Table 7.
(17)R2=30.13−4.23A+16.91B−0.70AB+0.22A2−0.73B2.Table 7
Analysis of variance table.
Source
Sum of squares
Mean square
F-value
p value
Prob>F
Model
91.47
18.29
35479.45
<0.0001
Significant
A
45.98
45.98
89173.59
<0.0001
B
44.67
44.67
86638.72
<0.0001
AB
0.20
0.20
392.73
<0.0001
A2
0.55
0.55
1067.20
<0.0001
B2
0.024
0.024
46.57
0.0002The modelF-value of 35479.45 implies the model is significant. There is only a 0.01% chance that an F-value this large could occur due to noise. Values of “Prob>F” less than 0.0500 indicate model terms are significant. In this case A, B, AB, A2, and B2 are significant model terms. And the significance order is A>B.The factor effect graph of variances for the M-H boundary is shown in Figure12. The M-H boundary is negatively correlated with principal stress difference and positively correlated with the tensile strength of the cemented joint. The result is consistent with the univariate analysis above.Figure 12
Factor effect graph of variances for the M-H boundary.


## 4.6.1. L-M Boundary

According to the previous analysis, for the L-M boundary, the influence of the tensile strength of joints is negligible. Therefore, in the analysis of the L-M boundary, three independent variables,A, B, and C, were selected—A: principal stress difference; B: inherent shear strength, τ0; and C: coefficient of friction, kf. The selection of variables and the experimental design are shown in Tables 4 and 5. Then, the RSM for the L-M boundary was conducted. The regression equation is shown in equation (16). And the analysis of variance is listed in Table 6.
(16)R1=10.73−3.04A+16.69B−30.68C−1.51AB+2.16AC+0.26BC+0.24A2+1.06B2−3.74C2.Table 6
Analysis of variance table.
Source
Sum of squares
Mean square
F-value
p value
Prob>F
Model
139.16
15.46
10080.89
<0.0001
Significant
A
31.69
31.69
20658.06
<0.0001
B
61.37
61.37
40010.96
<0.0001
C
43.90
43.90
28624.17
<0.0001
AB
0.94
0.94
611.84
<0.0001
AC
0.48
0.48
313.07
<0.0001
BC
4.5e-4
4.5e-4
0.29
0.5999
A2
0.67
0.67
437.30
<0.0001
B2
0.052
0.052
34.02
0.0002
C2
0.040
0.040
26.35
0.0004The modelF-value of 10080.89 implies the model is significant. There is only a 0.01% chance that an F-value this large could occur due to noise. Values of “Prob>F” less than 0.0500 indicate that model terms are significant. And the significance order is B>C>A.The factor effect graph shows the linear effect of changing the level of a single factor. It is constructed by predicting the responses for the low (-1) and high (+1) levels of a factor. As shown in Figure11, the L-M boundary is negatively correlated with principal stress difference and coefficient of friction and positively correlated with inherent shear strength. The result is consistent with the univariate analysis above.Figure 11
Factor effect graph of variances for the L-M boundary.
(a)
(b)
(c)


## 4.6.2. M-H Boundary

From the univariate analysis, it is known that the M-H boundary is not affected byτ0 and kf. Therefore, in the analysis of the M-H boundary, only two independent variables A and B are used for analysis. They are A: principal stress difference, and B: tensile strength of cemented joint, σt. The selection of variables and the experimental design are shown in Tables 4 and 5. Then, the RSM for the L-M boundary was conducted. The regression equation is shown in equation (17). And the analysis of variance is listed in Table 7.
(17)R2=30.13−4.23A+16.91B−0.70AB+0.22A2−0.73B2.Table 7
Analysis of variance table.
Source
Sum of squares
Mean square
F-value
p value
Prob>F
Model
91.47
18.29
35479.45
<0.0001
Significant
A
45.98
45.98
89173.59
<0.0001
B
44.67
44.67
86638.72
<0.0001
AB
0.20
0.20
392.73
<0.0001
A2
0.55
0.55
1067.20
<0.0001
B2
0.024
0.024
46.57
0.0002The modelF-value of 35479.45 implies the model is significant. There is only a 0.01% chance that an F-value this large could occur due to noise. Values of “Prob>F” less than 0.0500 indicate model terms are significant. In this case A, B, AB, A2, and B2 are significant model terms. And the significance order is A>B.The factor effect graph of variances for the M-H boundary is shown in Figure12. The M-H boundary is negatively correlated with principal stress difference and positively correlated with the tensile strength of the cemented joint. The result is consistent with the univariate analysis above.Figure 12
Factor effect graph of variances for the M-H boundary.


## 5. Conclusion

In order to analyze the mechanism of hydraulic fracturing fracture propagation in formations with cemented joints, the critical pressure of three fracture propagation modes under the intersection of hydraulic fracturing fracture and closed natural fracture is derived, and the parameter analysis is carried out. The main conclusions are as follows:(1)
The intersection angle can be divided into three parts according to the fracture propagation pattern, namely, low angle, medium angle, and high angle. Three types of behavior for a fracture are a cross-state at a low angle, turning at a middle angle, and extending forward at a high angle. The first two conditions are conducive to the formation of complex fracture networks
(2)
Univariate analysis was performed on each variable by the control variable method. The result shows that the increase ofτ0 and kf can improve the complexity of the fracture network after fracturing. And low differential stress conditions are conducive to the formation of complex fracture networks. Besides, the appropriate increase of the fracture tensile strength is conducive to fracture diversion
(3)
The relations between the two boundaries and different independent variables are obtained, and their regression formulas are obtained by fitting by RSM



---

*Source: 5406870-2019-06-27.xml*