#### Abstract

Defining the trajectory of hydraulic fractures crossing bedding planes and other fractures is a significant issue in determining the effectiveness of the stimulation. In this work, a damage evolution law is used to describe the initiation and propagation of the fracture. The model couples rock deformation and gas seepage using the finite element method and is validated against classical theoretical analysis. The simulation results define four basic intersection scenarios between the fluid-driven and preexisting fractures: (a) inserting—the hydraulic fracture inserts into a bedding plane and continues to propagate along it; (b) L-shaped crossing—the hydraulic fracture approaches the fracture/bedding plane then branches into the plane without crossing it; (c) T-shaped crossing—the hydraulic fracture approaches the fracture/bedding plane, branches into it, and crosses through it; (d) direct crossing—the hydraulic fracture crosses one or more bedding planes without branching into them. The intersection scenario changes from (a) → (b) → (c) → (d) in specimens with horizontal bedding planes when the stress ratio () increases from 0.2 to 5. Similarly, the intersection type changes from (d) → (c) → (a) with an increase in the bedding plane angle (0° → 90°). Stiffness of the bedding planes also exerts a significant influence on the propagation of hydraulic fractures. As the stiffness ratio increases from 0.1 to 0.4 and 0.8, the seepage area decreases from 22.2% to 41.8%, and the intersection type changes from a T-shaped crossing to a direct crossing.

#### 1. Introduction

The production of environmentally friendly unconventional gas has increased rapidly in recent years including shale gas, gas from tight sandstones, and coal bed methane [1–4], obviating the recovery of hydrocarbons from conventional reservoirs and coal [5, 6]. In the US, unconventional gas accounted for 46% of total production of natural gas in 2016 [7]. However, the permeability of these unconventional reservoirs is extremely low (less than 0.1 mD) [8, 9], necessitating massive hydraulic fracturing to improve production [10]. At present, water is the primary fracturing fluid due to its low cost and ready availability [11–15]. However, water-based fracturing may be constrained in water-deficient areas since fracturing treatments for a single reservoir may consume more than 227 m^{3} [16]. In addition, the disposal of the large amount of flowback water that is precharged with dissolved pollutants may adversely affect potable aquifers, induce seismicity [17], and pose a potential risk to public health [18]. Recently, laboratory experiments indicate that gas may be a superior candidate to water for fracturing. These results show that the breakdown pressures of gas fracturing are lower than those for water fracturing [19, 20]. Compared to fractures induced by water, those induced by gas have rougher surfaces, more tortuous paths, and a larger potential to form a complex fracture network, in turn potentially resulting in higher surface flow-transfer surface areas and reduced conductive flow lengths [10, 19, 21]. Additionally, the adsorption capacity of CO_{2} is 4-20 times of that of methane [22, 23], enabling the competitive replacement of methane by CO_{2} and the cosequestration of CO_{2}.

Shale is a sedimentary rock containing bedding planes [24], which exert a significant influence on the propagation of hydraulic fractures. Indeed, the complex fracture network created by hydraulic fracturing may, in part, result from the intersection scenarios between hydraulic fractures and bedding planes [25–28]. Thus, it is of vital importance to understand the behavior of hydraulic fractures intersecting bedding planes during propagation.

A variety of studies focus on the behavior of hydraulic fractures intersecting preexisting bedding planes, including theoretical analyses, laboratory experiments, and numerical simulations. Multiple crossing criteria have been proposed using different methods, such as linear elastic fracture mechanics [29], theory of dislocations [30], and systematic analysis [31].

Laboratory hydraulic fracturing experiments conducted on rock containing bedding planes investigate the controls on intersection behavior from different perspectives. Hydraulic fracturing experiments in prefractured shale indicate that hydraulic fractures tend to cross a preexisting fracture only under high differential stresses and at high angles of approach [32]. Three forms of intersections were obtained in experiments when hydraulic fractures propagate in rock with cemented natural fractures [33], including (1) the hydraulic fracture bypassing the natural fracture, (2) the hydraulic fracture arresting into and diverting along the natural fracture, a (3) a combination of bypass and diversion. Also, the impact of permeability, fluid flow rate, friction coefficient, stress coefficient, and form of the bedding planes are shown to impact the nature of the intersection and potential crossing [34–36].

Numerical simulation is an effective way to study the intersection relation between hydraulic fractures and bedding planes. Such approaches have explored the physics of fracture-inhomogeneity interactions, indicating that hydraulic fracture branching and diversion are the result of inhomogeneity [37, 38]. Investigations of the interaction between hydraulic fractures and bedding planes have utilized discrete element methods (DEM) [39, 40], displacement discontinuity methods (DDM) [41], and extended finite element methods (XFEM) [42]. However, the finite element method (FEM) combining with damage theory is also an effective way to simulate the fracture initiation and propagation. The intersection mechanism between the fluid-driven fracture and bedding planes using FEM is rarely reported.

In this work, a coupled hydraulic-mechanical model is proposed where a damage evolution law is employed. It is then solved by FEM using COMSOL and MATLAB and utilized to simulate the fracturing processes. The intersection scenarios between hydraulic fractures and bedding planes under several conditions are numerically researched.

#### 2. Governing Equations

We develop a damage-based hydraulic-mechanical model that follows the evolution of damage around a borehole with a propagating fluid-driven fracture. The model couples fluid and mechanical deformations to define the geometry of the resulting fluid-driven fracture.

##### 2.1. Rock Deformation Equation

According to the elastic theory, the constitutive equation considering the influence of pore pressure can be written as where is the stress tensor, is the shear modulus of rock, is the strain tensor, is Poisson’s ratio of rock, is the volumetric strain, is the Kronecker delta, is the Biot coefficient, and is the pore pressure.

Combining the modified constitutive equation, the geometric equation, and the equilibrium equation, the modified Navier-type equation can be written as where is the displacement vector and is the component of the body force.

##### 2.2. Gas Flow Equation

The governing equation for CO_{2} flow based on mass balance can be defined as
where is the CO_{2} mass per volume of rock, is the CO_{2} density, is the seepage velocity of CO_{2}, is the source of CO_{2}, and is the time variable.

According to Darcy’s law, the seepage velocity of CO_{2} can be defined as
where is the permeability of the shale rock and is the viscosity of CO_{2}.

The shale rock is assumed saturated with the injected CO_{2}; therefore, CO_{2} mass per volume of rock can be written as ( is the porosity of the rock). And CO_{2} will transfer from the gaseous to the supercritical state when the pressure reaches 7.43 MPa (the temperature is kept at 350 K). As shown in Figure 1, the density and viscosity of CO_{2} change dramatically when the phase change occurs. Base on the above, the first item of equation (3) can be induced as [13]
where is the compressibility coefficient of CO_{2} which can be calculated from Figure 1.

Substituting equations (4) and (5) into equation (3), the CO_{2} continuity equation is shown as

##### 2.3. Damage Evolution Law

A damage evolution law based on representative elemental volume (REV) is used in this study to describe the initiation and propagation of hydraulic fracture in numerical samples with bedding planes. The evolution of stress-strain of REV under uniaxial tension or compression is shown in Figure 2. Shear or tension damage is initiated when the stress state of a REV meets the Mohr-Coulomb criterion or the maximum tensile stress criterion, as expressed by where and are the first principal stress and third principal stress, respectively; and are the tensile strength and compressive strength of the REV, respectively; is the internal friction angle; and and are the threshold functions. Once a REV begins to get damage, the evolution of its stress-strain relation can be described by a nonlinear function as shown in Figure 2. A damage variable is used to represent the damage level and is defined as [43] where and are the first principal stain and the third principal strain, respectively, and and are the tensile strain and compressive strain, respectively.

For a damaged REV, the elastic modulus decreases but the permeability increases correspondingly as the damage variable increases. The evolution of elastic modulus and permeability with damage variable can be defined as follows: where and are the elastic modulus and the initial elastic modulus of a REV, respectively; is the initial permeability of a REV; and is the damage-permeability effect coefficient to define the influence of damage on permeability.

##### 2.4. Rock Heterogeneity

Shale is a kind of sedimentary rock which is heterogeneous. Previous studies indicated that the heterogeneity of rock plays an important role on the propagation of microcracks [44]. And many studies have shown that Weibull distribution can well characterize the heterogeneity of rock [13, 45]. In this work, the Weibull distribution is introduced to describe the heterogeneity of shale rock. Mechanical parameters of REVs (strength and elastic modulus) are assumed to satisfy the Weibull distribution function, and the probability density function is written as where represents the mechanical parameter of REVs, is the average value of the mechanical parameter, and is the heterogeneity coefficient.

#### 3. Verification and Implementation of the Proposed Numerical Model

We use a finite element method to solve the proposed numerical model and obtain numerical results, and we compare the numerical results with two classical analytical results to verify the effect of the proposed model.

##### 3.1. Model Implementation

The numerical model is established in the previous section. Due to the complex coupled relationship between solid mechanics field and fluid seepage field, it is difficult to obtain an analytical solution. Therefore, the finite element method (FEM) is adopted to solve these coupled equations. The primary procedures are summarized as shown in Figure 3: (a)After setting up the model geometry, the geometry is discretized into a series of REVs. Then, the initial mechanical parameters are defined at the REV scale and the boundary conditions are applied correspondingly(b)Numerical calculation is conducted by COMSOL Multiphysics at the initial load step. After calculation, the stress and strain of REVs are obtained for the following analysis(c)According to criterions (7) and (8), all REVs will be checked whether the damage occurs; this step is completed via MATLAB(d)The damage variable of damaged REVs can be calculated according to equation (9), then the elastic modulus and permeability of the damaged REVs are updated with equations (10) and (11)(e)Numerical simulation is conducted with the updated parameters, and the simulation results are compared with results of the former iteration step. If the damage area expands, steps (c)–(e) are repeated; otherwise, step (f) is applied(f)The boundary conditions are updated in the next load increment

##### 3.2. Model Verification

In this part, the proposed model for simulating the propagation of the hydraulic fracture is validated. There are two classical theoretical solutions for forecasting the breakdown pressure in terms of far-field stresses, tensile strength, and initial pore pressure. One is proposed by Hubbert and Willis [46] (the H-W solution) for impermeable rock whilst the other one is presented by Haimson and Fairhurst [47] (the H-F solution) for permeable rock. These two solutions can be expressed by where and are the breakdown pressures of the H-W solution and the H-F solution, respectively; is the tensile strength of the shale rock; and is the initial pore pressure of the shale rock.

In this part, the tensile strength of the shale rock is 6 MPa, the initial pore pressure of the rock is 1 MPa, the initial permeability of the rock is 10^{-18} m^{2}, the Biot coefficient is 0.1, and Poisson’s ratio of the rock is 0.225. The model geometry for verification is shown in Figure 4. The vertical in situ stress is kept at 30 MPa, while the horizontal in situ stress varies from 10 MPa to 30 MPa. And a tectonic stress ratio () is introduced. The breakdown pressures under different tectonic stress ratios obtained by the H-W solution, H-F solution, and numerical simulation are shown in Figure 5. The results showed that the simulation results are a little bit smaller than the H-W solutions but greater than the H-F solutions, since the permeability of shale rock in this part is close to the permeability adopted in the H-W solution. The similar results can be also found in simulations conducted by Lu et al. [48] and Zhang et al. [45]; this verifies the accuracy of the model coupling rock deformation and gas seepage and the numerical implementation, though no bedding planes are involved in this part.

#### 4. Numerical Settings

Shale rock naturally contains bedding planes with different directions owing to geological deposition and folding. These bedding planes in different directions play a significant role in the propagation of the hydraulic fracture. Thus, numerical samples with bedding planes in different directions are adopted in this work. As shown in Figure 6, the numerical sample is a 2D plane rectangle () with a borehole (0.04 m in diameter) in its center. The dotted lines in the numerical sample represent bedding planes of which the thickness has a uniform value of 2 mm. Bedding planes are distributed uniformly in the numerical sample, and the distance between two adjacent bedding planes is 0.02 m. It should be noted that is the angle between the bedding plane and the horizontal direction. As for boundary settings, loads applied on the left boundary and top boundary represent horizontal in situ stress and vertical in situ stress, respectively. And the right boundary and bottom boundary are set as rollers. All boundaries are set as no-flow boundaries except for the borehole, on which the fracturing fluid is injected with an injection rate of 0.0053 m^{3}/s. The parameters used in the simulations can be found in Table 1.

#### 5. Results and Discussion

The physical processes for CO_{2}-driven hydraulic fracture trajectories across a bedding plane are simulated. Besides, the evolution of intersection scenarios between hydraulic fractures and bedding planes under several conditions (stress ratio, bedding plane angle, and bedding plane stiffness) is obtained and discussed.

##### 5.1. Effect of Stress Ratio

###### 5.1.1. Simulation Results

A series of numerical simulation tests were conducted to study the impact of stress ratio on the propagation of hydraulic fractures. Eight stress ratios were used in this study which are (0.2), 2 MPa/5 MPa (0.4), 3 MPa/5 MPa (0.6), 4 MPa/5 MPa (0.8), 5 MPa/4 MPa (1.25), 5 MPa/3 MPa (1.67), 5 MPa/2 MPa (2.5), and 5 MPa/1 MPa (5). Bedding planes were set as horizontal for all numerical samples in this part, and boundary settings and mechanical parameters can be found in Section 4.

In Figure 7, the fracture initiation pressure means the injection pressure at the borehole when damage first occurs in numerical samples. It should be noted that the horizontal in situ stress is kept at 5 MPa while the vertical in situ stress increases from 1 MPa to 5 MPa. It can be seen from the results that the fracture initiation pressure increases from 2.4 MPa to 13.5 MPa as the stress ratio increases from 0.2 to 1. And the relation between the fracture initiation pressure and the stress ratio can be fitted by a linear function.

Figure 8 shows the distribution of hydraulic fractures in samples with horizontal bedding planes under different stress ratios. As for stress ratio , the main hydraulic fracture propagates horizontally, that is, the direction of the maximum principal stress. Also, hydraulic fractures are observed inserting into bedding planes because of their lower strength and elastic modulus compared with the rock matrix. No fracture propagates along the vertical direction. When the stress ratio , hydraulic fractures mainly propagate along the horizontal direction or insert into bedding planes. Some fractures are observed approaching the bedding plane and then branching into it or directly crossing through the bedding plane. As for , the propagation of hydraulic fractures in the horizontal direction is not as deep as the case of . Hydraulic fractures branching into and crossing through the bedding planes can be observed in the vertical direction. In the following cases, vertical stress remains at 5 MPa, while horizontal stress varies from 4 MPa to 1 MPa. For the stress ratio and , the intersection scenarios of hydraulic fractures and bedding planes are similar with the case of , containing hydraulic fractures crossing through bedding planes, approaching bedding planes, and then branching into and crossing through bedding planes. However, the propagation of hydraulic fractures in the horizontal direction is gradually restricted. As for and , a main vertical hydraulic fracture is formed along the direction of the maximum principal stress and almost no horizontal hydraulic fractures are observed in this condition. Based on the above results, it is indicated that the bedding plane and stress ratio have a significant influence on the distribution of hydraulic fractures. The horizontal radius and vertical radius of hydraulic fractures under different stress ratios is shown in Figure 9, which can reflect the fracturing area in a specific case. When the stress ratio increases from 0.2 to 5, the horizontal radius decreases first dramatically and then slowly from 0.08 m to 0.007 m, whilst the vertical radius increases first sharply and then gradually from 0.00243 m to 0.0622 m.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

###### 5.1.2. Intersection Scenario between the Hydraulic Fracture and Bedding Plane

Understanding the complexity of the fracture network is of vital importance for hydraulic fracturing design. The fracture network near the gas reservoirs is formed through propagation and combination of basic intersection scenarios. Thus, the basic intersection types between the hydraulic fracture and bedding plane are summarized in this subsection.

Based on the simulation results above, four types of intersection scenarios between hydraulic fractures and bedding planes are shown in Figure 10. (a) Inserting: the hydraulic fracture inserts into a bedding plane and continues to propagate along it, i.e., the hydraulic fracture is arrested by the bedding plane. (b) L-shaped crossing: the hydraulic fracture approaches the fracture/bedding plane then branches into the plane without crossing it. (c) T-shaped crossing: the hydraulic fracture approaches the fracture/bedding plane, branches into it, and crosses through it. (d) Direct crossing: the hydraulic fracture crosses one or more bedding planes without branching into them. It indicates that the intersection types vary from (a) → (b) → (c) → (d) with the increase of stress ratio in specimens with the horizontal bedding planes.

**(a)**

**(b)**

**(c)**

**(d)**

##### 5.2. Effect of Bedding Plane Angle

Numerical simulations were performed to study the effect of the bedding plane angle on the propagation of hydraulic fractures. Six kinds of specimens with different bedding plane angles (, , , , , and ) were used in this part. The horizontal in situ stress and vertical in situ stress were kept at 1 MPa and 5 MPa, respectively. And other boundary settings and mechanical parameters can be found in Section 4.

The distribution of hydraulic fractures in specimens with different bedding plane angles is shown in Figure 11. For the bedding plane angle , a vertical main fracture is formed crossing through bedding planes since the vertical in situ stress is five times of the horizontal in situ stress. When , the distribution of hydraulic fractures is similar with the former case (). The results indicate that the bedding planes with angle do not have an obvious influence on the propagation of hydraulic fractures. With the bedding plane angle increasing to 30°, the intersection type between the hydraulic fracture and bedding plane changes from direct crossing to T-shaped crossing. It should be noted that only a small amount of hydraulic fractures branch into the bedding planes. As for the bedding angle , the T-shaped crossing is the primary intersection scenario; more hydraulic fractures branch into the bedding planes compared with the case of . With the bedding plane angle increasing to 60°, the results show that the intersection scenario changes from T-shaped crossing to inserting; no crossing scenario is observed in this condition. When the bedding plane angle , hydraulic fractures propagate along the vertical direction or insert into bedding planes. The results indicate that the intersection type changes from (d) → (c) → (a) with the increase of the bedding plane angle .

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

##### 5.3. Effect of Bedding Plane Stiffness

Three numerical tests were performed to research the behavior of hydraulic fractures propagating in specimens with different levels of stiffness of bedding planes. The stiffness (elastic modulus ) of the rock matrix was fixed whilst three different levels of stiffness of bedding planes (elastic modulus ) were used in this part (, , and ). The bedding plane angle was 30° in the specimen, and the horizontal in situ stress and vertical in situ stress were 1 MPa and 5 MPa, respectively. Again, the boundary settings and the injection rate are the same as those defined in Section 4.

Figure 12 shows the distribution of hydraulic fractures in specimens with different levels of stiffness of bedding planes. For the case of , when the hydraulic fracture approaches the bedding plane, it is easy to branch into the bedding plane since the stiffness of the bedding plane is much smaller than that of the rock matrix. T-shaped crossing is formed between hydraulic fractures and bedding planes. Comparing the case of with the former one (), it can be found that the number of hydraulic fractures branching into bedding planes decreases. T-shaped crossing is also observed in this case. For the case of , a vertical main fracture is formed crossing bedding planes, and direct crossing is the intersection scenario between the hydraulic fracture and bedding plane. The results show that the intersection type changes from T-shaped crossing to direct crossing with the increase of the bedding plane stiffness. Besides, the bedding plane stiffness has an obvious impact on the development of the seepage area and acoustic emission (Figure 13). It should be noted that the acoustic emission is recorded from the number of damaged REVs per second. In all three cases, the seepage area increases slowly in the initiate stage (1 s–30 s) then increases rapidly. For the case of , the seepage area is 0.00418 m^{2} at time , and the highest acoustic emission is 655. When the stiffness ratio increases from 0.1 to 0.4 and 0.8, the seepage area decreases 22.2% and 41.8% to 0.00325 m^{2} and 0.00243 m^{2}, respectively, whilst the highest acoustic emission decreases to 550 and 364, respectively. The results indicate that hydraulic fractures are formed more easily, and a relatively more complex intersection scenario is obtained with a lower stiffness of the bedding plane.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

#### 6. Conclusions

Understanding the mechanism of the intersection scenario between hydraulic fractures and bedding planes is of vital importance for creating a complex fracture network and improving the recovery of unconventional resources. In this work, a coupled hydraulic-mechanical model is developed where a damage evolution law is used to describe the initiation and propagation of hydraulic fractures. This model is then used to conduct a series of numerical simulations to investigate the propagation of hydraulic fractures in specimens with bedding planes. The following conclusions can be obtained: (1)Stress ratio has a vital impact on the intersection scenario between a hydraulic fracture and a bedding plane. Four types of intersection scenarios are summarized based on the study: (a) inserting—the hydraulic fracture inserts into a bedding plane and continues to propagate along it; (b) L-shaped crossing—the hydraulic fracture approaches the fracture/bedding plane then branches into the plane without crossing it; (c) T-shaped crossing—the hydraulic fracture approaches the fracture/bedding plane, branches into it, and crosses through it; (d) direct crossing—the hydraulic fracture crosses one or more bedding planes without branching into them. The simulation results indicate that intersection types vary from (a) → (b) → (c) → (d) in specimens with horizontal bedding planes when the stress ratio increases from 0.2 to 5. Besides, the fracture initiation pressure increases from 2.4 MPa to 13.5 MPa while the stress ratio increases from 0.2 to 1(2)The bedding plane angle can also greatly affect the propagation of hydraulic fractures. The results show that the intersection type changes from (d) → (c) → (a) with the increase of the bedding plane angle (0° → 90°)(3)We also investigate the influence of bedding plane stiffness on the propagation of hydraulic fractures. The results indicate that the intersection type changes from the T-shaped crossing to the direct crossing with the increase of the bedding plane stiffness. When the stiffness ratio increases from 0.1 to 0.4 and 0.8, the seepage area decreases 22.2% and 41.8% and the highest acoustic emission decreases from 655 to 550 and 364, respectively

#### Data Availability

The numerical simulation data used to support the findings of this study are available from the corresponding author upon reasonable request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China (51804339).