# Analysis of Hydraulic Fracture Propagation Using a Mixed Mode and a Uniaxial Strain Model considering Geomechanical Properties in a Naturally Fractured Shale Reservoir

Authors: Youngho Jang; Gayoung Park; Seoyoon Kwon; Baehyun Min
Journal: Geofluids
 (2020)

---

## Abstract

This study proposes a hydraulic fracture propagation model with a mixed mode comprising opening and sliding modes to describe a complex fracture network in a naturally fractured shale gas formation. We combine the fracture propagation model with the mixed mode and the uniaxial strain model with tectonic impacts to calculate the stress distribution using geomechanical properties. A discrete fracture network is employed to realize the fracture network composed of natural and hydraulic fractures. We compare the fracture propagation behaviours of three cases representing the Barnett, Marcellus, and Eagle Ford shale gas formations. Sensitivity analysis is performed to investigate the effects of the geomechanical properties of the reservoir on the sliding mode’s contribution to the mixed mode. The numerical results highlight the significance of the mixed mode for the accurate assessment of fracture propagation behaviours in shale gas formations with high brittleness.


---

## Body

## 1. Introduction

Hydraulic fracturing using a horizontal well is an essential production technique in shale gas formations [1–5]. Various approaches have been developed to model hydraulic fracturing phenomena to obtain an accurate estimation of shale gas production—for example, extended and generalized finite elements methods (XFEM and GFEM, respectively) [6–9], displacement discontinuity methods (DDM) [10–12], boundary element methods (BEM) [13–15], phase-field formulations [16–22], and discrete fracture network (DFN) [23–25].The propagation of a hydraulic fracture in a naturally fractured reservoir is affected by the existence of natural fractures due to the alteration of the local stress distribution and the friction coefficients on the fracture plane [26]. The propagation yields complex fracture geometry since the hydraulic fracture evolves in a nontrivial way due to natural fractures and geomechanical characteristics of the reservoir [27]. Several experimental studies have investigated the interaction between hydraulic and natural fractures to comprehend fracture propagation [28–32]. Criteria have been proposed to assess whether a hydraulic fracture crosses a natural fracture at the intersection between hydraulic and natural fractures. Blanton [33] suggested a fracture-crossing criterion based on the Griffith theory [34] and performed experiments with varied intersection angles and stresses in order to study the interaction between hydraulic and natural fractures. In addition, Renshaw and Pollard [35] also proposed a fracture-crossing criterion according to orthogonal intersections based on linear elastic fracture mechanics [36]; later, this criterion was extended to nonorthogonal intersections [37]. Several numerical models have been designed to describe complex fracture networks incorporating the interaction between hydraulic and natural fractures. For example, Dahi-Taleghani and Olson [38] developed a two-dimensional XFEM model, and Weng et al. [39] simulated multiple hydraulic fractures and calculated fracture width in a naturally fractured reservoir using DDM coupled with a pseudo-three-dimensional equation.These previous studies of fracture-crossing criteria have mainly focused on mimicking the opening mode in fracture mechanics, allowing a hydraulic fracture to cross a natural fracture directly. However, the tip of a hydraulic fracture may slip along the face of a natural fracture up to the elastic region if both fractures are encountered nonorthogonally. If hydraulic fracture sliding occurs, the next step of hydraulic fracturing is either its propagation into the rock matrix or its arrest within the natural fracture. In this case, not only the opening mode but also the sliding mode becomes remarkable for fracture propagation. The significance of the sliding mode in fracture propagation has been studied analytically and numerically [40–44]. For instance, Jang et al. [44] suggested a fracture-crossing criterion with a mixed mode, which activated both the opening and sliding modes under the assumption that the principal stresses were given. In Jang et al. [44], the fracture-crossing criterion was derived using a linear superposition of both the opening and sliding modes. This criterion was applied to the fracture interaction in the DFN. To the best of our knowledge, researches implementing the mixed mode in the DFN are rare and need further investigation to mimic fracture behaviours more realistically with reflecting various geomechanical properties.It is essential to obtain an accurate stress distribution to estimate in situ stresses at the location of interest [45, 46] since the direction of fracture propagation in porous media depends on the stress distribution around a fracture tip. To evaluate the stress distribution, the critical properties are geomechanical properties such as Poisson’s ratio, Young’s modulus, and the ratio of tectonic strains [21, 47]. As geomechanical properties vary in unconventional shale reservoirs [48], the estimation of hydrocarbon production needs to reflect the variability of the geomechanical properties of such reservoirs. Several studies have assessed the role of the geomechanical properties in fracking [8, 49–52]. For example, a uniaxial strain model with tectonic impacts was developed to calculate in situ principal stresses from geomechanical properties [49]. Brittleness is a geomechanical component that is useful for representing the geomechanical characteristics of shale reservoirs [53]. Rickman et al. [50] introduced the brittleness index (BI), which is calculated using Poisson’s ratio and Young’s modulus. The BI is regarded as an efficient indicator to evaluate the fracturing potential of an oilfield [53]. Note that the success of shale gas production with hydraulic fracturing depends on the brittleness of shale since brittle shales with preexisting natural fractures are hydraulically fracture-prone in the opening and sliding modes [54].The purpose of this study is to show the significance of the sliding mode in fracking based on the mixed mode and DFN to accurately model hydraulic fracture propagation behaviours in a naturally fractured shale reservoir. In addition, various geomechanical properties are reflected to analyse an impact on the fracture interaction. We combine a fracture propagation model with the mixed mode [44] and a uniaxial strain model with tectonic impacts [49] to calculate the distribution of the maximum and minimum horizontal stresses using the geomechanical properties of the reservoir. DFN is employed to realize a fracture network composed of hydraulic and natural fractures. The propagation behaviours of hydraulic fractures with mixed mode are compared with those with opening mode in three major shale gas formations: Barnett, Marcellus, and Eagle Ford. Also, the results of the comparative study are analysed based on BI values of the shale gas formations. Finally, we perform a sensitivity analysis to investigate the effects of the geomechanical properties of the reservoir on the sliding mode’s contribution to the mixed mode.


## 2. Propagation of Hydraulic Fracturing in a Naturally Fractured Reservoir

### 2.1. Fracture Propagation and Modes I, II, and III in Fracture Mechanics

In past studies, fracture propagation has been analysed based on the linear elastic fracture mechanics theory [36] under the assumption that the material is isotropic, linear elastic, and has a quasistatic and isothermal deformation [55, 56]. Researchers have applied the linear elastic fracture mechanics theory to porous media [57–60]. This theory is valid only if an inelastic deformable zone at a fracture tip is small compared to the fracture size [61]. The stress field near the fracture tip can be computed using the theory of elasticity [62].Figure1 illustrates three modes of fracture propagation: Mode I, Mode II, and Mode III [63]. Mode I is the opening mode caused by tensile stress perpendicular to the fracture plane. Mode II is the sliding mode (i.e., the in-plane shear mode), which results from shear stress acting parallel to the fracture surface and perpendicular to the fracture front. Mode III is the tearing mode (i.e., the antiplane shear mode), which is induced by shear stress acting parallel to both the fracture surface and the fracture front. Fracture growth can be modelled using each mode or a combination of modes, which is referred to as a mixed mode [64]. Mode I leads to in-plane tensile propagation. A mixed mode composed of Mode I and Mode II (i.e., in-plane shear) results in bent fractures or sharp kinks, while a mixed mode consisting of Mode I and Mode III (i.e., antiplane shear) causes segmented fracture fronts. As fracture segmentation mostly occurs in unconsolidated or less-consolidated materials [65, 66], this study focuses on comparing the effects of pure Mode I and those of a mixed mode (Mode I and Mode II) on hydraulic fracture propagation in a naturally fractured shale reservoir.Figure 1
Three modes of fracture propagation: Mode I (opening mode), Mode II (sliding mode), and Mode III (tearing mode).


### 2.2. Fracture Crossing between Hydraulic and Natural Fractures

The presence of natural fractures in a shale reservoir affects the propagation of hydraulic fractures [27], thereby leading to the generation of a complex fracture network [38]. There are four interaction types when a hydraulic fracture encounters a natural fracture, as shown in Figure 2 [67]. First, a hydraulic fracture can directly cross a natural fracture without changing its original propagation direction (Figure 2(a)). Second, a hydraulic fracture can join a natural fracture and create a new fracturing path at the tip of the natural fracture (Figure 2(b)). Third, a hydraulic fracture can be diverted along a natural fracture and kink out at a weak point of the natural fracture (Figure 2(c)). Last, a hydraulic fracture can be arrested within a natural fracture (Figure 2(d)). For a hydraulic fracture to cross a natural fracture, two conditions must be satisfied. When the maximum principal stress reaches the tensile strength of porous media, a new fracture is initiated on the opposite side of the natural fracture interface. Also, there should be no shear slippage in the face of the natural fracture [68]. Properties of a material interface, such as frictional resistance and cohesion, control the material’s critical point to prevent slipping or allow fracture crossing [69].Figure 2
Four interaction types between a hydraulic fracture (HF) and a natural fracture (NF). (a) Direct cross. (b) Cross at a natural fracture tip. (c) Cross with an offset. (d) Arrest.
(a)(b)(c)(d)


## 2.1. Fracture Propagation and Modes I, II, and III in Fracture Mechanics

In past studies, fracture propagation has been analysed based on the linear elastic fracture mechanics theory [36] under the assumption that the material is isotropic, linear elastic, and has a quasistatic and isothermal deformation [55, 56]. Researchers have applied the linear elastic fracture mechanics theory to porous media [57–60]. This theory is valid only if an inelastic deformable zone at a fracture tip is small compared to the fracture size [61]. The stress field near the fracture tip can be computed using the theory of elasticity [62].Figure1 illustrates three modes of fracture propagation: Mode I, Mode II, and Mode III [63]. Mode I is the opening mode caused by tensile stress perpendicular to the fracture plane. Mode II is the sliding mode (i.e., the in-plane shear mode), which results from shear stress acting parallel to the fracture surface and perpendicular to the fracture front. Mode III is the tearing mode (i.e., the antiplane shear mode), which is induced by shear stress acting parallel to both the fracture surface and the fracture front. Fracture growth can be modelled using each mode or a combination of modes, which is referred to as a mixed mode [64]. Mode I leads to in-plane tensile propagation. A mixed mode composed of Mode I and Mode II (i.e., in-plane shear) results in bent fractures or sharp kinks, while a mixed mode consisting of Mode I and Mode III (i.e., antiplane shear) causes segmented fracture fronts. As fracture segmentation mostly occurs in unconsolidated or less-consolidated materials [65, 66], this study focuses on comparing the effects of pure Mode I and those of a mixed mode (Mode I and Mode II) on hydraulic fracture propagation in a naturally fractured shale reservoir.Figure 1
Three modes of fracture propagation: Mode I (opening mode), Mode II (sliding mode), and Mode III (tearing mode).


## 2.2. Fracture Crossing between Hydraulic and Natural Fractures

The presence of natural fractures in a shale reservoir affects the propagation of hydraulic fractures [27], thereby leading to the generation of a complex fracture network [38]. There are four interaction types when a hydraulic fracture encounters a natural fracture, as shown in Figure 2 [67]. First, a hydraulic fracture can directly cross a natural fracture without changing its original propagation direction (Figure 2(a)). Second, a hydraulic fracture can join a natural fracture and create a new fracturing path at the tip of the natural fracture (Figure 2(b)). Third, a hydraulic fracture can be diverted along a natural fracture and kink out at a weak point of the natural fracture (Figure 2(c)). Last, a hydraulic fracture can be arrested within a natural fracture (Figure 2(d)). For a hydraulic fracture to cross a natural fracture, two conditions must be satisfied. When the maximum principal stress reaches the tensile strength of porous media, a new fracture is initiated on the opposite side of the natural fracture interface. Also, there should be no shear slippage in the face of the natural fracture [68]. Properties of a material interface, such as frictional resistance and cohesion, control the material’s critical point to prevent slipping or allow fracture crossing [69].Figure 2
Four interaction types between a hydraulic fracture (HF) and a natural fracture (NF). (a) Direct cross. (b) Cross at a natural fracture tip. (c) Cross with an offset. (d) Arrest.
(a)(b)(c)(d)


## 3. Modelling of Hydraulic Fracture Propagation

Section3 describes the workflow of this study for fracture propagation modelling. We combined a previous model [44] and a uniaxial strain model with tectonic stress [49]. By assuming that the overburden stress is higher than the horizontal stress in the formation, this study assumes that fracture propagation occurs on the horizontal plane (i.e., the x-y plane) with a constant fracture height identical to the formation thickness. The porous medium is linearly elastic, isotropic except for principal stresses, and has a quasistatic and isothermal deformation. In reality, the growth of a fracture depends on the injection pressure or injection rate. For simplicity, a hydraulic fracture is allowed to grow up to the maximum fracture half-length set a priori with reference to [70]. In other words, the preset value is the stopping condition for numerical simulation. As this study does not solve the flow of fracturing fluid during fracking, no uncertainties such as fracturing fluid leak are taken into account. Modelling fluid-flow-driven fracturing propagation in the DFN is under investigation as our future work.


### 3.1. Variables for the Modelling of Hydraulic Fracture Propagation

Figure3 shows a schematic of hydraulic fracture propagation crossing a natural fracture with the model variables used in this study. To create multiple transverse hydraulic fractures to enhance hydrocarbon productivity, it is preferable to drill a horizontal well parallel to the direction of the minimum horizontal stress σh [71]. In this case, induced hydraulic fractures are normal to the direction of σh, and this normal vector is identical to the direction of the maximum horizontal stress σH [8]. As shown in Figure 3, in reality, the direction of the horizontal well drilling deviates from the direction of σh due to field constraints. As a result, the fracture initiates nonorthogonally from the horizontal well; in other words, the fracture propagation direction is not parallel to the direction of σH. Let α be the inclination angle between the originally intended fracture propagation direction and σH. Note that the principal horizontal stresses (σh and σH) and α are used to calculate the coupled stress field [σx, σy, τxy]T (Section 3.2.1).Figure 3
Schematic of hydraulic fracture propagation crossing a natural fracture.The next step is to determine the fracture initiation angleθi (Section 3.2.2). If the growing hydraulic fracture meets a preexisting natural fracture, the intersection angle between the two fractures β is obtained (Section 3.2.3). Here, β is the factor that is used to determine whether either fracture-crossing or sliding occurs at the intersection. If β is greater than the minimum fracture-crossing angle, fracture-crossing occurs with the fracture reinitiation angle γ (Section 3.2.4). Otherwise, the hydraulic fracture slides along the natural fracture until reaching a weak point within the natural fracture at which fracture reinitiation is possible.


### 3.2. Workflow for the Modelling of Hydraulic Fracture Propagation

The variables presented in Figure3 are calculated according to the workflow of this study (Figure 4). Detailed calculations can be found in the following subsections of the present section. The first step is to stochastically generate a natural fracture network using the DFN model developed by Jang et al. [72]. The coordinates of natural fractures and the horizontal well are stored. Also, the maximum fracture half-length amax is preset. The in situ stresses are calculated (Section 3.2.1). Secondly, a hydraulic fracture is initiated from a horizontal well using the fracture initiation criterion (Section 3.2.2). If the hydraulic fracture encounters a natural fracture, a decision is made using the fracture-crossing criterion to determine whether the hydraulic fracture crosses the natural fracture (Section 3.2.3). If the cross is activated, the model calculates the fracture reinitiation angle γ and proceeds the fracture propagation with this angle. Otherwise, the hydraulic fracture evolves within the natural fracture (Section 3.2.4). The processes from Section 3.2.2 to Section 3.2.4 are repeated until the sum of the lengths of all connected fractures becomes greater than or equal to amax. Lastly, grid mapping is conducted to transfer the fracture network data of the DFN model (e.g., the coordinates of fractures and the gridblock index of a fracture network). A total connected fracture volume (TCFV) is computed from the DFN model. Oil production can be estimated if a fluid-flow simulator is coupled with the DFN model through grid mapping. In this study, all source codes were written in the C programming language.Figure 4
Flowchart of the developed hydraulic fracture propagation model with mixed mode.


#### 3.2.1. Estimation of In Situ Stresses

For isotropic porous media under uniaxial strain with tectonic stress effects, the minimum horizontal stressσh and the maximum horizontal stress σH are calculated using Equations (1) and (2), respectively [49, 73]:
(1)σh=ν1−νσv−δ1p+δ1p+E1−ν2εh+νεH,(2)σH=ν1−νσv−δ1p+δ1p+E1−ν2εH+νεh,where ν is Poisson’s ratio, σv is the vertical stress, δ1 is Biot’s coefficient, p is the pore pressure, E is Young’s modulus, εh is the minimum tectonic strain, and εH is the maximum tectonic strain. On the right-hand side of each equation, the sum of the first two terms corresponds to the uniaxial strain condition while the third term corresponds to the tectonic stress effects.For a stress field near a fracture tip, the opening mode (Mode I) and sliding mode (Mode II) are linearly superimposed based on the theory of linear elastic fracture mechanics [36]. We let σx and σy denote the normal stresses along the x- and y-directions, respectively. τxy is the shear stress acting on the horizontal plane (i.e., x−y plane). In the Cartesian coordinate system, Equation (3) is used to calculate the coupled stress field [σx, σy, τxy]T by adding the distant in situ stresses (i.e., σH and σh) and the fracture-tip stresses (i.e., the sum of the second and third terms on the right hand side):
(3)σxσyτxy=σHσh0+KI2πrcosθ21−sinθ2sin3θ21+sinθ2sin3θ2sinθ2cos3θ2+KII2πrsinθ2−2+cosθ2cos3θ2cosθ2cos3θ2cotθ21−sinθ2sin3θ2,where r is the radial distance from the fracture tip, θ is the polar angle to the direction of the fracture tip, and KI and KII are the stress intensity factors for implementing the opening mode (Mode I) and sliding mode (Mode II), respectively. The subscript of each stress intensity factor indicates the mode type. Here, tensile stress is positive. KI and KII are computed as follows [74]:
(4)KI=P+δ1p−σHsin2α−σhcos2απa,(5)KII=σH−σh2sin2απa,where P is the fracturing pressure and a is the fracture half-length.


#### 3.2.2. Fracture Initiation Criterion

Once the coupled stress field is computed, the next step is to determine the fracture initiation angleθi. Based on the maximum tangential stress criterion [75], the hydraulic fracture is generated at the angle where the tangential stress is maximum as follows:
(6)θi=2tan−114KIKII−14KIKII2+8.


#### 3.2.3. Fracture Crossing Criterion with Mixed Mode

If the hydraulic fracture encounters a natural fracture nonorthogonally, the intersection angle between the two fracturesβ is calculated. The tip of the hydraulic fracture slides into the elastic region along the face of the natural fracture. Then, the tip propagates through the rock matrix as a function of the stress difference between the tip of the hydraulic fracture and the face of the natural fracture. Therefore, it is essential to consider the mixed mode involving the shear force acting along the face of the natural fracture.According to the linear friction law [35], fracture crossing with the mixed mode occurs at the intersection if the criterion below is satisfied:
(7)τβ1<S0−μσβ2,where τβ1 is the combined shear stress and σβ2 is the combined normal stress, as noted in Figure 3; S0 is the cohesion of the natural fracture face; and μ is the friction coefficient. If τβ1 at the face of the natural fracture is smaller than the effective stress at the tip of the hydraulic fracture S0−μσβ2, the hydraulic fracture crosses the natural fracture at the intersection, driven by the superposition of stresses at the tip of the hydraulic fracture and the interface of the natural fracture. In this study, the minimum fracture-crossing angle refers to the smallest intersection angle that satisfies the inequality (7).Asθ equals β at the intersection, τβ1 and σβ2 are calculated using Equations (8) and (9), respectively:
(8)τβ1=τtip,β1+τr,β1=σH−σh2−K1sinβ2sin3β2−K1tanαtanβ21+cosβ2cos3β2sin2β+K1sinβ2cos3β2+K1tanα−K1tanαsinβ2sin3β2cos2β,(9)σβ2=σtip,β2+σr,β2=σH+σh2+K1−K1tanαsinβ2−σH−σh2−K1sinβ2sin3β2−K1tanαtanβ21+cosβ2cos3β2cos2β−K1sinβ2cos3β2+K1tanα−K1tanαsinβ2sin3β2sin2β,where τtip,β1 is the shear stress at the tip of the hydraulic fracture, τr,β1 is the shear stress on the interface of the natural fracture, σtip,β2 is the normal stress at the tip of the hydraulic fracture, σr,β2 is the normal stress on the interface of the natural fracture, and K1 is the stress level (i.e., substitution variable) required to reinitiate a fracture. Note that K1 is the larger solution of the quadratic polynomial equation below:
(10)AK12+BK1+C=0,where
(11)A=cos2β2−3tan2αsin2β2+2tanαtanβ2cosβ+1,(12)B=−2σH−σh2sinβ2sin3β2−2T0+σH+σh−tanαtanβ22σH−σh2+2σH−σh2cosβ2cos3β2+2T0−σH−σh+cot2β2,(13)C=σH−σh22−T0−σH+σh22,where T0 is the tensile strength of porous media. Equations (8), (9), (10), (11), (12), and (13) are derived from the concept that fractures open when the maximum principal stress reaches the tensile strength. The equations are solved iteratively using subroutines described in Jang et al. [44].


#### 3.2.4. Fracture Reinitiation Criterion

If the hydraulic fracture crosses the natural fracture, the fracture reinitiation angleγ is calculated using the coupled stress field [σx, σy, τxy]T as follows:
(14)γ=12tan−12τxyσx−σy−90°.If fracture crossing does not occur and the sum of the lengths of all connected fractures is smaller thanamax, the fracture is reinitiated at the end of the natural fracture using Equation (6).


### 3.3. Model Validation

The performance of the developed model was tested using two validation cases. First, we evaluated the accuracy of the uniaxial strain model, which was integrated with the proposed model, using reference data obtained from the Cretaceous Travis Peak Formation in the East Texas basin [49]. This formation consists of fine- to very fine-grained sandstones, muddy sandstones, and sandy mudstones. The matrix porosity and matrix permeability are less than 8% and 1.0E-3 μm2, respectively. Table 1 presents the vertical stress σv, minimum horizontal stress σh, pore pressure p, Young’s modulus E, Poisson’s ratio ν, and the depth at which the parameters were measured. Figure 5 shows a scatter plot comparing the measured and predicted σh data; the predicted σh data were obtained by solving Equation (1). For this equation, strain data were acquired from the eastern midcontinent data presented in Dolinar [76]: the minimum and maximum strain values εh and εH were 8.0E-5 MPa and 3.7E-4 MPa, respectively. The prediction results reproduced the measured data sufficiently: the mean absolute deviation (MAD) was 2.920, the coefficient of determination (R2) was 0.901, and the root-mean-square error (RMSE) was 1.369.Table 1
Measured stresses and elastic properties used for the validation of the uniaxial strain model.
Depth [m]σv [MPa]σh [MPa]p [MPa]Depth [m]E [GPa]ν [-]253260.739.727.82533360.27265363.636.329.72653480.12284268.240.030.52835480.18289969.540.732.82898660.14301472.341.934.13012620.15303272.742.134.33030660.15Figure 5
Scatter plot comparing the measured and predicted values of the minimum horizontal stressσh obtained using the uniaxial strain model. The measured data were acquired from Thiercelin and Plumb [49].Second, the numerical results of the developed model were compared with the experimental results [29] to check whether a hydraulic fracture crossed a natural fracture. Six experiments were designed with various intersection angles β and various differences between principal stresses (σH−σh) under a fixed vertical stress σv, as shown in Table 2. The friction coefficient μ was 0.615, the tensile strength of porous media T0 was 4.05 MPa, and the cohesion S0 was zero. The model results agreed with all the experimental results.Table 2
A comparison of experimental and numerical results for the validation of the developed model in terms of fracture crossing at an intersection angleβ.
σv [MPa]σH−σh [MPa]β [°]Experimental resultNumerical result27.66.990CrossingCrossing0.790NoncrossingNoncrossing10.375CrossingCrossing1.475NoncrossingNoncrossing10.345NoncrossingNoncrossing1.445NoncrossingNoncrossing


### 3.4. Brittleness Index

In this study, the BI proposed by Rickman et al. [50] was used to quantitatively compare fracture propagation behaviours in shale gas formations in Section 4. The BI is defined as follows:
(15)BI=12100E−EminEmax−Emin+100ν−νmaxνmin−νmax%,where Emin and Emax are the minimum and maximum Young’s modulus in a given domain, and νmin and νmax are the minimum and maximum Poisson’s ratio in the same domain. The lower the value of ν is and the higher the value of E is, the more brittle the rock is. Thus, the effect of the sliding mode can be analysed using the BI.ν and E used in Equations (1), (2), and (15) are static properties. For practical purposes in shale [77], dynamic ν can be considered equal to static ν. Similar to Rickman et al. [50], this study converted the dynamic Young’s modulus ED into the static Young’s modulus E using Equation (16) [77] as follows:
(16)E=ED3.36742.042.


## 3.1. Variables for the Modelling of Hydraulic Fracture Propagation

Figure3 shows a schematic of hydraulic fracture propagation crossing a natural fracture with the model variables used in this study. To create multiple transverse hydraulic fractures to enhance hydrocarbon productivity, it is preferable to drill a horizontal well parallel to the direction of the minimum horizontal stress σh [71]. In this case, induced hydraulic fractures are normal to the direction of σh, and this normal vector is identical to the direction of the maximum horizontal stress σH [8]. As shown in Figure 3, in reality, the direction of the horizontal well drilling deviates from the direction of σh due to field constraints. As a result, the fracture initiates nonorthogonally from the horizontal well; in other words, the fracture propagation direction is not parallel to the direction of σH. Let α be the inclination angle between the originally intended fracture propagation direction and σH. Note that the principal horizontal stresses (σh and σH) and α are used to calculate the coupled stress field [σx, σy, τxy]T (Section 3.2.1).Figure 3
Schematic of hydraulic fracture propagation crossing a natural fracture.The next step is to determine the fracture initiation angleθi (Section 3.2.2). If the growing hydraulic fracture meets a preexisting natural fracture, the intersection angle between the two fractures β is obtained (Section 3.2.3). Here, β is the factor that is used to determine whether either fracture-crossing or sliding occurs at the intersection. If β is greater than the minimum fracture-crossing angle, fracture-crossing occurs with the fracture reinitiation angle γ (Section 3.2.4). Otherwise, the hydraulic fracture slides along the natural fracture until reaching a weak point within the natural fracture at which fracture reinitiation is possible.


## 3.2. Workflow for the Modelling of Hydraulic Fracture Propagation

The variables presented in Figure3 are calculated according to the workflow of this study (Figure 4). Detailed calculations can be found in the following subsections of the present section. The first step is to stochastically generate a natural fracture network using the DFN model developed by Jang et al. [72]. The coordinates of natural fractures and the horizontal well are stored. Also, the maximum fracture half-length amax is preset. The in situ stresses are calculated (Section 3.2.1). Secondly, a hydraulic fracture is initiated from a horizontal well using the fracture initiation criterion (Section 3.2.2). If the hydraulic fracture encounters a natural fracture, a decision is made using the fracture-crossing criterion to determine whether the hydraulic fracture crosses the natural fracture (Section 3.2.3). If the cross is activated, the model calculates the fracture reinitiation angle γ and proceeds the fracture propagation with this angle. Otherwise, the hydraulic fracture evolves within the natural fracture (Section 3.2.4). The processes from Section 3.2.2 to Section 3.2.4 are repeated until the sum of the lengths of all connected fractures becomes greater than or equal to amax. Lastly, grid mapping is conducted to transfer the fracture network data of the DFN model (e.g., the coordinates of fractures and the gridblock index of a fracture network). A total connected fracture volume (TCFV) is computed from the DFN model. Oil production can be estimated if a fluid-flow simulator is coupled with the DFN model through grid mapping. In this study, all source codes were written in the C programming language.Figure 4
Flowchart of the developed hydraulic fracture propagation model with mixed mode.


### 3.2.1. Estimation of In Situ Stresses

For isotropic porous media under uniaxial strain with tectonic stress effects, the minimum horizontal stressσh and the maximum horizontal stress σH are calculated using Equations (1) and (2), respectively [49, 73]:
(1)σh=ν1−νσv−δ1p+δ1p+E1−ν2εh+νεH,(2)σH=ν1−νσv−δ1p+δ1p+E1−ν2εH+νεh,where ν is Poisson’s ratio, σv is the vertical stress, δ1 is Biot’s coefficient, p is the pore pressure, E is Young’s modulus, εh is the minimum tectonic strain, and εH is the maximum tectonic strain. On the right-hand side of each equation, the sum of the first two terms corresponds to the uniaxial strain condition while the third term corresponds to the tectonic stress effects.For a stress field near a fracture tip, the opening mode (Mode I) and sliding mode (Mode II) are linearly superimposed based on the theory of linear elastic fracture mechanics [36]. We let σx and σy denote the normal stresses along the x- and y-directions, respectively. τxy is the shear stress acting on the horizontal plane (i.e., x−y plane). In the Cartesian coordinate system, Equation (3) is used to calculate the coupled stress field [σx, σy, τxy]T by adding the distant in situ stresses (i.e., σH and σh) and the fracture-tip stresses (i.e., the sum of the second and third terms on the right hand side):
(3)σxσyτxy=σHσh0+KI2πrcosθ21−sinθ2sin3θ21+sinθ2sin3θ2sinθ2cos3θ2+KII2πrsinθ2−2+cosθ2cos3θ2cosθ2cos3θ2cotθ21−sinθ2sin3θ2,where r is the radial distance from the fracture tip, θ is the polar angle to the direction of the fracture tip, and KI and KII are the stress intensity factors for implementing the opening mode (Mode I) and sliding mode (Mode II), respectively. The subscript of each stress intensity factor indicates the mode type. Here, tensile stress is positive. KI and KII are computed as follows [74]:
(4)KI=P+δ1p−σHsin2α−σhcos2απa,(5)KII=σH−σh2sin2απa,where P is the fracturing pressure and a is the fracture half-length.


### 3.2.2. Fracture Initiation Criterion

Once the coupled stress field is computed, the next step is to determine the fracture initiation angleθi. Based on the maximum tangential stress criterion [75], the hydraulic fracture is generated at the angle where the tangential stress is maximum as follows:
(6)θi=2tan−114KIKII−14KIKII2+8.


### 3.2.3. Fracture Crossing Criterion with Mixed Mode

If the hydraulic fracture encounters a natural fracture nonorthogonally, the intersection angle between the two fracturesβ is calculated. The tip of the hydraulic fracture slides into the elastic region along the face of the natural fracture. Then, the tip propagates through the rock matrix as a function of the stress difference between the tip of the hydraulic fracture and the face of the natural fracture. Therefore, it is essential to consider the mixed mode involving the shear force acting along the face of the natural fracture.According to the linear friction law [35], fracture crossing with the mixed mode occurs at the intersection if the criterion below is satisfied:
(7)τβ1<S0−μσβ2,where τβ1 is the combined shear stress and σβ2 is the combined normal stress, as noted in Figure 3; S0 is the cohesion of the natural fracture face; and μ is the friction coefficient. If τβ1 at the face of the natural fracture is smaller than the effective stress at the tip of the hydraulic fracture S0−μσβ2, the hydraulic fracture crosses the natural fracture at the intersection, driven by the superposition of stresses at the tip of the hydraulic fracture and the interface of the natural fracture. In this study, the minimum fracture-crossing angle refers to the smallest intersection angle that satisfies the inequality (7).Asθ equals β at the intersection, τβ1 and σβ2 are calculated using Equations (8) and (9), respectively:
(8)τβ1=τtip,β1+τr,β1=σH−σh2−K1sinβ2sin3β2−K1tanαtanβ21+cosβ2cos3β2sin2β+K1sinβ2cos3β2+K1tanα−K1tanαsinβ2sin3β2cos2β,(9)σβ2=σtip,β2+σr,β2=σH+σh2+K1−K1tanαsinβ2−σH−σh2−K1sinβ2sin3β2−K1tanαtanβ21+cosβ2cos3β2cos2β−K1sinβ2cos3β2+K1tanα−K1tanαsinβ2sin3β2sin2β,where τtip,β1 is the shear stress at the tip of the hydraulic fracture, τr,β1 is the shear stress on the interface of the natural fracture, σtip,β2 is the normal stress at the tip of the hydraulic fracture, σr,β2 is the normal stress on the interface of the natural fracture, and K1 is the stress level (i.e., substitution variable) required to reinitiate a fracture. Note that K1 is the larger solution of the quadratic polynomial equation below:
(10)AK12+BK1+C=0,where
(11)A=cos2β2−3tan2αsin2β2+2tanαtanβ2cosβ+1,(12)B=−2σH−σh2sinβ2sin3β2−2T0+σH+σh−tanαtanβ22σH−σh2+2σH−σh2cosβ2cos3β2+2T0−σH−σh+cot2β2,(13)C=σH−σh22−T0−σH+σh22,where T0 is the tensile strength of porous media. Equations (8), (9), (10), (11), (12), and (13) are derived from the concept that fractures open when the maximum principal stress reaches the tensile strength. The equations are solved iteratively using subroutines described in Jang et al. [44].


### 3.2.4. Fracture Reinitiation Criterion

If the hydraulic fracture crosses the natural fracture, the fracture reinitiation angleγ is calculated using the coupled stress field [σx, σy, τxy]T as follows:
(14)γ=12tan−12τxyσx−σy−90°.If fracture crossing does not occur and the sum of the lengths of all connected fractures is smaller thanamax, the fracture is reinitiated at the end of the natural fracture using Equation (6).


## 3.2.1. Estimation of In Situ Stresses

For isotropic porous media under uniaxial strain with tectonic stress effects, the minimum horizontal stressσh and the maximum horizontal stress σH are calculated using Equations (1) and (2), respectively [49, 73]:
(1)σh=ν1−νσv−δ1p+δ1p+E1−ν2εh+νεH,(2)σH=ν1−νσv−δ1p+δ1p+E1−ν2εH+νεh,where ν is Poisson’s ratio, σv is the vertical stress, δ1 is Biot’s coefficient, p is the pore pressure, E is Young’s modulus, εh is the minimum tectonic strain, and εH is the maximum tectonic strain. On the right-hand side of each equation, the sum of the first two terms corresponds to the uniaxial strain condition while the third term corresponds to the tectonic stress effects.For a stress field near a fracture tip, the opening mode (Mode I) and sliding mode (Mode II) are linearly superimposed based on the theory of linear elastic fracture mechanics [36]. We let σx and σy denote the normal stresses along the x- and y-directions, respectively. τxy is the shear stress acting on the horizontal plane (i.e., x−y plane). In the Cartesian coordinate system, Equation (3) is used to calculate the coupled stress field [σx, σy, τxy]T by adding the distant in situ stresses (i.e., σH and σh) and the fracture-tip stresses (i.e., the sum of the second and third terms on the right hand side):
(3)σxσyτxy=σHσh0+KI2πrcosθ21−sinθ2sin3θ21+sinθ2sin3θ2sinθ2cos3θ2+KII2πrsinθ2−2+cosθ2cos3θ2cosθ2cos3θ2cotθ21−sinθ2sin3θ2,where r is the radial distance from the fracture tip, θ is the polar angle to the direction of the fracture tip, and KI and KII are the stress intensity factors for implementing the opening mode (Mode I) and sliding mode (Mode II), respectively. The subscript of each stress intensity factor indicates the mode type. Here, tensile stress is positive. KI and KII are computed as follows [74]:
(4)KI=P+δ1p−σHsin2α−σhcos2απa,(5)KII=σH−σh2sin2απa,where P is the fracturing pressure and a is the fracture half-length.


## 3.2.2. Fracture Initiation Criterion

Once the coupled stress field is computed, the next step is to determine the fracture initiation angleθi. Based on the maximum tangential stress criterion [75], the hydraulic fracture is generated at the angle where the tangential stress is maximum as follows:
(6)θi=2tan−114KIKII−14KIKII2+8.


## 3.2.3. Fracture Crossing Criterion with Mixed Mode

If the hydraulic fracture encounters a natural fracture nonorthogonally, the intersection angle between the two fracturesβ is calculated. The tip of the hydraulic fracture slides into the elastic region along the face of the natural fracture. Then, the tip propagates through the rock matrix as a function of the stress difference between the tip of the hydraulic fracture and the face of the natural fracture. Therefore, it is essential to consider the mixed mode involving the shear force acting along the face of the natural fracture.According to the linear friction law [35], fracture crossing with the mixed mode occurs at the intersection if the criterion below is satisfied:
(7)τβ1<S0−μσβ2,where τβ1 is the combined shear stress and σβ2 is the combined normal stress, as noted in Figure 3; S0 is the cohesion of the natural fracture face; and μ is the friction coefficient. If τβ1 at the face of the natural fracture is smaller than the effective stress at the tip of the hydraulic fracture S0−μσβ2, the hydraulic fracture crosses the natural fracture at the intersection, driven by the superposition of stresses at the tip of the hydraulic fracture and the interface of the natural fracture. In this study, the minimum fracture-crossing angle refers to the smallest intersection angle that satisfies the inequality (7).Asθ equals β at the intersection, τβ1 and σβ2 are calculated using Equations (8) and (9), respectively:
(8)τβ1=τtip,β1+τr,β1=σH−σh2−K1sinβ2sin3β2−K1tanαtanβ21+cosβ2cos3β2sin2β+K1sinβ2cos3β2+K1tanα−K1tanαsinβ2sin3β2cos2β,(9)σβ2=σtip,β2+σr,β2=σH+σh2+K1−K1tanαsinβ2−σH−σh2−K1sinβ2sin3β2−K1tanαtanβ21+cosβ2cos3β2cos2β−K1sinβ2cos3β2+K1tanα−K1tanαsinβ2sin3β2sin2β,where τtip,β1 is the shear stress at the tip of the hydraulic fracture, τr,β1 is the shear stress on the interface of the natural fracture, σtip,β2 is the normal stress at the tip of the hydraulic fracture, σr,β2 is the normal stress on the interface of the natural fracture, and K1 is the stress level (i.e., substitution variable) required to reinitiate a fracture. Note that K1 is the larger solution of the quadratic polynomial equation below:
(10)AK12+BK1+C=0,where
(11)A=cos2β2−3tan2αsin2β2+2tanαtanβ2cosβ+1,(12)B=−2σH−σh2sinβ2sin3β2−2T0+σH+σh−tanαtanβ22σH−σh2+2σH−σh2cosβ2cos3β2+2T0−σH−σh+cot2β2,(13)C=σH−σh22−T0−σH+σh22,where T0 is the tensile strength of porous media. Equations (8), (9), (10), (11), (12), and (13) are derived from the concept that fractures open when the maximum principal stress reaches the tensile strength. The equations are solved iteratively using subroutines described in Jang et al. [44].


## 3.2.4. Fracture Reinitiation Criterion

If the hydraulic fracture crosses the natural fracture, the fracture reinitiation angleγ is calculated using the coupled stress field [σx, σy, τxy]T as follows:
(14)γ=12tan−12τxyσx−σy−90°.If fracture crossing does not occur and the sum of the lengths of all connected fractures is smaller thanamax, the fracture is reinitiated at the end of the natural fracture using Equation (6).


## 3.3. Model Validation

The performance of the developed model was tested using two validation cases. First, we evaluated the accuracy of the uniaxial strain model, which was integrated with the proposed model, using reference data obtained from the Cretaceous Travis Peak Formation in the East Texas basin [49]. This formation consists of fine- to very fine-grained sandstones, muddy sandstones, and sandy mudstones. The matrix porosity and matrix permeability are less than 8% and 1.0E-3 μm2, respectively. Table 1 presents the vertical stress σv, minimum horizontal stress σh, pore pressure p, Young’s modulus E, Poisson’s ratio ν, and the depth at which the parameters were measured. Figure 5 shows a scatter plot comparing the measured and predicted σh data; the predicted σh data were obtained by solving Equation (1). For this equation, strain data were acquired from the eastern midcontinent data presented in Dolinar [76]: the minimum and maximum strain values εh and εH were 8.0E-5 MPa and 3.7E-4 MPa, respectively. The prediction results reproduced the measured data sufficiently: the mean absolute deviation (MAD) was 2.920, the coefficient of determination (R2) was 0.901, and the root-mean-square error (RMSE) was 1.369.Table 1
Measured stresses and elastic properties used for the validation of the uniaxial strain model.
Depth [m]σv [MPa]σh [MPa]p [MPa]Depth [m]E [GPa]ν [-]253260.739.727.82533360.27265363.636.329.72653480.12284268.240.030.52835480.18289969.540.732.82898660.14301472.341.934.13012620.15303272.742.134.33030660.15Figure 5
Scatter plot comparing the measured and predicted values of the minimum horizontal stressσh obtained using the uniaxial strain model. The measured data were acquired from Thiercelin and Plumb [49].Second, the numerical results of the developed model were compared with the experimental results [29] to check whether a hydraulic fracture crossed a natural fracture. Six experiments were designed with various intersection angles β and various differences between principal stresses (σH−σh) under a fixed vertical stress σv, as shown in Table 2. The friction coefficient μ was 0.615, the tensile strength of porous media T0 was 4.05 MPa, and the cohesion S0 was zero. The model results agreed with all the experimental results.Table 2
A comparison of experimental and numerical results for the validation of the developed model in terms of fracture crossing at an intersection angleβ.
σv [MPa]σH−σh [MPa]β [°]Experimental resultNumerical result27.66.990CrossingCrossing0.790NoncrossingNoncrossing10.375CrossingCrossing1.475NoncrossingNoncrossing10.345NoncrossingNoncrossing1.445NoncrossingNoncrossing


## 3.4. Brittleness Index

In this study, the BI proposed by Rickman et al. [50] was used to quantitatively compare fracture propagation behaviours in shale gas formations in Section 4. The BI is defined as follows:
(15)BI=12100E−EminEmax−Emin+100ν−νmaxνmin−νmax%,where Emin and Emax are the minimum and maximum Young’s modulus in a given domain, and νmin and νmax are the minimum and maximum Poisson’s ratio in the same domain. The lower the value of ν is and the higher the value of E is, the more brittle the rock is. Thus, the effect of the sliding mode can be analysed using the BI.ν and E used in Equations (1), (2), and (15) are static properties. For practical purposes in shale [77], dynamic ν can be considered equal to static ν. Similar to Rickman et al. [50], this study converted the dynamic Young’s modulus ED into the static Young’s modulus E using Equation (16) [77] as follows:
(16)E=ED3.36742.042.


## 4. Results and Discussion

In order to analyse the effects of the mixed mode on fracture propagation behaviours, we applied the developed model to fractured reservoirs mimicking the Marcellus, Barnett, and Eagle Ford shale formations. In addition, the results of the comparative study are analysed based on the BI values of the shale gas formations. Finally, we performed a sensitivity analysis to investigate the effects of the geomechanical properties of the reservoir on the contribution of the sliding mode to the mixed mode.


### 4.1. Experimental Setting

Figure6 depicts a base system of DFN, and Table 3 summarizes the input data (i.e., reservoir properties and well configuration) of the base system. This system has σh in the W-E direction. A four-stage (A, B, C, and D) hydraulic fracturing was completed for a horizontal well with an α of 15°. The length of the horizontal well is 473 m, amax is 61 m, μ is 0.7, and δ1 is 1. Table 4 presents the geomechanical properties of the Marcellus, Barnett, and Eagle Ford shale formations. For ν and ED, the range and the mode with the highest frequency were acquired from Yenugu [78], resulting in three different BI values. Since these properties are dynamic properties obtained from well logging, ED values were converted into E values using Equation (16). Representative values of εH and εh were obtained from Dolinar [76]. For a hydraulic fracture to initiate on the opposite side of a natural fracture face, the maximum principal stress must be equal to the rock’s tensile strength T0. Here, the values of T0 were constants of 3, 4, and 6 MPa for the Marcellus, Barnett, and Eagle Ford formations, respectively.Figure 6
Horizontal well with four fracturing stages (A, B, C, and D) in a naturally fractured reservoir.Table 3
Input data used for fracture propagation analysis.
ParameterValueParameterValueReservoir volume [m3]762×366×9Top depth [m]762Length of horizontal well [m]473Number of fracturing stages4amax [m]61Fracture permeability [μm2]1μ [dimensionless]0.7δ1 [dimensionless]1Table 4
Geomechanical properties of the Marcellus, Barnett, and Eagle Ford shale formations.
Shale formationStatistical parameterν [-]ED [GPa]εH/εh [-]T0 [MPa]MarcellusRange0.11–0.2816–281.76–2.823Mode0.17182.17BarnettRange0.05–0.2627–483.00–6.634Mode0.2454.60Eagle FordRange0.22–0.3734–623.00–6.636Mode0.3354.60Although this study utilises field data acquired from references (Table4), the domain shown in Figure 6 is synthetic. For a fair comparison, all the three shale formation cases use the same domain. As a consequence, the same amax is given to each case.


### 4.2. Comparison of Fracture Propagation Behaviours in Three Shale Formations

Using the developed model with the mixed mode, we compare fracture propagation behaviours in the three shale formations using the mode values of the geomechanical properties presented in Table4. The inclination angle α of 15° yielded a fracture initiation angle θi of 7.5°. The simulation results with the mixed mode were also compared with the simulation results with the opening mode to highlight the effects of the mixed mode on fracture propagation.Figure7 shows the simulation results for the evolved fracture networks with and without the sliding mode for the three shale formations. All shale formations were given an identical initial natural fracture network. We monitored whether the hydraulic fracture initiated from Stage A crossed the natural fracture, where their intersection is pointed out in a dashed circle in Figure 7. At the intersection in the dashed circle, β was 69°. Then, we calculated the minimum fracture-crossing angles using the inequality (7). In Figures 7(a)–7(f), the calculated minimum fracture-crossing angles were 58°, 73°, 85°, 53°, 65°, and 80°, respectively. Based on these results, adopting either an opening mode or a mixed mode resulted in fracture-crossing at the intersection for Marcellus and Eagle Ford cases. Interestingly, for the Barnett case with the mixed mode (Figure 7(b)), the hydraulic fracture did not cross the natural fracture as β of 69° was less than the minimum fracture-crossing angle of 73°. The hydraulic fracture grew at the tip of the connected natural fracture. However, when adopting the opening mode, the hydraulic fracture crossed the natural fracture at their intersection as the minimum fracture-crossing angle decreased to 65° (Figure 7(e)). This result highlights the significance of implementing the mixed mode for the accurate modelling of fracture propagation.Figure 7
Simulated fracture network systems in three shale formations. (a–c) show the results generated using the mixed mode model for the Marcellus, Barnett, and Eagle Ford formations, respectively; (d–f) show the results generated using the opening-mode model for the same three formations, respectively.
(a)(b)(c)(d)(e)(f)Based on the BI values of the three shale formations, we examined the effects of fracturing modes on the TCFV. When using the mode values forν and E, the BI values were 57.14%, 40.69%, and 25.12% for the Barnett, Marcellus, and Eagle Ford formations, respectively. Here, TCFV means the total fracture volume connected to the horizontal well. Figure 8 summarizes the TCFV estimates with and without the sliding mode. The TCFV values obtained using the mixed mode were 4,815, 4,107, and 2,764 m3 for the Barnett, Marcellus, and Eagle Ford formations, respectively. For all the formations, ignoring the sliding mode led to the overestimation of the TCFV values. The TCFV values obtained using the opening mode were 5,903, 4,581, and 2,908 m3 for the Barnett, Marcellus, and Eagle Ford formations, respectively. The greater the BI value was, the greater the TCFV was since fracture-crossing occurred easily. In other words, if the value of BI was relatively high such as in the Barnett shale formation, the fracture network resulted in a large deformation longitudinally rather than transversely. On the other hand, if the value of BI was relatively small such as in the Eagle Ford shale formation, a large deformation in the transverse direction was captured. Meanwhile, the difference in TCFV with and without the sliding mode increased as the BI value increased. This is because the influence of the sliding mode on the fracture propagation increases with increasing BI. In the next section, a sensitivity analysis was performed to examine if the above results could be extended to other shale reservoirs.Figure 8
Total connected fracture volume (TCFV) obtained using mixed and opening modes.


### 4.3. Effects of Geomechanical Properties on Fracture Propagation

Sensitivity analysis was performed to investigate the effects of three reservoir geomechanical properties (i.e.,ν, E, and εH/εh) on fracture propagation behaviours. In order to examine the effects of the uniaxial strain model (see Equations (1) and (2)) on the minimum fracture-crossing angle, ν, E, and εH/εh were selected as sensitivity analysis factors.


#### 4.3.1. Effects of Poisson’s Ratio

Figure9 depicts the minimum fracture-crossing angles obtained with various values of Poisson’s ratio ν in the range given in Table 4, where ν indicates the distribution of the rock deformation direction. The increase in ν increased the strain rate, thereby decreasing the minimum fracture-crossing angle. In other words, fracture-crossing is more likely to occur as ν increases.Figure 9
The effect of Poisson’s ratio on the minimum fracture-crossing angle.As shown in Figures9 and 10, all three shale formations showed a similar trend in terms of the minimum fracture-crossing angle and relative difference: the difference between the two calculated minimum fracture-crossing angles of the mixed mode and the opening mode divided by the minimum fracture-crossing angle with the mixed mode. However, their magnitudes were different. Figure 10 shows the relationship between ν and the relative difference. The relative difference increased as ν increased. For the Barnett shale, the relative difference was 1.7% when ν was at the lower limit of 0.05, and the relative difference increased to 13.2% when ν was at the upper limit of 0.26. This trend was also observed in the other two shale formations. These results were due to the interaction between the shear modulus and ν: because the shear modulus is inversely proportional to ν, the shear modulus decreases as ν increases. As a result, in the mixed mode, the influence of the sliding mode on the interaction between fractures became more significant as ν increased.Figure 10
The effect of Poisson’s ratio on the relative difference.


#### 4.3.2. Effects of Young’s Modulus

We also analysed the fracture propagation behaviour considering the Young’s modulusE, which indicates the stiffness of the rock formation. Figure 11 shows the minimum fracture-crossing angles obtained with various values of E for the three shale formations. As E increased, the minimum fracture-crossing angle increased. As a result, fracture-crossing was difficult to occur.Figure 11
The effect of Young’s modulus on the minimum fracture-crossing angle.Figure12 shows the trend of the relative difference. For example, in the Barnett shale, when E was at the lower limit of 27 GPa, the relative difference was 9.8%. At the upper limit of 48 GPa, the relative difference decreased to 2.7%. As the shear modulus is proportional to E, an increase in E increases the shear modulus. As a result, for a fixed stress value, the influence of the sliding mode increased when E decreased. As a result, in the mixed mode, the influence of the sliding mode on the interaction between fractures became more significant as E decreased.Figure 12
The effect of Young’s modulus on the relative difference.


#### 4.3.3. Effects of Tectonic Strain Anisotropy

Moreover, we investigated the effects of tectonic strain anisotropy on fracture propagation since shale formations typically have a large strain anisotropy, which affects the direction of fracture propagation. Figure13 shows the minimum fracture-crossing angle obtained for various values of the strain ratio εH/εh. As the strain ratio increased, the minimum fracture-crossing angle decreased and fracture-crossing was more likely to occur. As shown in Figure 14, the relative difference increased as εH/εh increased. In the Barnett shale, when εH/εh was 3.00, the relative difference was 1.3%. The relative difference increased to 9.6% when εH/εh was 6.63. This trend existed because the rock became brittle through the increase in the shear stress and shear strain as εH/εh increased.Figure 13
The effect of the tectonic strain ratio on the minimum fracture-crossing angle.Figure 14
The effect of tectonic strain ratio on the relative difference.


## 4.1. Experimental Setting

Figure6 depicts a base system of DFN, and Table 3 summarizes the input data (i.e., reservoir properties and well configuration) of the base system. This system has σh in the W-E direction. A four-stage (A, B, C, and D) hydraulic fracturing was completed for a horizontal well with an α of 15°. The length of the horizontal well is 473 m, amax is 61 m, μ is 0.7, and δ1 is 1. Table 4 presents the geomechanical properties of the Marcellus, Barnett, and Eagle Ford shale formations. For ν and ED, the range and the mode with the highest frequency were acquired from Yenugu [78], resulting in three different BI values. Since these properties are dynamic properties obtained from well logging, ED values were converted into E values using Equation (16). Representative values of εH and εh were obtained from Dolinar [76]. For a hydraulic fracture to initiate on the opposite side of a natural fracture face, the maximum principal stress must be equal to the rock’s tensile strength T0. Here, the values of T0 were constants of 3, 4, and 6 MPa for the Marcellus, Barnett, and Eagle Ford formations, respectively.Figure 6
Horizontal well with four fracturing stages (A, B, C, and D) in a naturally fractured reservoir.Table 3
Input data used for fracture propagation analysis.
ParameterValueParameterValueReservoir volume [m3]762×366×9Top depth [m]762Length of horizontal well [m]473Number of fracturing stages4amax [m]61Fracture permeability [μm2]1μ [dimensionless]0.7δ1 [dimensionless]1Table 4
Geomechanical properties of the Marcellus, Barnett, and Eagle Ford shale formations.
Shale formationStatistical parameterν [-]ED [GPa]εH/εh [-]T0 [MPa]MarcellusRange0.11–0.2816–281.76–2.823Mode0.17182.17BarnettRange0.05–0.2627–483.00–6.634Mode0.2454.60Eagle FordRange0.22–0.3734–623.00–6.636Mode0.3354.60Although this study utilises field data acquired from references (Table4), the domain shown in Figure 6 is synthetic. For a fair comparison, all the three shale formation cases use the same domain. As a consequence, the same amax is given to each case.


## 4.2. Comparison of Fracture Propagation Behaviours in Three Shale Formations

Using the developed model with the mixed mode, we compare fracture propagation behaviours in the three shale formations using the mode values of the geomechanical properties presented in Table4. The inclination angle α of 15° yielded a fracture initiation angle θi of 7.5°. The simulation results with the mixed mode were also compared with the simulation results with the opening mode to highlight the effects of the mixed mode on fracture propagation.Figure7 shows the simulation results for the evolved fracture networks with and without the sliding mode for the three shale formations. All shale formations were given an identical initial natural fracture network. We monitored whether the hydraulic fracture initiated from Stage A crossed the natural fracture, where their intersection is pointed out in a dashed circle in Figure 7. At the intersection in the dashed circle, β was 69°. Then, we calculated the minimum fracture-crossing angles using the inequality (7). In Figures 7(a)–7(f), the calculated minimum fracture-crossing angles were 58°, 73°, 85°, 53°, 65°, and 80°, respectively. Based on these results, adopting either an opening mode or a mixed mode resulted in fracture-crossing at the intersection for Marcellus and Eagle Ford cases. Interestingly, for the Barnett case with the mixed mode (Figure 7(b)), the hydraulic fracture did not cross the natural fracture as β of 69° was less than the minimum fracture-crossing angle of 73°. The hydraulic fracture grew at the tip of the connected natural fracture. However, when adopting the opening mode, the hydraulic fracture crossed the natural fracture at their intersection as the minimum fracture-crossing angle decreased to 65° (Figure 7(e)). This result highlights the significance of implementing the mixed mode for the accurate modelling of fracture propagation.Figure 7
Simulated fracture network systems in three shale formations. (a–c) show the results generated using the mixed mode model for the Marcellus, Barnett, and Eagle Ford formations, respectively; (d–f) show the results generated using the opening-mode model for the same three formations, respectively.
(a)(b)(c)(d)(e)(f)Based on the BI values of the three shale formations, we examined the effects of fracturing modes on the TCFV. When using the mode values forν and E, the BI values were 57.14%, 40.69%, and 25.12% for the Barnett, Marcellus, and Eagle Ford formations, respectively. Here, TCFV means the total fracture volume connected to the horizontal well. Figure 8 summarizes the TCFV estimates with and without the sliding mode. The TCFV values obtained using the mixed mode were 4,815, 4,107, and 2,764 m3 for the Barnett, Marcellus, and Eagle Ford formations, respectively. For all the formations, ignoring the sliding mode led to the overestimation of the TCFV values. The TCFV values obtained using the opening mode were 5,903, 4,581, and 2,908 m3 for the Barnett, Marcellus, and Eagle Ford formations, respectively. The greater the BI value was, the greater the TCFV was since fracture-crossing occurred easily. In other words, if the value of BI was relatively high such as in the Barnett shale formation, the fracture network resulted in a large deformation longitudinally rather than transversely. On the other hand, if the value of BI was relatively small such as in the Eagle Ford shale formation, a large deformation in the transverse direction was captured. Meanwhile, the difference in TCFV with and without the sliding mode increased as the BI value increased. This is because the influence of the sliding mode on the fracture propagation increases with increasing BI. In the next section, a sensitivity analysis was performed to examine if the above results could be extended to other shale reservoirs.Figure 8
Total connected fracture volume (TCFV) obtained using mixed and opening modes.


## 4.3. Effects of Geomechanical Properties on Fracture Propagation

Sensitivity analysis was performed to investigate the effects of three reservoir geomechanical properties (i.e.,ν, E, and εH/εh) on fracture propagation behaviours. In order to examine the effects of the uniaxial strain model (see Equations (1) and (2)) on the minimum fracture-crossing angle, ν, E, and εH/εh were selected as sensitivity analysis factors.


### 4.3.1. Effects of Poisson’s Ratio

Figure9 depicts the minimum fracture-crossing angles obtained with various values of Poisson’s ratio ν in the range given in Table 4, where ν indicates the distribution of the rock deformation direction. The increase in ν increased the strain rate, thereby decreasing the minimum fracture-crossing angle. In other words, fracture-crossing is more likely to occur as ν increases.Figure 9
The effect of Poisson’s ratio on the minimum fracture-crossing angle.As shown in Figures9 and 10, all three shale formations showed a similar trend in terms of the minimum fracture-crossing angle and relative difference: the difference between the two calculated minimum fracture-crossing angles of the mixed mode and the opening mode divided by the minimum fracture-crossing angle with the mixed mode. However, their magnitudes were different. Figure 10 shows the relationship between ν and the relative difference. The relative difference increased as ν increased. For the Barnett shale, the relative difference was 1.7% when ν was at the lower limit of 0.05, and the relative difference increased to 13.2% when ν was at the upper limit of 0.26. This trend was also observed in the other two shale formations. These results were due to the interaction between the shear modulus and ν: because the shear modulus is inversely proportional to ν, the shear modulus decreases as ν increases. As a result, in the mixed mode, the influence of the sliding mode on the interaction between fractures became more significant as ν increased.Figure 10
The effect of Poisson’s ratio on the relative difference.


### 4.3.2. Effects of Young’s Modulus

We also analysed the fracture propagation behaviour considering the Young’s modulusE, which indicates the stiffness of the rock formation. Figure 11 shows the minimum fracture-crossing angles obtained with various values of E for the three shale formations. As E increased, the minimum fracture-crossing angle increased. As a result, fracture-crossing was difficult to occur.Figure 11
The effect of Young’s modulus on the minimum fracture-crossing angle.Figure12 shows the trend of the relative difference. For example, in the Barnett shale, when E was at the lower limit of 27 GPa, the relative difference was 9.8%. At the upper limit of 48 GPa, the relative difference decreased to 2.7%. As the shear modulus is proportional to E, an increase in E increases the shear modulus. As a result, for a fixed stress value, the influence of the sliding mode increased when E decreased. As a result, in the mixed mode, the influence of the sliding mode on the interaction between fractures became more significant as E decreased.Figure 12
The effect of Young’s modulus on the relative difference.


### 4.3.3. Effects of Tectonic Strain Anisotropy

Moreover, we investigated the effects of tectonic strain anisotropy on fracture propagation since shale formations typically have a large strain anisotropy, which affects the direction of fracture propagation. Figure13 shows the minimum fracture-crossing angle obtained for various values of the strain ratio εH/εh. As the strain ratio increased, the minimum fracture-crossing angle decreased and fracture-crossing was more likely to occur. As shown in Figure 14, the relative difference increased as εH/εh increased. In the Barnett shale, when εH/εh was 3.00, the relative difference was 1.3%. The relative difference increased to 9.6% when εH/εh was 6.63. This trend existed because the rock became brittle through the increase in the shear stress and shear strain as εH/εh increased.Figure 13
The effect of the tectonic strain ratio on the minimum fracture-crossing angle.Figure 14
The effect of tectonic strain ratio on the relative difference.


## 4.3.1. Effects of Poisson’s Ratio

Figure9 depicts the minimum fracture-crossing angles obtained with various values of Poisson’s ratio ν in the range given in Table 4, where ν indicates the distribution of the rock deformation direction. The increase in ν increased the strain rate, thereby decreasing the minimum fracture-crossing angle. In other words, fracture-crossing is more likely to occur as ν increases.Figure 9
The effect of Poisson’s ratio on the minimum fracture-crossing angle.As shown in Figures9 and 10, all three shale formations showed a similar trend in terms of the minimum fracture-crossing angle and relative difference: the difference between the two calculated minimum fracture-crossing angles of the mixed mode and the opening mode divided by the minimum fracture-crossing angle with the mixed mode. However, their magnitudes were different. Figure 10 shows the relationship between ν and the relative difference. The relative difference increased as ν increased. For the Barnett shale, the relative difference was 1.7% when ν was at the lower limit of 0.05, and the relative difference increased to 13.2% when ν was at the upper limit of 0.26. This trend was also observed in the other two shale formations. These results were due to the interaction between the shear modulus and ν: because the shear modulus is inversely proportional to ν, the shear modulus decreases as ν increases. As a result, in the mixed mode, the influence of the sliding mode on the interaction between fractures became more significant as ν increased.Figure 10
The effect of Poisson’s ratio on the relative difference.


## 4.3.2. Effects of Young’s Modulus

We also analysed the fracture propagation behaviour considering the Young’s modulusE, which indicates the stiffness of the rock formation. Figure 11 shows the minimum fracture-crossing angles obtained with various values of E for the three shale formations. As E increased, the minimum fracture-crossing angle increased. As a result, fracture-crossing was difficult to occur.Figure 11
The effect of Young’s modulus on the minimum fracture-crossing angle.Figure12 shows the trend of the relative difference. For example, in the Barnett shale, when E was at the lower limit of 27 GPa, the relative difference was 9.8%. At the upper limit of 48 GPa, the relative difference decreased to 2.7%. As the shear modulus is proportional to E, an increase in E increases the shear modulus. As a result, for a fixed stress value, the influence of the sliding mode increased when E decreased. As a result, in the mixed mode, the influence of the sliding mode on the interaction between fractures became more significant as E decreased.Figure 12
The effect of Young’s modulus on the relative difference.


## 4.3.3. Effects of Tectonic Strain Anisotropy

Moreover, we investigated the effects of tectonic strain anisotropy on fracture propagation since shale formations typically have a large strain anisotropy, which affects the direction of fracture propagation. Figure13 shows the minimum fracture-crossing angle obtained for various values of the strain ratio εH/εh. As the strain ratio increased, the minimum fracture-crossing angle decreased and fracture-crossing was more likely to occur. As shown in Figure 14, the relative difference increased as εH/εh increased. In the Barnett shale, when εH/εh was 3.00, the relative difference was 1.3%. The relative difference increased to 9.6% when εH/εh was 6.63. This trend existed because the rock became brittle through the increase in the shear stress and shear strain as εH/εh increased.Figure 13
The effect of the tectonic strain ratio on the minimum fracture-crossing angle.Figure 14
The effect of tectonic strain ratio on the relative difference.


## 5. Conclusions

This study analysed the fracture propagation behaviours in shale gas reservoirs through simulations considering the effects of the mixed mode and the geomechanical properties of the reservoir mimicking three shale gas formations: Marcellus, Barnett, and Eagle Ford. The fracture propagation model was implemented with the mixed mode that linearly superimposed the opening mode and sliding mode and was combined with the uniaxial strain model. In our simulation results, the propagation results differed depending on the shale formation and fracture mode. We examined the effects of fracturing modes on the TCFV based on the BI value of each shale formation. In particular, the results indicate that neglecting the sliding mode might lead to an overestimation of the productivity in shale gas reservoirs with high BI values. In addition, the results of sensitivity analysis for various geomechanical properties revealed that the influence of the sliding mode increased when Poisson’s ratio increased, strain anisotropy increased, and Young’s modulus decreased. The sensitivity analysis implies that the proposed approach is extendable to other fields. In conclusion, our simulation results highlighted the significance of implementing the sliding mode in the mixed mode for the accurate modelling of fracture propagation in shale reservoirs with high brittleness.



---

*Source: 6690848-2020-12-23.xml*