# Effect of Joint Geometrical Parameters on Hydraulic Fracture Network Propagation in Naturally Jointed Shale Reservoirs

Authors: Zhaohui Chong; Qiangling Yao; Xuehua Li
Journal: Geofluids
 (2018)

---

## Abstract

The presence of a significant amount of discontinuous joints results in the inhomogeneous nature of the shale reservoirs. The geometrical parameters of these joints exert effects on the propagation of a hydraulic fracture network in the hydraulic fracturing process. Therefore, mechanisms of fluid injection-induced fracture initiation and propagation in jointed reservoirs should be well understood to unleash the full potential of hydraulic fracturing. In this paper, a coupled hydromechanical model based on the discrete element method is developed to explore the effect of the geometrical parameters of the joints on the breakdown pressure, the number and proportion of hydraulic fractures, and the hydraulic fracture network pattern generated in shale reservoirs. The microparameters of the matrix and joint used in the shale reservoir model are calibrated through the physical experiment. The hydraulic parameters used in the model are validated through comparing the breakdown pressure derived from numerical modeling against that calculated from the theoretical equation. Sensitivity analysis is performed on the geometrical parameters of the joints. Results demonstrate that the HFN pattern resulting from hydraulic fracturing can be roughly divided into four types, i.e., crossing mode, tip-to-tip mode, step path mode, and opening mode. Asβ (joint orientation with respect to horizontal principal stress in plane) increases from 0° to 15° or 30°, the hydraulic fracture network pattern changes from tip-to-tip mode to crossing mode, followed by a gradual decrease in the breakdown pressure and the number of cracks. In this case, the hydraulic fracture network pattern is controlled by both γ (joint step angle) and β. When β is 45° or 60°, the crossing mode gains dominance, and the breakdown pressure and the number of cracks reach the lowest level. In this case, the HFN pattern is essentially dependent on β and d (joint spacing). As β reaches 75° or 90°, the step path mode is ubiquitous in all shale reservoirs, and the breakdown pressure and the number of the cracks both increase. In this case, β has a direct effect on the HFN pattern. In shale reservoirs with the same β, either decrease in k (joint persistency) and e (joint aperture) or increase in d leads to the increase in the breakdown pressure and the number of cracks. It is also found that changes in d and e result in the variation in the proportion of different types of hydraulic fractures. The opening mode of the hydraulic fracture network pattern is observed when e increases to 1.2 × 10−2 m.


---

## Body

## 1. Introduction

Hydraulic fracturing is widely used in the petroleum industry for boosting the yield of oil/gas wells [1–3]. It was later introduced to the mining industry to weaken the hard roof so as to avoid rock burst hazard [4]. Proved by research results in recent years, hydraulic fracturing is an efficient stimulation approach that can significantly enhance the permeability [5–7]. In the past studies, shale reservoirs are regarded as dual-porosity medium composed of shale matrix and fracture (joint) network [8, 9]. Through hydraulic fracturing, the high-pressure water is able to penetrate the shale matrix and joint, thus not only releasing the adsorbed gas but also creating a fracture network for gas to transport [10, 11]. However, the uncertainty regarding the composition of the shale reservoirs compounded by the complex geological setting makes it extremely difficult to control and predict fracture initiation and propagation in hydraulic fracturing [12]. Therefore, it is of great significance to gain a good understanding of fracture initiation and crack propagation, which can guide hydraulic fracturing operation in naturally jointed shale reservoirs.A series of models were established by the predecessors to investigate the correlation between the hydraulic fracture (HF) and the naturally occurring discontinuities such as natural fracture (NF), joint, fissure, and cleat. Researchers and scholars also put forward some criterions regarding the relation between the HF and NF. Renshaw and Pollard [13] proposed a relatively simple criterion to determine the propagation direction of the HF after it traverses through the NF that is normal to it. Also, this criterion was proved effective by a following laboratory experiment. Gu and Weng [14] and Gu et al. [15] developed this criterion and employed an updated criterion to investigate the correlation between the nonorthogonal HF and NF. To date, the correlation between HF and NF has been divided mainly into 6 types [13, 16, 17]: (1) HF directly crossing NF, (2) HF intersecting NF, (3) HF crossing NF with an offset, (4) HF arrested by NF, (5) HF shear slipping along NF, and (6) HF branching or turning at the end of NF.To simulate HF initiation and propagation in jointed rock mass, scholars made huge efforts. Xu et al. [18], Meyer and Bazan [19] put forth a wire-mesh model and discrete fracture network (DFN) model. In these pseudo-3D models, the hydraulic fracture network (HFN) is made up of two sets of cracks. One is the cracks that extend along the direction of the principal stress, and the other is the cracks that are in parallel with the direction of the principal stress by a fixed interval. Weng et al. [20] put forward an unconventional pseudo-3D model that considers several important characteristics of fracturing treatment simulation such as transport of proppant, behavior of non-Newtonian fluid, and permeability parameters of the wellbore. More importantly, the proposed model is able to simulate HFN resulting from the interaction of HF and joints. Olson and Taleghani [21] presented a 2D plane-strain extended finite element method (XFEM), which can model the formation of fracture wing and makes it possible for HF to change the propagation direction. The boundary element method [22] is also introduced to study the interaction between the joints and NF by analyzing the tensile and shear cracks. All the above methods are continuum-based, which can generate arbitrary fracturing paths without remeshing.Recent years have witnessed the popularity of the discrete element method (DEM) in the study of HF propagation. In this method, the reservoir is generated by a great number of assembled bonded particles (in particle flow code (PFC )) [23] or by a great number of transformable blocks (in universal distinct element code (UDEC) and three-dimensional distinct element code (3DEC)) [24, 25]. The joints between the adjacent particles or blocks morphed into the fluid network which provides conditions for HF propagation [26, 27]. The advantage of DEM over the continuum-based method is that a great number of preexisting discontinuities can be added into the modeled reservoir and the mechanical behavior of these discontinuities can be modeled in the formation of the HFN in jointed rock mass.In actuality, shale reservoirs contain a great amount of joints, exhibiting typical inhomogeneous behavior. Figure1 is a site photo showing the Longmaxi shale reservoir in Chongqing, China, which contains a massive amount of joints. These joints occur in a laminar manner. The effect of these joints on the formation of HFN has already been captured by physical experiments [28]. In addition, the activation of these joints induced by HF initiation and propagation has been extensively studied [29, 30]. However, many challenges still lie ahead in the study of the interaction between HF and joints as well as the formation of HFN in naturally jointed shale reservoirs. In the process of performing horizontal well fracturing in shale reservoirs, a massive amount of HFs interact with the joints, damaging the integrity of the shale reservoirs further. In this case, the geometrical parameters of the joints have a pronounced effect on HFN propagation.Figure 1
(a) A shale outcrop exposed area and (b) the enlarged view of joint planes in the Longmaxi reservoirs.In this study, a DEM-based coupled hydromechanical model is developed to explore the effect of the geometrical parameters of the joints on HF initiation and propagation as well as the HFN pattern. The parameters of the matrix and joint used in the shale reservoir model are calibrated through physical experiment, and the hydraulic parameters used in the model are validated through comparing the breakdown pressure derived from numerical modeling against that calculated from the theoretical equation. In addition, the calibrated model is also employed to perform sensitivity analysis on the geometrical parameters of the joints such as joint orientation angleβ, joint persistency k, joint spacing d, joint step angle γ, and joint aperture e.


## 2. Simulation Methodology

### 2.1. Parallel Bonded Model

The PFC, a simulation code based on the DEM is an effective numerical tool for simulating the behavior of jointed rock mass. In PFC, the rock matrix is modeled as an assembly of circular particles. Though PFC is based on DEM, it can still be employed to model the deformation characteristics and mechanical behavior of a continuous medium. The parallel bonded model (PBM) is used to define and simulate the particle contacts as it can deliver more realistic results in modeling rock fracture.Figure2(a) schematizes the micromechanical behavior of the particle contacts. At particle contacts, two sets of springs provide the normal and shear stiffness [31]. As shown in Figures 2(b) and 2(c), the first set consists of two springs that stand for normal and shear stiffness, which are denoted as kn and ks, respectively. In the model, the particles are idealized as a rigid body that is unbreakable and undeformable. However, the particles can overlap if the compressive force is great enough. The second set of springs provides the parallel bond normal and shear stiffness denoted as k¯n and k¯s, which are uniformly distributed over a disc-shaped cross section lying on the contact plane. In addition to the normal and shear stiffness, the parallel bonds also have normal and shear strength. A parallel bond is active over a finite area and can resist the rotation of particles. When the normal or shear stress exceeds the bond strength, the parallel bond is broken and bond stiffness is removed. Once the parallel bond is broken under shear stress, the shear strength falls to its residual value whereas upon breakage of a bond under tensile stress, the tensile strength is set to zero. The residual shear strength depends on the normal force and the friction coefficient of particles at contacts [32]. Microcracks are produced in PBM once the bonds fail due to external loading, and macrocracks are generated by the propagation and coalescence of these microcracks.Figure 2
The parallel bond model (PBM) in PFC. (a) Schematic diagram illustrating the micromechanical behavior of the particle contacts; (b) deformation and failure mechanism of the bond under tensile and shear stress after [33]; (c) micromechanical parameters of the bond at contacts.
(a)
(b)
(c)


### 2.2. Smooth Joint Model (SJM)

Traditionally, joints are modeled in PFC by either removing the particles that fall within a specific range or assigning a low coefficient of friction to particles at certain contacts. However, this kind of treatment can result in an inherent microscale roughness of the joint surface and leads to an unrealistic increase in shear strength and dilation during shearing along the shear plane. To overcome this problem, a smooth joint model (SJM) was proposed [33]. A typical example of a two-dimensional SJM is shown in Figure 3. In this model, after the joint plane is defined, a smooth joint (SJ) is assigned to the contacts between the particles. The centers of these particles lie on the opposite side of the defined joint plane. At these contacts, the parallel bonds are removed and the smooth joints are defined in a direction parallel with the joint plane.Figure 3
A typical SJ contact model in two dimensions after [33].Forces and relative displacement can be expressed in local coordinates of the joint plane:(1)F=Fntj+Fs,(2)U=Unn^j+Us.The strength of a SJ contact follows the Coulomb sliding model. The maximum possible magnitude of the shear forceFs∗ depends upon the coefficient of friction μj:
(3)Fs∗=Fn⋅μj,The value of the shear force is updated through recomputation at the end of each time step. IfFs<Fs∗, then Fs=Fs′. The forces are updated as follows:
(4)Fn≔Fn+k¯njAΔUn,(5)Fs′=Fs−k¯sjAΔUs.IfFs′>Fs∗, then sliding occurs at the contacts between the particles. The forces are updated through
(6)Fs=Fs∗,(7)Fn≔Fn+ΔUs∗tanφj⋅k¯njA=Fn+Fs′−Fs∗k¯sj⋅k¯njtanφj.For the SJ, the gap value between the particles can be defined. However, the SJ can only be activated when the gap value diminishes to zero. Therefore, by defining a specific gap value, joints with non-zero aperture at the initial can be modeled.


### 2.3. Coupled Hydromechanical Model

In PFC, the algorithm governing fluid flow is developed out of a network flow model. The fluid flow model is based on two important assumptions: (1) the flow domain is characterized by the polygons formed by straight lines linking the center point of the adjacent particles in a random plane, and (2) the flow domains are connected by the flow channels between the adjacent particles.Figure4 shows a schematic diagram of this fluid flow model. Due to the pressure difference in each flow domain, the fluid flows from one flow domain to another. This process is assumed as one in which the fluid passes in between two parallel plates in laminar flow. The flow process is governed by the cubic law. Therefore, the flow volume Q is expressed as [34]
(8)Q=a312μ⋅PA−PBL.Figure 4
DEM used in coupled hydromechanical simulation; (a) fluid domains (blue polygons), centers of domains (red circles), and flow path (red lines) constituting the hydraulic network; (b) distribution of hydraulic pressure; (c) domain and flow channel; (d) hydraulic pressure build-up; (e) change of aperture.In the fluid flow process, the pressure incrementΔP in the reservoir can be calculated with [35]
(9)ΔP=KfVd∑Q⋅Δt−ΔVd.Figure4(b) illustrates how pore pressure is formed within a rock matrix, where the size of the yellow spot is in direct proportion to the magnitude of the pore pressure. The changes to the stress state resulting from bond breakage at particle contacts are determined by the bond strength and Mohr-Coulomb failure criterion. For an open fracture, the permeability is assumed to be infinite. The fluid flow algorithm is based on two assumptions that fluid domains are connected by fluid flow channels (Figures 4(c) and 4(d)).The fluid pressure applied to the contact between the particles results in the deformation of the reservoir. Therefore, the hydraulic aperture in the proposed model in this study is not the actual aperture of the flow path. It is an empirical value that characterizes the permeability of the rock under different confining pressures. The aperture of the flow path changes after the fluid is injected into the model (Figure4(e)). In the fluid flow process, the aperture of the flow path is correlated with the normal stress at the contact. The relation is expressed as
(10)a=ainf+azeo−ainfeξσn.Equation (10) characterizes the relation between a and ainf in hydraulic fracturing. As σn tends to infinity, a gets close to ainf. In (8), Q applies specifically to the fluid that flows from one flow domain to another. It is a microscopic quantity. Therefore, it cannot be used to calculate the permeability coefficient of the rock, which is a macroscopic quantity. It can be obtained as follows:
(11)k=112V∑pipesLa3.As the fluid pressure exceeds the tensile or shear strength of the bonded particles, cracks start to initiate between the adjacent particles. By this time, the aperture of the flow path is assumed to be infinite. This is an irreversible process, meaning the bonding between the particles loses its efficacy or is no longer effective. However, it should also be noted thatainf is not left to be infinite. Instead, it is still set to a value to represent a different flow state.Once the cracks are initiated, the infinite setting might undermine the stability of the simulation, delivering an adverse impact on the modeling results. This is because the sudden formation of the cracks results in the instantaneous flow of the fluids between the fluid passages and causes momentary changes in fluid pressure. In this paper, after the breakage of the particles, the fluid pressurePf′ is assumed to be the mean value of the pressures in the two flow domains prior to particle debonding [36, 37]:
(12)Pf′=PfA+PfB2.


## 2.1. Parallel Bonded Model

The PFC, a simulation code based on the DEM is an effective numerical tool for simulating the behavior of jointed rock mass. In PFC, the rock matrix is modeled as an assembly of circular particles. Though PFC is based on DEM, it can still be employed to model the deformation characteristics and mechanical behavior of a continuous medium. The parallel bonded model (PBM) is used to define and simulate the particle contacts as it can deliver more realistic results in modeling rock fracture.Figure2(a) schematizes the micromechanical behavior of the particle contacts. At particle contacts, two sets of springs provide the normal and shear stiffness [31]. As shown in Figures 2(b) and 2(c), the first set consists of two springs that stand for normal and shear stiffness, which are denoted as kn and ks, respectively. In the model, the particles are idealized as a rigid body that is unbreakable and undeformable. However, the particles can overlap if the compressive force is great enough. The second set of springs provides the parallel bond normal and shear stiffness denoted as k¯n and k¯s, which are uniformly distributed over a disc-shaped cross section lying on the contact plane. In addition to the normal and shear stiffness, the parallel bonds also have normal and shear strength. A parallel bond is active over a finite area and can resist the rotation of particles. When the normal or shear stress exceeds the bond strength, the parallel bond is broken and bond stiffness is removed. Once the parallel bond is broken under shear stress, the shear strength falls to its residual value whereas upon breakage of a bond under tensile stress, the tensile strength is set to zero. The residual shear strength depends on the normal force and the friction coefficient of particles at contacts [32]. Microcracks are produced in PBM once the bonds fail due to external loading, and macrocracks are generated by the propagation and coalescence of these microcracks.Figure 2
The parallel bond model (PBM) in PFC. (a) Schematic diagram illustrating the micromechanical behavior of the particle contacts; (b) deformation and failure mechanism of the bond under tensile and shear stress after [33]; (c) micromechanical parameters of the bond at contacts.
(a)
(b)
(c)


## 2.2. Smooth Joint Model (SJM)

Traditionally, joints are modeled in PFC by either removing the particles that fall within a specific range or assigning a low coefficient of friction to particles at certain contacts. However, this kind of treatment can result in an inherent microscale roughness of the joint surface and leads to an unrealistic increase in shear strength and dilation during shearing along the shear plane. To overcome this problem, a smooth joint model (SJM) was proposed [33]. A typical example of a two-dimensional SJM is shown in Figure 3. In this model, after the joint plane is defined, a smooth joint (SJ) is assigned to the contacts between the particles. The centers of these particles lie on the opposite side of the defined joint plane. At these contacts, the parallel bonds are removed and the smooth joints are defined in a direction parallel with the joint plane.Figure 3
A typical SJ contact model in two dimensions after [33].Forces and relative displacement can be expressed in local coordinates of the joint plane:(1)F=Fntj+Fs,(2)U=Unn^j+Us.The strength of a SJ contact follows the Coulomb sliding model. The maximum possible magnitude of the shear forceFs∗ depends upon the coefficient of friction μj:
(3)Fs∗=Fn⋅μj,The value of the shear force is updated through recomputation at the end of each time step. IfFs<Fs∗, then Fs=Fs′. The forces are updated as follows:
(4)Fn≔Fn+k¯njAΔUn,(5)Fs′=Fs−k¯sjAΔUs.IfFs′>Fs∗, then sliding occurs at the contacts between the particles. The forces are updated through
(6)Fs=Fs∗,(7)Fn≔Fn+ΔUs∗tanφj⋅k¯njA=Fn+Fs′−Fs∗k¯sj⋅k¯njtanφj.For the SJ, the gap value between the particles can be defined. However, the SJ can only be activated when the gap value diminishes to zero. Therefore, by defining a specific gap value, joints with non-zero aperture at the initial can be modeled.


## 2.3. Coupled Hydromechanical Model

In PFC, the algorithm governing fluid flow is developed out of a network flow model. The fluid flow model is based on two important assumptions: (1) the flow domain is characterized by the polygons formed by straight lines linking the center point of the adjacent particles in a random plane, and (2) the flow domains are connected by the flow channels between the adjacent particles.Figure4 shows a schematic diagram of this fluid flow model. Due to the pressure difference in each flow domain, the fluid flows from one flow domain to another. This process is assumed as one in which the fluid passes in between two parallel plates in laminar flow. The flow process is governed by the cubic law. Therefore, the flow volume Q is expressed as [34]
(8)Q=a312μ⋅PA−PBL.Figure 4
DEM used in coupled hydromechanical simulation; (a) fluid domains (blue polygons), centers of domains (red circles), and flow path (red lines) constituting the hydraulic network; (b) distribution of hydraulic pressure; (c) domain and flow channel; (d) hydraulic pressure build-up; (e) change of aperture.In the fluid flow process, the pressure incrementΔP in the reservoir can be calculated with [35]
(9)ΔP=KfVd∑Q⋅Δt−ΔVd.Figure4(b) illustrates how pore pressure is formed within a rock matrix, where the size of the yellow spot is in direct proportion to the magnitude of the pore pressure. The changes to the stress state resulting from bond breakage at particle contacts are determined by the bond strength and Mohr-Coulomb failure criterion. For an open fracture, the permeability is assumed to be infinite. The fluid flow algorithm is based on two assumptions that fluid domains are connected by fluid flow channels (Figures 4(c) and 4(d)).The fluid pressure applied to the contact between the particles results in the deformation of the reservoir. Therefore, the hydraulic aperture in the proposed model in this study is not the actual aperture of the flow path. It is an empirical value that characterizes the permeability of the rock under different confining pressures. The aperture of the flow path changes after the fluid is injected into the model (Figure4(e)). In the fluid flow process, the aperture of the flow path is correlated with the normal stress at the contact. The relation is expressed as
(10)a=ainf+azeo−ainfeξσn.Equation (10) characterizes the relation between a and ainf in hydraulic fracturing. As σn tends to infinity, a gets close to ainf. In (8), Q applies specifically to the fluid that flows from one flow domain to another. It is a microscopic quantity. Therefore, it cannot be used to calculate the permeability coefficient of the rock, which is a macroscopic quantity. It can be obtained as follows:
(11)k=112V∑pipesLa3.As the fluid pressure exceeds the tensile or shear strength of the bonded particles, cracks start to initiate between the adjacent particles. By this time, the aperture of the flow path is assumed to be infinite. This is an irreversible process, meaning the bonding between the particles loses its efficacy or is no longer effective. However, it should also be noted thatainf is not left to be infinite. Instead, it is still set to a value to represent a different flow state.Once the cracks are initiated, the infinite setting might undermine the stability of the simulation, delivering an adverse impact on the modeling results. This is because the sudden formation of the cracks results in the instantaneous flow of the fluids between the fluid passages and causes momentary changes in fluid pressure. In this paper, after the breakage of the particles, the fluid pressurePf′ is assumed to be the mean value of the pressures in the two flow domains prior to particle debonding [36, 37]:
(12)Pf′=PfA+PfB2.


## 3. Model Steps and Microparameter Calibration

### 3.1. Model Steps

Using DEM, Wang et al. [38] performed a hydraulic fracturing study on five typical coal mass models, i.e., intact coal, layered jointed coal, vertical jointed coal, orthogonal jointed coal, and synthetic jointed coal, and compared the HFN pattern produced in different models. A conceptual model for a jointed rock mass with multiple planes of nonpersistent joints is shown in Figure 5. The dimension of the model is 1.0 × 1.0 m, and the particle size falls within a range of 60 to 88 mm.Figure 5
Joint geometrical model for hydraulic fracturing simulation.To create a coupled hydromechanical model in PFC, a confined square frame with four sets of walls is produced in the first place. Then, the shale reservoir model is established by the following steps: (1) Generate particles and arrange them randomly to fill the square frame. (2) Stop the moving particles until they are not embedded into each other. (3) Eliminate floating particles to ensure that there is no great dispersion in the magnitude of the hydraulic aperture. (4) Assemble SJM in specified zones and assemble PBM in the remaining zones. (5) Apply 10 MPa and 5 MPa to the model in the horizontal and vertical direction, respectively. (6) Excavate a borehole with a diameter of 50 mm in the center of the model and run the model until reaching equilibrium. (7) Form the flow domains in the model using an independent subprogram.In this study, the SJM involves microparameters such as joint orientation angleβ, joint persistency k, joint spacing d, joint step angle γ, and joint aperture e. Figure 5 clearly illustrates the above microparameters. The joint orientation angle is defined as the angle with respect to horizontal principal stress in plane SH, which changes from 0° to 90° at an interval of 15°. Joint persistency k is defined by the joint length Lj and the rock bridge length Lr as
(13)k=LjLj+Lr.Joint spacingd is defined as the straight-line distance between the centerline of the adjacent parallel joints. The joint step angle γ is defined as the angle of the line connecting the end of the present joint set with the starting point of the next joint set. Joint aperture e is defined as the width of the SJ.


### 3.2. Microparameter Calibration

The microparameters of the jointed rock formations modeled in PFC cannot be directly obtained from the physical experiment. However, these microparameters can be calibrated by the physical tests performed on specimens with various joint orientation angles. Firstly, Brazilian disc tests were carried out on the samples. Secondly, sensitivity analyses were conducted on the mechanical parameters of the joint, including joint strength, ratio of cohesion to tensile strength, normal stiffness, shear stiffness, and angle range. The authors’ previous studies show that [39] (1) the strength of the SJ in the proposed model that represents the joints exerts a main effect on the strength of the reservoir; (2) the elastic modulus of the reservoir is related to the contact stiffness of the joints; and (3) the ratio of the SJ cohesion to its tensile strength mainly affects the number of the cracks generated in the reservoir, thus having further effects on the failure mode of the rock mass. Through trial and error, the authors finally obtain the microparameters that replicate the macromechanical behavior of the reservoir, which are tabulated in Tables 1 and 2.Table 1
Microparameters used in shale matrix and joint plane after calibration.
Parameters
Nonbrittle mineral
Brittle mineral
Bond
Particle properties
Density (kg·m−3)
1750
2500
N/A
Particle radiusRmin (mm)
0.25
0.25
N/A
Young’s modulus/GPa
36.2
68.2
N/A
Normal-to-shear stiffness ratiokn/ks
0.8
0.8
N/A
Friction coefficientμ
0.67
0.67
N/A
Contact properties
Normal stiffnessk¯n (GPa/m)
72,400
136,400
78,000
Shear stiffnessk¯s (GPa/m)
90,500
130,500
97,500
Tensile strengthσ¯ (MPa)
25.6
32.8 ± 3
16.5 ± 5
Cohesionc¯ (MPa)
30.3
37.6 ± 3
18.6 ± 5Table 2
Microparameters used in SJ contact for simulated specimen after calibration.
Parameters
Values
Smooth joint model (SJM)
Normal stiffnesssj_kn (GPa/m)
54,000
Shear stiffnesssj_ks (GPa/m)
68,000
Tensile strengthσ¯s (MPa)
12.5 ± 2
Cohesionc¯s (MPa)
14.5 ± 1
Friction angleφs (°)
0To validate the accuracy of the microparameters, we carried out laboratory tests. The angle between the joint plane and the vertical direction is assumed asθ. Figure 6 summarizes the failure mode of the specimens tested at an angle ranging from 0° to 90° and presents the photos showing the distribution of the microcracks and macrocracks. Figure 7 presents the comparison results regarding failure load and splitting modulus in the physical experiment and numerical simulation. The yellow stands for the joint and the gray stands for shale matrix in figures of the second column. The comparison shows that the simulation results show good agreement with the experimental results. It should be noted that while performing the calibration on the microparameters of the shale reservoirs, we split the shale matrix into two types, i.e., brittle mineral and nonbrittle mineral, to have a model that delivers more precise results. For more particulars, readers can refer to the authors’ previous publication [39].Figure 6
Comparison of failure modes between experimental data and DEM simulation data. Legend of microscopic failure mode.Figure 7
Comparison of (a) failure load in Brazilian tests and (b) splitting modulus from experimental data and DEM simulation data.
(a)
(b)In Figure6, the thickness of the core is 25 mm, which is rather small. In physical experiments, the specimens used in the Brazilian disc test are transversely isotropic shale, meaning that the specimens do not exhibit anisotropic characteristics along the z-direction. The specimens in the Brazilian disc test fail under quasistatic conditions. In addition, the shale is very brittle. Therefore, from fracture initiation to failure, this process takes a very short time. The process for the fractures to evolve gradually is virtually nonexistent. In a nutshell, it is safe to say that the shale is isotropic along the z-direction.It should be noted that the 3D version of PFC is also able to model the above physical experiment. However, using PFC3D to model this experiment is very time-consuming and uneconomical. In addition, it is difficult to control the loading and boundary conditions if using PFC3D to perform the modeling. Instead, PFC2D is easy to use and can better illustrate the evolution of the fractures. Therefore, PFC2D is employed in this study.A great number of scholars and researchers have validated the feasibility of modeling the flow dynamics in solid materials by using DEM [40–42]. However, to model hydraulic fracturing using DEM, not only should the microparameters of the matrix and joint in the shale reservoir be determined but also the microparameters governing fluid flow should be taken into account. Therefore, we performed calibration on the hydraulic fracturing breakdown pressure in a homogeneous shale reservoir. To date, no method has been available to calibrate the fluid parameters in a coupled hydromechanical model. However, Haimson and Fairhurst [43] and Fairhurst [44] put forward a regression equation for calculating the breakdown pressure
(14)σbp=3Sh−SH−P0+σt,where σbp is the breakdown pressure of the reservoir, P0 is the initial pore pressure, and σt is the tensile strength of the reservoir.To validate the reliability of using the bonded particle method to model hydraulic fracturing, we set the initial pore pressure as 0 in this study. The tensile strength of the homogeneous shale reservoir is selected as 10.76 MPa. For the confining pressure in theX-direction (SH), it is 20 and 10 MPa; the confining pressure in the Y-direction (Sh). It varies from 5 to 10 MPa. The microparameters governing fluid flow in the model are listed in Table 3. The comparison between the numerical simulation and the values derived from the theoretical equation under different confining pressure ratios are shown in Figure 8. Results demonstrate that simulation results fit well with the theoretical values regardless of a certain disagreement. The disagreement can be reasonably explained. Equation (14) is derived based on elastic mechanics that assumes the rock mass as elastic, isotropic, and homogeneous. However, this assumption is against the reality.Table 3
Computational parameters of hydraulic properties after calibration.
Fluid parameters
Unit
Values
Injection rate
m3·s−1
3.0 × 10−6
Fluid bulk modulus (Kf)
GPa
2.35
Fluid dynamic viscosity (μ)
Pa·s
1.40 × 10−4
Initial fluid aperture (azeo)
m
2.85 × 10−6
Infinite fluid aperture (ainf)
m
2.85 × 10−7Figure 8
Comparison of breakdown pressure under different confining pressure ratios obtained by the theoretical equation and numerical simulation.


## 3.1. Model Steps

Using DEM, Wang et al. [38] performed a hydraulic fracturing study on five typical coal mass models, i.e., intact coal, layered jointed coal, vertical jointed coal, orthogonal jointed coal, and synthetic jointed coal, and compared the HFN pattern produced in different models. A conceptual model for a jointed rock mass with multiple planes of nonpersistent joints is shown in Figure 5. The dimension of the model is 1.0 × 1.0 m, and the particle size falls within a range of 60 to 88 mm.Figure 5
Joint geometrical model for hydraulic fracturing simulation.To create a coupled hydromechanical model in PFC, a confined square frame with four sets of walls is produced in the first place. Then, the shale reservoir model is established by the following steps: (1) Generate particles and arrange them randomly to fill the square frame. (2) Stop the moving particles until they are not embedded into each other. (3) Eliminate floating particles to ensure that there is no great dispersion in the magnitude of the hydraulic aperture. (4) Assemble SJM in specified zones and assemble PBM in the remaining zones. (5) Apply 10 MPa and 5 MPa to the model in the horizontal and vertical direction, respectively. (6) Excavate a borehole with a diameter of 50 mm in the center of the model and run the model until reaching equilibrium. (7) Form the flow domains in the model using an independent subprogram.In this study, the SJM involves microparameters such as joint orientation angleβ, joint persistency k, joint spacing d, joint step angle γ, and joint aperture e. Figure 5 clearly illustrates the above microparameters. The joint orientation angle is defined as the angle with respect to horizontal principal stress in plane SH, which changes from 0° to 90° at an interval of 15°. Joint persistency k is defined by the joint length Lj and the rock bridge length Lr as
(13)k=LjLj+Lr.Joint spacingd is defined as the straight-line distance between the centerline of the adjacent parallel joints. The joint step angle γ is defined as the angle of the line connecting the end of the present joint set with the starting point of the next joint set. Joint aperture e is defined as the width of the SJ.


## 3.2. Microparameter Calibration

The microparameters of the jointed rock formations modeled in PFC cannot be directly obtained from the physical experiment. However, these microparameters can be calibrated by the physical tests performed on specimens with various joint orientation angles. Firstly, Brazilian disc tests were carried out on the samples. Secondly, sensitivity analyses were conducted on the mechanical parameters of the joint, including joint strength, ratio of cohesion to tensile strength, normal stiffness, shear stiffness, and angle range. The authors’ previous studies show that [39] (1) the strength of the SJ in the proposed model that represents the joints exerts a main effect on the strength of the reservoir; (2) the elastic modulus of the reservoir is related to the contact stiffness of the joints; and (3) the ratio of the SJ cohesion to its tensile strength mainly affects the number of the cracks generated in the reservoir, thus having further effects on the failure mode of the rock mass. Through trial and error, the authors finally obtain the microparameters that replicate the macromechanical behavior of the reservoir, which are tabulated in Tables 1 and 2.Table 1
Microparameters used in shale matrix and joint plane after calibration.
Parameters
Nonbrittle mineral
Brittle mineral
Bond
Particle properties
Density (kg·m−3)
1750
2500
N/A
Particle radiusRmin (mm)
0.25
0.25
N/A
Young’s modulus/GPa
36.2
68.2
N/A
Normal-to-shear stiffness ratiokn/ks
0.8
0.8
N/A
Friction coefficientμ
0.67
0.67
N/A
Contact properties
Normal stiffnessk¯n (GPa/m)
72,400
136,400
78,000
Shear stiffnessk¯s (GPa/m)
90,500
130,500
97,500
Tensile strengthσ¯ (MPa)
25.6
32.8 ± 3
16.5 ± 5
Cohesionc¯ (MPa)
30.3
37.6 ± 3
18.6 ± 5Table 2
Microparameters used in SJ contact for simulated specimen after calibration.
Parameters
Values
Smooth joint model (SJM)
Normal stiffnesssj_kn (GPa/m)
54,000
Shear stiffnesssj_ks (GPa/m)
68,000
Tensile strengthσ¯s (MPa)
12.5 ± 2
Cohesionc¯s (MPa)
14.5 ± 1
Friction angleφs (°)
0To validate the accuracy of the microparameters, we carried out laboratory tests. The angle between the joint plane and the vertical direction is assumed asθ. Figure 6 summarizes the failure mode of the specimens tested at an angle ranging from 0° to 90° and presents the photos showing the distribution of the microcracks and macrocracks. Figure 7 presents the comparison results regarding failure load and splitting modulus in the physical experiment and numerical simulation. The yellow stands for the joint and the gray stands for shale matrix in figures of the second column. The comparison shows that the simulation results show good agreement with the experimental results. It should be noted that while performing the calibration on the microparameters of the shale reservoirs, we split the shale matrix into two types, i.e., brittle mineral and nonbrittle mineral, to have a model that delivers more precise results. For more particulars, readers can refer to the authors’ previous publication [39].Figure 6
Comparison of failure modes between experimental data and DEM simulation data. Legend of microscopic failure mode.Figure 7
Comparison of (a) failure load in Brazilian tests and (b) splitting modulus from experimental data and DEM simulation data.
(a)
(b)In Figure6, the thickness of the core is 25 mm, which is rather small. In physical experiments, the specimens used in the Brazilian disc test are transversely isotropic shale, meaning that the specimens do not exhibit anisotropic characteristics along the z-direction. The specimens in the Brazilian disc test fail under quasistatic conditions. In addition, the shale is very brittle. Therefore, from fracture initiation to failure, this process takes a very short time. The process for the fractures to evolve gradually is virtually nonexistent. In a nutshell, it is safe to say that the shale is isotropic along the z-direction.It should be noted that the 3D version of PFC is also able to model the above physical experiment. However, using PFC3D to model this experiment is very time-consuming and uneconomical. In addition, it is difficult to control the loading and boundary conditions if using PFC3D to perform the modeling. Instead, PFC2D is easy to use and can better illustrate the evolution of the fractures. Therefore, PFC2D is employed in this study.A great number of scholars and researchers have validated the feasibility of modeling the flow dynamics in solid materials by using DEM [40–42]. However, to model hydraulic fracturing using DEM, not only should the microparameters of the matrix and joint in the shale reservoir be determined but also the microparameters governing fluid flow should be taken into account. Therefore, we performed calibration on the hydraulic fracturing breakdown pressure in a homogeneous shale reservoir. To date, no method has been available to calibrate the fluid parameters in a coupled hydromechanical model. However, Haimson and Fairhurst [43] and Fairhurst [44] put forward a regression equation for calculating the breakdown pressure
(14)σbp=3Sh−SH−P0+σt,where σbp is the breakdown pressure of the reservoir, P0 is the initial pore pressure, and σt is the tensile strength of the reservoir.To validate the reliability of using the bonded particle method to model hydraulic fracturing, we set the initial pore pressure as 0 in this study. The tensile strength of the homogeneous shale reservoir is selected as 10.76 MPa. For the confining pressure in theX-direction (SH), it is 20 and 10 MPa; the confining pressure in the Y-direction (Sh). It varies from 5 to 10 MPa. The microparameters governing fluid flow in the model are listed in Table 3. The comparison between the numerical simulation and the values derived from the theoretical equation under different confining pressure ratios are shown in Figure 8. Results demonstrate that simulation results fit well with the theoretical values regardless of a certain disagreement. The disagreement can be reasonably explained. Equation (14) is derived based on elastic mechanics that assumes the rock mass as elastic, isotropic, and homogeneous. However, this assumption is against the reality.Table 3
Computational parameters of hydraulic properties after calibration.
Fluid parameters
Unit
Values
Injection rate
m3·s−1
3.0 × 10−6
Fluid bulk modulus (Kf)
GPa
2.35
Fluid dynamic viscosity (μ)
Pa·s
1.40 × 10−4
Initial fluid aperture (azeo)
m
2.85 × 10−6
Infinite fluid aperture (ainf)
m
2.85 × 10−7Figure 8
Comparison of breakdown pressure under different confining pressure ratios obtained by the theoretical equation and numerical simulation.


## 4. Results and Analysis

To understand the geometrical parameters of the joints on the breakdown pressure, the number of hydraulic fractures, and the pattern of the generated HFN, we performed a series of sensitivity analyses using the calibrated model. In each study group, only the factor of interest is different and the other factors such as joint plane configuration and hydraulic parameters are kept the same. The fluid injection time is kept constant at 10 s. In this study, the factors of interest include joint orientation angleβ, joint persistency k, joint spacing d, joint step angle γ, and joint aperture e.The coupled hydromechanical model in this study can calculate the contact force between the particles, which is used to extrapolate the equivalent continuum stress. As the amount of the injection fluid increases, the bond between the particles breaks, resulting in HF. The cracks in the formation can be used for quantitative or qualitative analysis. To evaluate the macrocrack through monitoring the accumulative quantity of the microcracks has been proven effective by many scholars and researchers [45, 46]. Using an independent subprogram, various types of hydraulic cracks can be identified and recorded, i.e., tensile hydraulic fracture caused by shale matrix (THFSM), shear hydraulic fracture caused by shale matrix (SHFSM), tensile hydraulic fracture caused by joint plane (THFJP), and shear hydraulic fracture caused by joint plane (SHFJP).In a homogeneous shale reservoir, HF in general propagates along the direction of the maximum principal stress. However, this is not applicable in jointed shale reservoirs. The modeling results demonstrate that the geometrical parameters of the joints have a significant effect on the HFN pattern formed in jointed shale reservoirs.In general, after hydraulic fracturing, the HFN pattern in jointed shale reservoirs fall into four distinct types, i.e., crossing mode, tip-to-tip mode, step path mode and opening mode. Figure9 illustrates the four types of HFN patterns characterized by crack distribution and stress concentration, where more dark arrows indicate greater stress concentration. Figure 10 presents the evolution of the four types of cracks with the injection time.
(1)
Crossing mode. This HFN pattern results from the stress concentration near the shale matrix between the joints and the slight stress between the rock bridges (Figure 9(a)). In the meantime, a small amount of shear HF initiates from the internal zones of the joints and propagates to the boundary along the joint orientation. In this scenario, HF no longer propagates in a direction normal to the principal stress. Instead, it propagates in a direction nearly vertical to the joints. In this HFN pattern, THFSM gains the dominance over other types of cracks, accounting for about 60% in number. As some joints are activated during HF propagation, the HF generated in the joints accounts for 30% in number (Figure 10(a)).
(2)
Tip-to-tip mode. This HFN pattern results from the stress concentration near the joint tips and the presence of a low stress zone between the rock bridges, as shown in Figure 9(b). The tensile HF initiates near the joint tips and coalesces with the HF initiated near other joint tips. The shear HF generates and propagates to the boundary along the joint. In this pattern, the number of THFSM further outweighs that of the other three HFs, which is about two times of the number of the cracks combined (Figure 10(b)).
(3)
Step path mode. At the beginning of fluid injection, a slight amount of HF initiates near the joint tips or within the internal zone. However, as more fluids are injected into the borehole, HF does not traverse the joints or propagate along the joint tip. Instead, along the joints, it initiates from the other end of the joints, resulting in joint slip (Figure 9(c)). In this pattern, the stress concentration zone forms in a direction parallel with the principal stress. Due to the propagation of HF along the joint, the number of SHFJP increases. Nonetheless, throughout the whole injection process, THFSM still maintains its role as the most dominant crack type (Figure 10(c)).
(4)
Opening mode. As joint aperture increases, the hydraulic pressure induced by fluid injection is completely absorbed by the joints. As a result, more fluid flows into the joints, generating THFJP. Stress concentration occurs only near the borehole, thus incapacitating the generation of HF in the shale matrix (Figure 9(d)). Due to the inherent low strength of the distant joints, the perturbation caused by the injection pressure generates HF. Therefore, SHFJP grows significantly in number, usurping the dominance of the THFSM (Figure 10(d)).Figure 9
Contact force distribution in coupled hydromechanical model: (a) crossing mode, (b) tip-to-tip mode, (c) step path mode, and (d) opening mode.
(a)
(b)
(c)
(d)Figure 10
Hydraulic fracture (HF) evolution for (a) crossing mode, (b) tip-to-tip mode, (c) step path mode, and (d) opening mode.
(a)
(b)
(c)
(d)


### 4.1. Effect of Joint Orientation Angle on Hydraulic Fracture Network Propagation

To investigate the effect ofβ, we keep other microparameters constant, i.e., k=0.7, d=120mm, γ=117°, and e=1.2×10−4m. Then, we performed numerical simulation with β changing from 0° (horizontal joint) to 90° (vertical joint) at an interval of 15°. The modeling results are shown in Figures 11–13, which separately present the effect of β on the breakdown pressure, the number of HF, and the HFN pattern. The four colors in Figure 13 represent the tensile and shear cracks generated in the shale matrix and the joint. To study to what extent the breakdown pressure in jointed shale reservoir decreases as opposed to that in homogeneous shale reservoir, the breakdown pressure σbp in jointed shale reservoir is normalized to the breakdown pressure σHo in a homogeneous shale reservoir, which is 13.22 MPa.Figure 11
Effect of joint orientation (β) on breakdown pressure.During hydraulic fracturing,σbp reaches its maximum value when β is 90° (vertical joint). While β is 0°, the joint is in a direction parallel with the maximum principal stress. In this case, stress concentrates near the joint tip, and the tensile HF initiates near the joint tips and coalesces with the HF initiated near other joint tips. In the meantime, SHFJP is generated in the joints. Therefore, this HFN pattern is characterized by tip-to-tip mode. As can be seen from Figure 12(a), in this case, the number of HF reaches a peak, amounting to about 2000. When β=15° and 30°, the joint is angled with respect to the direction of the maximum principal stress. In this case, stress concentrates near the joint tip and in the contacting areas between the shale matrix. Moreover, the breakdown pressure declines as shown in Figure 11. The generated HF forms a pattern that is combined by tip-to-tip mode and crossing mode. As β increases to 45° or 60°, stress concentrates within the shale matrix. Under the combined action of the joint orientation and the maximum principal stress, tensile HF propagates in a direction nearly vertical to the joint, forming a HFN pattern characterized by the crossing mode. In this pattern, HF firstly initiates in the internal zones of the joints, followed by the propagation and coalesce of HF between the joints. Finally, HF penetrates the joints. In this case, the breakdown pressure σbp is roughly about 30% of the breakdown pressure σHo. According to the Coulomb criterion, the jointed shale reservoir has its minimum tensile strength when β varies from 45° to 60°. In addition, the previous studies have validated that the breakdown pressure is positively correlated with the tensile strength of the shale reservoir [47]. Therefore, theoretical solution and numerical simulation deliver a consistent result regarding the breakdown pressure, which further validates the correctness of the microparameters used in this study.Figure 12
Effect of joint orientation (β) on hydraulic fracture (HF): (a) number and (b) percentage.
(a)
(b)Figure 13
Effect of joint orientation (β) on the propagation mode of the hydraulic fracture network.Figure 14
Effect of joint persistency (k) on breakdown pressure.Figure 15
Effect of joint persistency (k) on the number and percentage of different types of hydraulic fractures.Figure 16
Effect of joint persistency (k) on the propagation mode of the hydraulic fracture network.Figure 17
Effect of joint spacing (d) on breakdown pressure.Figure 18
Effect of joint spacing (d) on the number and percentage of different types of hydraulic fracture (HF).Figure 19
Effect of joint spacing (d) on the propagation mode of the hydraulic fracture network (HFN).Figure 20
Effect of the joint step angle (γ) on breakdown pressure.Figure 21
Effect of the joint step angle (γ) on the number and percentage of different types of hydraulic fracture (HF).Figure 22
Effect of the joint step angle (γ) on the propagation mode of the hydraulic fracture network (HFN).Figure 23
Effect of the joint aperture (e) on breakdown pressure.Figure 24
Effect of the joint aperture (e) on the number and percentage of different types of hydraulic fracture (HF).Figure 25
Effect of the joint aperture (e) on the propagation mode of the hydraulic fracture network (HFN).Afterβ reaches 75°, stress is less concentrated near the joint tip. Therefore, the breakdown pressure increases and is greater than the breakdown pressure while β is 0°. When HF interacts with the joints, the joints are activated. As a result, HF initiates at the other end of the joint and propagates. However, when β=75° or 90°, cracks do not cover a wide area in the shale reservoir except for areas surrounding the injection borehole. Therefore, the number of HF is smaller than that when β is 0° or 15°. However, it should be noted that the proportion of different types of HF remains unchanged despite the changes in the number of HF with β (Figure 12(b)).The modeling results demonstrate thatβ exerts a significant influence on the breakdown pressure, the number of HF, and the HFN pattern. Therefore, β will also be considered while investigating the effect of other microparameters of the joints.


### 4.2. Effect of Joint Persistency on Hydraulic Fracture Network Propagation

To study the effect ofk, we keep other microparameters constant except β (d=120mm, γ=117°, e=1.2×10−4m, and β changes from 0° to 90° at an interval of 15°). k is changed by varying the length of the joint while the length of the rock bridge is kept unchanged (Lr=60mm). Figures 14–16 present the effect of k on the breakdown pressure, the number of HF, and the HFN pattern. Results show that as k increases, β has a more pronounced effect on the breakdown pressure and the HFN pattern, as shown in Figures 14 and 16. k exerts a more apparent effect on the number of HF when β varies from 45° to 60° as shown in Figure 15.As shown in Figure16, the HF pattern in all the shale reservoirs is characterized by tip-to-tip mode when β is 0°. The number of HF varies slightly from one shale reservoir to another. In this scenario, an increase in k leads to a decrease in breakdown pressure (Figure 14). This is due mainly to the aggravation of stress concentration near the joint tips. In fact, stress has already concentrated around the joint tips even when k is 0.6. Therefore, the increase in k further aggravates stress concentration near the joint tips. This is why the HFN pattern remains unchanged while the breakdown pressure diminishes. As β increases to 15° and 30°, with a small k (k=0.6 or 0.7), the HFN pattern formed is a combination of the tip-to-tip mode and the crossing mode. However, the HFN pattern is still characterized by the tip-to-tip mode when k changes to 0.8. As a larger k results in a greater joint density, stress is more easily concentrated around the joint tips. Therefore, it is easier for HF to initiate at the joint tips and to propagate further (Figure 16). As a result, the breakdown pressure significantly drops (Figure 14). According to the theory of fracture mechanics, greater joint length results in greater stress intensity factor. Therefore, greater stress concentration occurs around the joint tips. This is why as k increases the breakdown pressure drops.Asβ increases to 45° or 60°, the HFN pattern is characterized by the crossing mode. an increase in k leads to an increase in the number of HF. In addition, the HF generated in joints accounts for a greater proportion, with SHFJP in particular. When k reaches the extreme values (k=0.8 and k=0.6), the difference in the number of HF amounts to over 400 and HF propagates in a direction nearly vertical to the joint, as shown in Figure 15. As a result, stress concentration gradually spreads to rock bridge areas, causing a difference in the breakdown pressure and the number of HF. Stress is less concentrated around the joint tips when β is equal to or greater than 75°. Stress is more evenly distributed across the rock bridge area. Therefore, k does not exert an apparent effect on the breakdown pressure. Its effect on the number of HF also weakens. The HFN pattern converts from the crossing mode to the step path mode.Therefore, we can conclude that (1) when keeping the joint orientation angle and other parameters unchanged, an increase ink leads to a decrease in the breakdown pressure and the number of HF; (2) an increase in k causes no perceptible changes to the proportion of HF; and (3) though k has an effect on the HFN pattern, the effect is not as significant as that of β.


### 4.3. Effect of Joint Spacing on Hydraulic Fracture Network Propagation

To study the effect ofd, we keep other microparameters constant except β (k=0.7, γ=117°, e=1.2×10−4m, and β changes from 0° to 90° at an interval of 15°). In this study, d is selected as 60, 120, and 180 mm, respectively. Figures 17–19 present the effect of d on the breakdown pressure, the number of HF, and the HFN pattern. Results show that as d increases, breakdown pressure grows as well and β has a more pronounced effect on the number of HF. When d increases to the largest value (d=180mm), after hydraulic fracturing, the HFN pattern in all the reservoirs is characterized by the step path mode.As shown in Figure17, when β is 0°, the increase in d weakens the effect of joints on stress concentration in the shale reservoir. Therefore, the breakdown pressure grows gradually. However, due to the decline in the number of joints, it is more difficult for the HF to propagate and coalesce. As a result, the number of the HF in the reservoir as a whole diminishes. As β increases to 15° and 30°, the HFN pattern varies from one reservoir to another. This indicates that d exerts a more obvious effect than k (compare Figures 16–19). When d is small (d=60mm), the stress concentration zones around the joint tips overlap. As a result, it is easy for HF to propagate along the joint tip, thus creating a HFN pattern characterized by the tip-to-tip mode. As d reaches its maximum (d=180mm), stress concentration zones around the joint tips are isolated from each other. As a result, after HF initiates at the joint tips, it propagates along the joints with a low strength. Finally, some zones show a HFN pattern characterized by the step path mode. This indicates that, with a small β, an increase in d results in changes to the HFN pattern. The HFN pattern converts from tip-to-tip mode to step path mode.Asβ increases to 45° or 60°, a different d leads to different HFN patterns. This is a combined result of the changes in d and β. Due to the variation of β, the breakdown pressure bottoms out at the same d. In the meantime, crack propagation gains momentum and the joints collapse before they are activated (Figure 19). When β gets larger (β=75° or 90°), the effect of d on the number of HF and the HFN pattern is attenuated. β is the primary reason for the ubiquitous presence of the step path mode in all the shale reservoirs. However, changes in d induce changes to the strength of the shale reservoir. As a result, the breakdown pressure still increases with increasing d. When β is 90° and d is 180 mm, the breakdown pressure in the jointed shale reservoir is over 0.9 times of that in the homogeneous shale reservoir.Therefore, we can conclude that (1) when keeping the joint orientation angle and other parameters unchanged, an increase ind leads to a significant increase in the breakdown pressure and (2) when β falls within the range of 45°~60°, variation in d can cause changes to the number of HF and the HFN pattern. An increase in d can lead to reduction of the proportion of HF induced by joint activation.


### 4.4. Effect of Joint Step Angle on Hydraulic Fracture Network Propagation

Joint step angleγ characterizes the relative position between the parallel joints (Figure 5). To study the effect of γ, we keep other microparameters constant except β (k=0.7, d=120mm, and e=1.2×10−4m. β varies from 0° to 90° at an interval of 15°). In this study, γ is selected as 60°, 90°, and 117°, respectively. Figures 20–22 present the effect of γ on the breakdown pressure, the number of HF, and the HFN pattern. As can be seen from Figure 20, γ has a negligible effect on breakdown pressure. The reason behind this imperceptible effect is that joint strength and joint density are not relevant to γ.As shown in Figure21, when β is 0°, 15°, and 30°, there is a variation in the total number of HF. However, the proportion of different types of HF remains virtually unchanged. This is due to the change in HFN pattern in the reservoir. Following the injection of the fluids into the borehole, rocks surrounding the borehole start to fracture. Though joint strength and joint density remain unchanged, a different β results in changes to the HF initiation locality and propagation direction. Therefore, different HFN patterns are generated (Figure 22).As can be seen from Figure22, when β is equal to or greater than 45°, it has a dominant effect on the HFN pattern. At a constant β, the HFN pattern does not vary with γ any more. Stress still concentrates around the joint tips and within the central section of the joints, irrespective of the changes of γ. The number of HF as well as the proportion of different types of HF are not sensitive to changes of γ. It should be noted that when β is 60°, there is variation in the number of HF (Figure 21). This is due to the low strength of the reservoir when β is 60°. In fact, the reservoir has minimum strength when β is 60°. Therefore, cracks easily initiate, coalesce, and propagate in shale matrices and joints.Therefore, we conclude that (1) whenβ is not greater than 30°, γ and β both have an effect on HFN pattern; (2) when β is not smaller than 30°, β has a dominant effect on the HFN pattern; and (3) at the same β, γ has a negligible effect on breakdown pressure.


### 4.5. Effect of Joint Aperture on Hydraulic Fracture Network Propagation

To study the effect ofe, we keep other microparameters constant except β (k=0.7, d=120mm, and γ=117°, with β changing from 0° to 90° at an interval of 15°). In this study, e is selected as 1.2 × 10−4, 1.2 × 10−3, and 1.2 × 10−2 m, respectively. Figures 23–25 present the effect of e on the breakdown pressure, the number of HF, and the HFN pattern. Results show that as e increases from 1.2 × 10−4 m to 1.2 × 10−3 m, breakdown pressure decreases, but not in an obvious manner, and there are no significant changes to the HFN pattern. Due to the increase in e, the shear HF grows, with its proportion in HF increasing by about 10%.However, ase increases to 1.2 × 10−2 m, significant changes take place in the breakdown pressure, the number of HF, and the HFN pattern. In this scenario, β does not have a dominant effect on breakdown pressure any more. The normalized breakdown pressure falls within a range of 0.15~0.25. In addition, the proportion of tensile HF decreases by about 25%. The above changes are attributable to the opening mode generated after hydraulic fracturing. Due to the increase in e, after injection, most of the fluids flow into the joints and only a slight amount of shale matrices surrounding the injection borehole are fractured. As the injection continues, distant joints start to slip and shift and HF is generated. As for joints that are more distant, it is less likely for them to be activated.Therefore, we can conclude that (1) when keeping the joint orientation angle and other parameters unchanged, a decrease ine leads to an increase in the breakdown pressure and the number of the cracks and (2) changes of e result in the variation in the proportion of different types of hydraulic fractures.


## 4.1. Effect of Joint Orientation Angle on Hydraulic Fracture Network Propagation

To investigate the effect ofβ, we keep other microparameters constant, i.e., k=0.7, d=120mm, γ=117°, and e=1.2×10−4m. Then, we performed numerical simulation with β changing from 0° (horizontal joint) to 90° (vertical joint) at an interval of 15°. The modeling results are shown in Figures 11–13, which separately present the effect of β on the breakdown pressure, the number of HF, and the HFN pattern. The four colors in Figure 13 represent the tensile and shear cracks generated in the shale matrix and the joint. To study to what extent the breakdown pressure in jointed shale reservoir decreases as opposed to that in homogeneous shale reservoir, the breakdown pressure σbp in jointed shale reservoir is normalized to the breakdown pressure σHo in a homogeneous shale reservoir, which is 13.22 MPa.Figure 11
Effect of joint orientation (β) on breakdown pressure.During hydraulic fracturing,σbp reaches its maximum value when β is 90° (vertical joint). While β is 0°, the joint is in a direction parallel with the maximum principal stress. In this case, stress concentrates near the joint tip, and the tensile HF initiates near the joint tips and coalesces with the HF initiated near other joint tips. In the meantime, SHFJP is generated in the joints. Therefore, this HFN pattern is characterized by tip-to-tip mode. As can be seen from Figure 12(a), in this case, the number of HF reaches a peak, amounting to about 2000. When β=15° and 30°, the joint is angled with respect to the direction of the maximum principal stress. In this case, stress concentrates near the joint tip and in the contacting areas between the shale matrix. Moreover, the breakdown pressure declines as shown in Figure 11. The generated HF forms a pattern that is combined by tip-to-tip mode and crossing mode. As β increases to 45° or 60°, stress concentrates within the shale matrix. Under the combined action of the joint orientation and the maximum principal stress, tensile HF propagates in a direction nearly vertical to the joint, forming a HFN pattern characterized by the crossing mode. In this pattern, HF firstly initiates in the internal zones of the joints, followed by the propagation and coalesce of HF between the joints. Finally, HF penetrates the joints. In this case, the breakdown pressure σbp is roughly about 30% of the breakdown pressure σHo. According to the Coulomb criterion, the jointed shale reservoir has its minimum tensile strength when β varies from 45° to 60°. In addition, the previous studies have validated that the breakdown pressure is positively correlated with the tensile strength of the shale reservoir [47]. Therefore, theoretical solution and numerical simulation deliver a consistent result regarding the breakdown pressure, which further validates the correctness of the microparameters used in this study.Figure 12
Effect of joint orientation (β) on hydraulic fracture (HF): (a) number and (b) percentage.
(a)
(b)Figure 13
Effect of joint orientation (β) on the propagation mode of the hydraulic fracture network.Figure 14
Effect of joint persistency (k) on breakdown pressure.Figure 15
Effect of joint persistency (k) on the number and percentage of different types of hydraulic fractures.Figure 16
Effect of joint persistency (k) on the propagation mode of the hydraulic fracture network.Figure 17
Effect of joint spacing (d) on breakdown pressure.Figure 18
Effect of joint spacing (d) on the number and percentage of different types of hydraulic fracture (HF).Figure 19
Effect of joint spacing (d) on the propagation mode of the hydraulic fracture network (HFN).Figure 20
Effect of the joint step angle (γ) on breakdown pressure.Figure 21
Effect of the joint step angle (γ) on the number and percentage of different types of hydraulic fracture (HF).Figure 22
Effect of the joint step angle (γ) on the propagation mode of the hydraulic fracture network (HFN).Figure 23
Effect of the joint aperture (e) on breakdown pressure.Figure 24
Effect of the joint aperture (e) on the number and percentage of different types of hydraulic fracture (HF).Figure 25
Effect of the joint aperture (e) on the propagation mode of the hydraulic fracture network (HFN).Afterβ reaches 75°, stress is less concentrated near the joint tip. Therefore, the breakdown pressure increases and is greater than the breakdown pressure while β is 0°. When HF interacts with the joints, the joints are activated. As a result, HF initiates at the other end of the joint and propagates. However, when β=75° or 90°, cracks do not cover a wide area in the shale reservoir except for areas surrounding the injection borehole. Therefore, the number of HF is smaller than that when β is 0° or 15°. However, it should be noted that the proportion of different types of HF remains unchanged despite the changes in the number of HF with β (Figure 12(b)).The modeling results demonstrate thatβ exerts a significant influence on the breakdown pressure, the number of HF, and the HFN pattern. Therefore, β will also be considered while investigating the effect of other microparameters of the joints.


## 4.2. Effect of Joint Persistency on Hydraulic Fracture Network Propagation

To study the effect ofk, we keep other microparameters constant except β (d=120mm, γ=117°, e=1.2×10−4m, and β changes from 0° to 90° at an interval of 15°). k is changed by varying the length of the joint while the length of the rock bridge is kept unchanged (Lr=60mm). Figures 14–16 present the effect of k on the breakdown pressure, the number of HF, and the HFN pattern. Results show that as k increases, β has a more pronounced effect on the breakdown pressure and the HFN pattern, as shown in Figures 14 and 16. k exerts a more apparent effect on the number of HF when β varies from 45° to 60° as shown in Figure 15.As shown in Figure16, the HF pattern in all the shale reservoirs is characterized by tip-to-tip mode when β is 0°. The number of HF varies slightly from one shale reservoir to another. In this scenario, an increase in k leads to a decrease in breakdown pressure (Figure 14). This is due mainly to the aggravation of stress concentration near the joint tips. In fact, stress has already concentrated around the joint tips even when k is 0.6. Therefore, the increase in k further aggravates stress concentration near the joint tips. This is why the HFN pattern remains unchanged while the breakdown pressure diminishes. As β increases to 15° and 30°, with a small k (k=0.6 or 0.7), the HFN pattern formed is a combination of the tip-to-tip mode and the crossing mode. However, the HFN pattern is still characterized by the tip-to-tip mode when k changes to 0.8. As a larger k results in a greater joint density, stress is more easily concentrated around the joint tips. Therefore, it is easier for HF to initiate at the joint tips and to propagate further (Figure 16). As a result, the breakdown pressure significantly drops (Figure 14). According to the theory of fracture mechanics, greater joint length results in greater stress intensity factor. Therefore, greater stress concentration occurs around the joint tips. This is why as k increases the breakdown pressure drops.Asβ increases to 45° or 60°, the HFN pattern is characterized by the crossing mode. an increase in k leads to an increase in the number of HF. In addition, the HF generated in joints accounts for a greater proportion, with SHFJP in particular. When k reaches the extreme values (k=0.8 and k=0.6), the difference in the number of HF amounts to over 400 and HF propagates in a direction nearly vertical to the joint, as shown in Figure 15. As a result, stress concentration gradually spreads to rock bridge areas, causing a difference in the breakdown pressure and the number of HF. Stress is less concentrated around the joint tips when β is equal to or greater than 75°. Stress is more evenly distributed across the rock bridge area. Therefore, k does not exert an apparent effect on the breakdown pressure. Its effect on the number of HF also weakens. The HFN pattern converts from the crossing mode to the step path mode.Therefore, we can conclude that (1) when keeping the joint orientation angle and other parameters unchanged, an increase ink leads to a decrease in the breakdown pressure and the number of HF; (2) an increase in k causes no perceptible changes to the proportion of HF; and (3) though k has an effect on the HFN pattern, the effect is not as significant as that of β.


## 4.3. Effect of Joint Spacing on Hydraulic Fracture Network Propagation

To study the effect ofd, we keep other microparameters constant except β (k=0.7, γ=117°, e=1.2×10−4m, and β changes from 0° to 90° at an interval of 15°). In this study, d is selected as 60, 120, and 180 mm, respectively. Figures 17–19 present the effect of d on the breakdown pressure, the number of HF, and the HFN pattern. Results show that as d increases, breakdown pressure grows as well and β has a more pronounced effect on the number of HF. When d increases to the largest value (d=180mm), after hydraulic fracturing, the HFN pattern in all the reservoirs is characterized by the step path mode.As shown in Figure17, when β is 0°, the increase in d weakens the effect of joints on stress concentration in the shale reservoir. Therefore, the breakdown pressure grows gradually. However, due to the decline in the number of joints, it is more difficult for the HF to propagate and coalesce. As a result, the number of the HF in the reservoir as a whole diminishes. As β increases to 15° and 30°, the HFN pattern varies from one reservoir to another. This indicates that d exerts a more obvious effect than k (compare Figures 16–19). When d is small (d=60mm), the stress concentration zones around the joint tips overlap. As a result, it is easy for HF to propagate along the joint tip, thus creating a HFN pattern characterized by the tip-to-tip mode. As d reaches its maximum (d=180mm), stress concentration zones around the joint tips are isolated from each other. As a result, after HF initiates at the joint tips, it propagates along the joints with a low strength. Finally, some zones show a HFN pattern characterized by the step path mode. This indicates that, with a small β, an increase in d results in changes to the HFN pattern. The HFN pattern converts from tip-to-tip mode to step path mode.Asβ increases to 45° or 60°, a different d leads to different HFN patterns. This is a combined result of the changes in d and β. Due to the variation of β, the breakdown pressure bottoms out at the same d. In the meantime, crack propagation gains momentum and the joints collapse before they are activated (Figure 19). When β gets larger (β=75° or 90°), the effect of d on the number of HF and the HFN pattern is attenuated. β is the primary reason for the ubiquitous presence of the step path mode in all the shale reservoirs. However, changes in d induce changes to the strength of the shale reservoir. As a result, the breakdown pressure still increases with increasing d. When β is 90° and d is 180 mm, the breakdown pressure in the jointed shale reservoir is over 0.9 times of that in the homogeneous shale reservoir.Therefore, we can conclude that (1) when keeping the joint orientation angle and other parameters unchanged, an increase ind leads to a significant increase in the breakdown pressure and (2) when β falls within the range of 45°~60°, variation in d can cause changes to the number of HF and the HFN pattern. An increase in d can lead to reduction of the proportion of HF induced by joint activation.


## 4.4. Effect of Joint Step Angle on Hydraulic Fracture Network Propagation

Joint step angleγ characterizes the relative position between the parallel joints (Figure 5). To study the effect of γ, we keep other microparameters constant except β (k=0.7, d=120mm, and e=1.2×10−4m. β varies from 0° to 90° at an interval of 15°). In this study, γ is selected as 60°, 90°, and 117°, respectively. Figures 20–22 present the effect of γ on the breakdown pressure, the number of HF, and the HFN pattern. As can be seen from Figure 20, γ has a negligible effect on breakdown pressure. The reason behind this imperceptible effect is that joint strength and joint density are not relevant to γ.As shown in Figure21, when β is 0°, 15°, and 30°, there is a variation in the total number of HF. However, the proportion of different types of HF remains virtually unchanged. This is due to the change in HFN pattern in the reservoir. Following the injection of the fluids into the borehole, rocks surrounding the borehole start to fracture. Though joint strength and joint density remain unchanged, a different β results in changes to the HF initiation locality and propagation direction. Therefore, different HFN patterns are generated (Figure 22).As can be seen from Figure22, when β is equal to or greater than 45°, it has a dominant effect on the HFN pattern. At a constant β, the HFN pattern does not vary with γ any more. Stress still concentrates around the joint tips and within the central section of the joints, irrespective of the changes of γ. The number of HF as well as the proportion of different types of HF are not sensitive to changes of γ. It should be noted that when β is 60°, there is variation in the number of HF (Figure 21). This is due to the low strength of the reservoir when β is 60°. In fact, the reservoir has minimum strength when β is 60°. Therefore, cracks easily initiate, coalesce, and propagate in shale matrices and joints.Therefore, we conclude that (1) whenβ is not greater than 30°, γ and β both have an effect on HFN pattern; (2) when β is not smaller than 30°, β has a dominant effect on the HFN pattern; and (3) at the same β, γ has a negligible effect on breakdown pressure.


## 4.5. Effect of Joint Aperture on Hydraulic Fracture Network Propagation

To study the effect ofe, we keep other microparameters constant except β (k=0.7, d=120mm, and γ=117°, with β changing from 0° to 90° at an interval of 15°). In this study, e is selected as 1.2 × 10−4, 1.2 × 10−3, and 1.2 × 10−2 m, respectively. Figures 23–25 present the effect of e on the breakdown pressure, the number of HF, and the HFN pattern. Results show that as e increases from 1.2 × 10−4 m to 1.2 × 10−3 m, breakdown pressure decreases, but not in an obvious manner, and there are no significant changes to the HFN pattern. Due to the increase in e, the shear HF grows, with its proportion in HF increasing by about 10%.However, ase increases to 1.2 × 10−2 m, significant changes take place in the breakdown pressure, the number of HF, and the HFN pattern. In this scenario, β does not have a dominant effect on breakdown pressure any more. The normalized breakdown pressure falls within a range of 0.15~0.25. In addition, the proportion of tensile HF decreases by about 25%. The above changes are attributable to the opening mode generated after hydraulic fracturing. Due to the increase in e, after injection, most of the fluids flow into the joints and only a slight amount of shale matrices surrounding the injection borehole are fractured. As the injection continues, distant joints start to slip and shift and HF is generated. As for joints that are more distant, it is less likely for them to be activated.Therefore, we can conclude that (1) when keeping the joint orientation angle and other parameters unchanged, a decrease ine leads to an increase in the breakdown pressure and the number of the cracks and (2) changes of e result in the variation in the proportion of different types of hydraulic fractures.


## 5. Conclusions

In this study, a discrete element method-based coupled hydromechanical model is developed in particle flow code to explore the effect of the geometrical parameters of the joints on the breakdown pressure, the number of hydraulic fracture, and the hydraulic fracture network pattern in jointed shale reservoirs. The hydraulic fracture network pattern in the shale reservoir resulting from hydraulic fracturing can be roughly divided into four types, i.e., crossing mode, tip-to-tip mode, step path mode, and opening mode.Whenβ is 0°, except for e=1.2×10−2m, most of the shale reservoirs have a hydraulic fracture network pattern that is characterized by tip-to-tip mode after hydraulic fracturing. When β increases to 15° or 30°, the hydraulic fracture network pattern changes from tip-to-tip mode to crossing mode, as stress concentrates near the joint tip and in the contacting areas between the shale matrices. In this scenario, γ and β both have an effect on the hydraulic fracture network pattern. Due to the change in the hydraulic fracture network pattern, the breakdown pressure and the number of cracks gradually decrease.Whenβ is 45° or 60°, the crossing mode gains dominance, and the hydraulic fracture network pattern depends mainly on β and d. In addition, the magnitude of the breakdown pressure and the number of the cracks reach the lowest level. As β reaches 75° or 90°, the step path mode is ubiquitous, and the magnitude of the breakdown pressure and the number of the cracks both increase. The breakdown pressure increases and is greater than the breakdown pressure when β is 0°. In this scenario, the hydraulic fracture network pattern is largely dependent on β.Under the sameβ, either decrease in k and e or increase in d leads to the rise in the breakdown pressure and the number of the cracks. It is also found that changes in d and e result in the variation in the proportion of different types of hydraulic fractures. In particular, the opening mode of the hydraulic fracture network is observed when e increases to 1.2 × 10−2 m.



---

*Source: 1852604-2018-10-25.xml*