Numerical Investigation on the Hydraulic Fracture Evolution of Jointed Shale

Joints are a common structure of heterogeneous shale rock masses, and in situ stress is the environment in which heterogeneous rock masses can be found.-e existence of joint plane and confining pressure difference influences the physical properties of shale and propagation of fractures. In this study, jointed shale specimens were simulated under different confining pressures to explore the failure patterns and fracture propagation behavior of hydraulic fracturing. Different from the common research of hydraulic fracturing on signal parallel joint rock mass, the simulations in this study considered three points (parallel joint, multi-dip angle joint, and no-joint points). -e effects of the single-dip angle joint, multi-dip angle joint, and confining pressure difference on the hydraulic fracture evolution and stress evolution of the jointed shale were studied comprehensively. -e confining pressure difference coefficient proposed in this study was used to accurately describe the confining pressure difference. Results indicate that the larger is the confining pressure difference, the stronger is the control of the maximum principal stress on fracture evolution; by contrast, the smaller is the confining pressure difference, the stronger is the control of the joint plane on fracture evolution. Under the same confining pressure difference, the hydraulic fracture propagates more easily along the small dip angle joint plane. As the value of the confining pressure difference coefficient moves closer to zero, the hydraulic fracture propagates randomly, the tensile stress region around the fracture tip widens, and the joint planes fractured by tensile increase. -is study can offer valuable guidance to the design of unconventional reservoir reconstruction.


Introduction
Shale is a kind of unconventional reservoir with strong heterogeneity and low permeability [1,2]. As a common structure in shale reservoirs, joints are expected to affect the mechanical properties of the shale rock mass, and the hydraulic fracture evolution in the shale may differ as opposed to those in other conventional rock masses. In situ stress, a type of long-term stress of shale reservoirs, is a function of time and space. Different in situ stresses affect the mechanical properties of rock masses during diagenesis and the initiation and propagation of fractures during fracturing. Hydraulic fracturing is an important consideration in ensuring production efficiency, a measure that is commonly used in petroleum engineering [2,3]. Unconventional oil and gas production may eventually replace depleted conventional oil and gas resources; thus, unconventional reservoir reconstruction, especially shale reservoir reconstruction, has become one of the popular research topics in recent years [4][5][6]. Knowing how a complex fracture network forms and how this network can improve reservoir permeability is an important research component of unconventional reservoir reconstruction [7,8]. Furthermore, discontinuities dominate the geometry, deformation modulus, strength, failure behavior, and permeability of rocks, and the existence of joints promotes the formation of complex fracture networks [9]. erefore, the study regarding the growth process of hydraulic fractures of jointed shale under different confining pressures is essential, as this topic is one of the important aspects of shale reservoir reconstruction research. Its investigation can also help to determine efficient fracturing scenarios and production designs in field operations.
Many scholars have recently investigated the hydraulic fracture evolution of jointed rocks and established that they can positively contribute to unconventional reservoir reconstruction. Gong et al. [10], who theoretically analyzed and adopted the flow-stress-damage coupled approach to investigate the propagation of hydraulic fractures under the joint impact of bedding planes and natural fractures, found that the main influencing factors of fracture propagation are the differences in horizontal principal stress. Chung et al. [11], who used the discrete element method to evaluate the influence of sedimentary beddings and preexisting joints on the fracturing response, found that the joint characteristics (i.e., orientation, aperture, and healing conditions) could affect the behavior of shale formation. Heng et al. [12] conducted threepoint bending tests of notched cylindrical specimens with different bedding orientations to investigate the influence of bedding planes on hydraulic fracture evolution. Men et al. [9] numerically simulated the failure patterns and fracture propagation behavior of jointed rocks in hydraulic fracturing and found that the maximum principal stress and joint plane are the factors that control fracture propagation in the global and in local scales, respectively.
In this study, a simulation investigation of the influence of single-dip angle joint, multi-dip angle joint, and confining pressure difference coefficient on the hydraulic fracture evolution and stress evolution of jointed shale was conducted using the rock failure process analysis (RFPA) system developed by Professor C.A. Tang in 1995. As a numerical simulation code, the RFPA code can simulate the progressive failure process of materials, and it has been widely used and recognized by scholars worldwide [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28]. e simulation results of jointed shale on hydraulic fracturing can also provide relevant reference about the formation of complex fracture networks in jointed shale.

Brief Introduction of the RFPA-Flow Code
Rock is a typical heterogeneous material with randomly distributed defects and mesodamage [29], and most codes cannot simulate this characteristic of the rock. Developed for heterogeneous materials (e.g., rock or concrete), the RFPA code is based on finite element and statistical damage theory, and it can simulate the failure process of heterogeneous and permeable geomaterials. In the simulation of heterogeneity and random distribution of material defects, the mesoscopic elements are assumed to be isotropic and homogeneous, and their mechanical properties (e.g., Young's modulus, Poisson's ratio, and strength properties, among others) are assumed to be linear in the RFPA model ( Figure 1) [30]. ese linear mesoscopic elements are also statistically distributed (e.g., normal, Poisson's, and Weibull distributions). Weibull distribution is commonly used to describe the strength of brittle materials (e.g., rock) [31]. However, the mechanical properties of the macromodel are nonlinear (Figure 1). In RFPA, fournode isoparametric elements are used as the basic elements. Moreover, in simulations, mesoscopic elements are defined as damaged (the damage element is shown in black) if they satisfy the strength criterion (e.g., double-shear, Hoek-Brown, or Mohr-Coulomb criterion). As stress increases, more and more elements ultimately fail. When these elements connect with each another, macrofractures are formed, and a failure process of the materials is reached (Figure 2).

Fundamental Equations in the RFPA-Flow
Code. In the RFPA model, φ(a) is defined as the statistical distribution density of the mechanical parameter of the mesoscopic elements (equation (1)) [32] and used to describe the heterogeneity of materials. e variable a represents the mechanical parameters (e.g., Young's modulus, Poisson's ratio, tensional strength, and compressive strength) of the elements, a 0 represents the scale parameters related to the average values of the mechanical parameters, and m is the heterogeneity index used to describe the heterogeneity of the solid materials [33,34]. A high m value indicates a concentrated distribution of the element's mechanical properties (i.e., a highly homogeneous material), whereas a low m value means dispersed distribution (i.e., a highly heterogeneous material) ( Figure 3). For heterogeneous rock materials, the m value is generally between 1 and 5. Figure 3 shows the discrete degree of mechanical parameters of the mesoelements in the RFPA model when the m value is 1 to 5 in the Weibull distribution.
e geological medium (rock) in the RFPA-Flow code is assumed to be fully saturated with fluid and that the fluid flow is governed by Darcy's law. e code can also be used to explore the effects of stress and damage on the permeability of the heterogeneous medium. e coupled process of stress and seepage in the deforming rock mass is governed by Biot's theory of consolidation. e basic formulations of the fluid-solid coupling in the RFPA code that are used in the analysis are as follows [27].
Equilibrium equation: Geometrical equation: Constitutive equation: Seepage equation: where ρ is the volume density; σ ij is the total stress; ε ij is the total strain and ε v is the volume strain; U i is the displacement; p is pore pressure; λ is the Lame coefficient; δ ij is the Kronecker constant; G is the modulus of shear deformation; Q is Biot's constant; K is the permeability coefficient; and α is the coefficient of pore pressure. e permeability of the heterogeneous material changes continuously in the failure process of the heterogeneous medium.
is condition can be attributed to the varying fluid flow behaviors between that in the propagating fractures and that in the preexisting fractures. e initiation of new fractures, propagation of preexisting fractures, coupling relationship of flow-stress-damage in the two kinds of fractures, and variations of rock permeability in the failure process can therefore also be considered in the RFPA code.
In elastic state, the relationship of the stress and permeability coefficient is described by where K 0 is the initial permeability coefficient; σ ′ is the effective stress; and b is the coupling parameter. As damage occurs, the permeability coefficient of mesoelement will change. To describe the relationship of permeability on the stress and the damage, the following coupling equation is introduced: where ξ(ξ ≥ 1.0) is a mutation coefficient of permeability accounting for the increase in permeability once the element reaches the damage state and β is the coupling parameter reflecting the influence of stress on the permeability coefficient. e values of the coefficients ξ, α, and β are determined experimentally, and they vary with the stress state and damage evolution of the heterogeneous medium. In the RFPA code, 0 ≤ α < 1 and ξ � 1 are assumed for the mesoscopic element in the elastic stage. Once damage occurs, the permeability of the mesoscopic element increases significantly. Furthermore, ξ � 5 is assumed for the damaged element and α � 1, ξ � 100 is assumed for the fully damaged element.

Establishment of RFPA Model.
In the RFPA code, a rectangular model (i.e., an irregular model can be imported by another code) is built and divided into plane four-node isoparametric elements. In this manner, the accuracy requirements can be ensured according to the computer's configuration. e material parameters of the matrix are given in the established rectangular matrix model. e establishment of the rock formation, coal seam, joint, fracture, and hole on the matrix is programmed by drawing different shapes (e.g., rectangular, line, circular, or arbitrary shapes) on the matrix and assigning the corresponding mechanical properties to the elements according to a prescribed value range of the shape. In this study, the joint is established by drawing lines and assigning the joint properties to the elements on the lines.

Numerical Model Set
e simulations were divided into three groups as a means of investigating the influence of joint plane and maximum principal stress on hydraulic fracture evolution ( Figure 4). Two of the samples were jointed specimens, whereas one of them was a common shale specimen. e first group entailed a hydraulic fracturing simulation of the parallel jointed shale under varying confining pressure differences (σ H �10 MPa, σ V � 5 MPa; σ H � 10 MPa, σ V � 8 MPa; and σ H � 10 MPa, σ V � 10 MPa). As shown in Figure 4(a), the joint dip angle α is the angle between the joint plane and the maximum principal stress direction (horizontal direction). e α values of the seven models in the first group were 0°, 15°, 30°, 45°, 60°, 75°, and 90°. e second group involved the simulation of the shale specimen without joints under varying confining pressure differences Figure 4(b)). Finally, as the study also covered the influence of joint and maximum principal stress on fracture evolution in the presence of multi-dip angle joints, the third group entailed the simulation of multi-dip angle joint models subjected to hydraulic fracturing under varying confining differences (Figure 4(c)). e model size in the simulations was 540 mm × 540 mm, and a wellbore with a radius of 25 mm was assumed to be located at the center of the specimen. e model was discretized into 72,900 elements. e thickness of the joint was 2 mm (single element), and the distance between the adjacent joint planes was 60 mm. erefore, the joint was established by assigning a given set of joint properties to the elements on lines. An increasing water pressure of 4 MPa for the initial value, followed by 0.1 MPa increments per step, was inputted into the wellbore. A planestrain component was employed for calculation. e list of rock and joint parameters used in the simulations is shown in Table 1 [35].

Numerical Results of the Parallel Joint Models.
e parallel joint simulations of the seven joint dip angles with different confining pressures were divided into three parts. e 21-specimen composite with varying joint dip angles and confining pressure differences was simulated to investigate the influences of hydraulic fracture evolution. e findings showed that the initiation and propagation of fractures significantly changed when the joint dip angle increased and the confining pressure difference decreased (Figures 5-7). e models with the large confining pressure difference are shown in Figure 5. When the joint dip angle was between 0°to 30°, the hydraulic fracture was initiated and then it propagated initially in the maximum principal stress direction, and then it propagated further along the joint after the fracture was connected to the joint. Under this scenario, the weak joint plane controlled the fracture evolution. When the joint dip angle was increased to 45°, the hydraulic fracture propagated partly along and partly across the joint plane, and the fracture propagated in the maximum principal stress direction macroscopically. erefore, the maximum principal stress and the weak joint plane both controlled the fracture evolution. When the joint dip angle was further increased (60°to 90°), the joint plane had no effect on fracture evolution. e hydraulic fracture was initiated and then propagated in the maximum principal stress direction. e maximum principal stress controlled the fracture evolution.
As shown in Figure 6, the fracture evolution of the same seven parallel joint models eventually changes under the small confining pressure difference. e influence scope of the joint in this group was larger than that of the specimens with a large confining pressure difference. When the joint dip angle was between 0°and 45°, the hydraulic fracture was initiated in the maximum principal stress direction and then propagated in a straight path along the joint plane. When the joint dip angle was increased to 60°, the fracture propagated partly across and partly along the joint. When the joint dip angle was further increased (75°to 90°), the hydraulic fracture was initiated and then propagated at a path close to the horizontal direction, and the maximum principal stress essentially controlled the fracture evolution. e last two simulations (dip angles of 75°and 90°) can be further compared using Figures 5 and 6. e fracture in Figure 6 is not as horizontal as that in Figure 5, which means that the influence of maximum principal stress on fracture evolution is decreased as the confining pressure difference decrease. Figure 7 shows the simulation results of the same seven parallel joint models under the condition of no-confining pressure difference. In the absence of maximum principal stress, the joint controlled the fracture evolution in all of the joint dip angle models. When the joint dip angle was between 0°and 90°, the fracture propagated along the joint plane near the wellbore, and the joint plane fully controlled the fracture evolution.

Numerical Results of the No-Joint Models.
ree shale specimens without joint planes under different confining pressures were simulated to investigate the influence of confining pressure difference on fracture evolution. e fracture dip angle changed when the confining pressure difference decreased (Figure 8). When the model was simulated under a larger confining pressure difference, the fracture was initiated and then propagated in the maximum stress direction. Here, the maximum principal stress controlled the fracture evolution, and the jagged fracture extended in the horizontal direction in the global scale (Figure 8(a)). When the confining pressure difference was small, the fracture deviated by a small angle from the horizontal direction, and the influence of maximum principal stress decreased when the confining pressure difference decreased (Figure 8(b)). As for the model with no-joint plane and confining pressure difference, the fracture was initiated and then propagated randomly (Figure 8(c)).
Hydraulic fractures are generally formed by tension. Here, the tensile stress (green region) developed around the fracture tip during fracture propagation and eventually pulled apart more rocks because of the increase in water pressure. When the confining pressure difference decreased, the tensile stress region widened, and the secondary fractures increased in number (Figures 8(d)-8(f )). e confining pressure difference coefficient μ was used to accurately describe the confining pressure difference under varying confining pressures.
where μ is the confining pressure difference coefficient (0 ≥ μ ≤ 1), σ H is the horizontal in situ stress, and σ V is the vertical in situ stress. e closer is μ to 0, the smaller is the  Figure 5: Pore pressure field of models with large confining pressure difference (σ H �10 MPa,

Advances in Materials Science and Engineering
confining pressure difference; the closer is μ to 1, the larger is the confining pressure difference. e confining pressure coefficient values of the three no-joint models as computed by equation (2) were 0.5, 0.2, and 0.
As shown in Figure 9, the breakdown pressure linearly decreases as the confining pressure difference coefficient increases. e larger is the coefficient μ, the easier the specimen becomes fractured, and the greater is the influence of the maximum principal stress.
As shown in Figures 10(a) and 10(b), the fractures propagate along the 0°and 15°joint plane, which is the area with the smallest dip angle joint in each side of the well. en, as shown in Figure 10(c), the fracture propagates along the 75°and 90°joint plane as the confining pressure difference continues to decrease to 0, and the change in fracture evolution becomes apparent. e absence of a confining pressure difference indicates the absence of both maximum principal stress and dip angle. e fracture randomly propagates along the joints, as depicted in Figure 10(c), because the effect of all joints on fracture evolution is the same. In summary, hydraulic fractures will propagate first along the smallest dip angle joint in the presence of a confining pressure difference, whereas such fractures will propagate randomly along the joint plane in the absence of a confining pressure difference.
In the figures, the minimum principal stress represents the condition in which fractures occur because of tensile stress (green region) during hydraulic fracturing, the confining pressure difference is small, and the tensile region is large. Here, the fracture propagates along the 0°and 15°joint plane in both Figures 10(d) and 10(e), but more failure elements can be observed on the 30°joint in Figure 10(e) than those in Figure 10 Figure 10(f ). is finding can be attributed to the 30°joint plane that is closer to the tensile region in Figure 10(e) than that in Figure 10(d), and two fractured joints can be observed in the tensile region in Figure 10(f ). e joint in the tensile region will be fractured more easily.

Comparison of Numerical and Experimental Results
A numerical jointed shale model is established according to an experiment to prove the reliability of the RFPA-Flow code [36]. e dimension of the specimen is 300 mm × 300 mm, and the radius of the well is 12 mm. e model is discretized into 300 × 300 elements. e confining pressure values of σ H and σ h are 5 and 3 MPa, respectively. e parameters used in the simulation are shown in Table 2. Moreover, the planestrain calculation is used in the simulation. Figure 11 presents a comparison of the fracture evolution between the numerical and experimental results. e confining pressure difference coefficient of the model is 0.4 according to equation (8). Furthermore, as previously discussed in Section 4, the maximum principal stress and joint plane can both influence the fracture evolution. e numerical results are shown in Figure 11. e fracture is initiated and propagates in the maximum principal stress direction, and then it propagates along the joint when the main fracture connects with the joint plane. e plane-strain calculation is used in the simulation, but the main fracture in the numerical results depicts it to be neither positioned across the joint nor forms a fracture through the model. is condition differs from those in the experimental 9    results. e difference can be attributed to the limitation of 2D simulation. However, the distribution of tensile stress (green area in the minimum principal stress photograph in Figure 11) suggests that the main hydraulic fracture will continue to propagate across the joint plane in the maximum principal stress direction, but the premise is the model is stable. Overall, the fracture evolution is essentially consistent between the simulated and experimental results, and they further verify the findings of the previous study described in Section 4 and the reliability of the RFPA-Flow code.

Conclusion
(1) In situations in which the joint plane and the maximum principal stress both exist, the joint plane controls the fracture evolution when the joint dip angle is small, while the maximum principal stress controls the fracture evolution when the joint dip angle is large. (2) e smaller is the confining pressure difference, the lower is the influence of the maximum principal stress on the fracture evolution, and the greater is the influence of the joint on the fracture evolution. (3) When the confining pressure difference coefficient is 0, the fracture propagates randomly along the joint, and the joint fully controls the fracture initiation and propagation. When the confining pressure difference coefficient close to 1 (the maximum value is 0.5 in this study), the fracture propagates in the horizontal     Advances in Materials Science and Engineering direction, and the maximum principal stress fully controls the fracture initiation and propagation. (4) e smaller is the confining pressure difference, the larger is the tensile stress region around the fracture tip, and the more joints and secondary fractures will be fractured by tensile stress during hydraulic fracturing.

Data Availability
e data used to support the findings of this study can be provided by the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.