CO2-Driven Hydraulic Fracturing Trajectories across a Preexisting Fracture

Beijing Research Institute of Uranium Geology (BRIUG), 10, Xiao-Guan-Dong-Li, P.O. Box 9818, Beijing 100029, China CNNC Key Laboratory on Geological Disposal of High-Level Radioactive Waste, 10, Xiao-Guan-Dong-Li, P.O. Box 9818, Beijing 100029, China State Key Laboratory for Geomechanics & Deep Underground Engineering, China University of Mining & Technology, Xuzhou, 221116 Jiangsu, China Department of Energy and Mineral Engineering, EMS Energy Institute and G3 Center, Pennsylvania State University, University Park PA 16802, USA


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][2][3][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][12][13][14][15]. However, waterbased 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][26][27][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][35][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.

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.

Rock Deformation Equation.
According to the elastic theory, the constitutive equation considering the influence of pore pressure can be written as where σ ij is the stress tensor, G is the shear modulus of rock, ε ij is the strain tensor, ν is Poisson's ratio of rock, ε v = ε 11 + ε 22 + ε 33 is the volumetric strain, δ ij is the Kronecker delta, η is the Biot coefficient, and p 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 u is the displacement vector and f i is the component of the body force.

Gas Flow Equation.
The governing equation for CO 2 flow based on mass balance can be defined as where m is the CO 2 mass per volume of rock, ρ g is the CO 2 density, q g is the seepage velocity of CO 2 , Q m is the source of CO 2 , and t is the time variable. According to Darcy's law, the seepage velocity of CO 2 can be defined as where k is the permeability of the shale rock and μ g 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 m = ρ g ϕ (ϕ 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 c = ð1/ρ g Þð∂ρ g /∂pÞ 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 σ 1 and σ 3 are the first principal stress and third principal stress, respectively; f t0 and f c0 are the tensile strength and compressive strength of the REV, respectively; φ is the internal friction angle; and H 1 and H 2 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 D is used to represent the damage level and is defined as [43] where ε 1 and ε 3 are the first principal stain and the third principal strain, respectively, and ε t and ε c 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 D can be defined as follows: where E and E 0 are the elastic modulus and the initial elastic modulus of a REV, respectively; k 0 is the initial permeability of a REV; and α k is the damage-permeability effect coefficient to define the influence of damage on permeability.

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 z represents the mechanical parameter of REVs, z 0 is the average value of the mechanical parameter, and ξ is the heterogeneity coefficient.

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.  Check whether damage occurs in REVs by equations (7) and (8)?
Update the parameters Check whether the damage area is extended?

Add a load increment
Check whether there is sample failure?      Geofluids (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 P HW and P HF are the breakdown pressures of the H-W solution and the H-F solution, respectively; σ t is the tensile strength of the shale rock; and p 0 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 σ y is kept at 30 MPa, while the horizontal in situ stress varies from 10 MPa to 30 MPa. And a tectonic stress ratio β (β = σ y /σ x ) 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.

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 (0:2 m × 0:2 m) 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

Geofluids
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.

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). 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 β = 0:2, the main hydraulic fracture propagates horizontally, that is, the direction of the maximum principal stress. Also, hydraulic fractures are observed However, the propagation of hydraulic fractures in the horizontal direction is gradually restricted. As for β = 2:5 and β = 5, 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.

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.

Geofluids
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. 8 Geofluids 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 α = 0°, 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 α = 15°, the distribution of hydraulic fractures is similar with the former case (α = 0°). The results indicate that the bedding planes with angle α = 15°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 α = 45°, the T-shaped crossing is the primary intersection scenario; more hydraulic fractures branch into the bedding planes compared with the case of α = 30°. With the bedding plane angle increasing to 60°, the results show that the intersection scenario changes from Tshaped crossing to inserting; no crossing scenario is observed in this condition. When the bedding plane angle α = 90°, 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 α.

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 E 2 ) of the rock matrix was fixed whilst three different levels of stiffness of bedding planes (elastic modulus E 1 ) were used in this part (E 1 /E 2 = 1/10, E 1 /E 2 = 4/10, and E 1 /E 2 = 8/10). 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 E 1 /E 2 = 1/10, when the hydraulic fracture approaches the bedding plane, it is easy to branch into the 9 Geofluids 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 E 1 /E 2 = 4/10 with the former one (E 1 /E 2 = 1/10), 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 E 1 /E 2 = 8/10, 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 Tshaped 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 E 1 /E 2 = 1/10, the seepage area is 0.00418 m 2 at time t = 40 s, and the highest acoustic emission is 655. When the stiffness ratio E 1 /E 2 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.

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  Geofluids 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 E 1 /E 2 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.