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

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.

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]. Sev-eral experimental studies have investigated the interaction between hydraulic and natural fractures to comprehend fracture propagation [28][29][30][31][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 fracturecrossing 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-threedimensional 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][41][42][43][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][50][51][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.

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][58][59][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]. Figure 1 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.

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

Modelling of Hydraulic Fracture Propagation
Section 3 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.

Variables for the Modelling of Hydraulic Fracture
Propagation. Figure 3 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).

Geofluids
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 fracturecrossing or sliding occurs at the intersection. If β is greater than the minimum fracture-crossing angle, fracturecrossing 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.

Workflow for the Modelling of Hydraulic Fracture
Propagation. The variables presented in Figure 3 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 a max 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 fracturecrossing 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 a max . 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 fluidflow simulator is coupled with the DFN model through grid mapping. In this study, all source codes were written in the C programming language.
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]: 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  4 Geofluids 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 xand 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): where r is the radial distance from the fracture tip, θ is the polar angle to the direction of the fracture tip, and K I and K II 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. K I and K II are computed as follows [74]: where P is the fracturing pressure and a is the fracture halflength.

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: 3.2.3. Fracture Crossing Criterion with Mixed Mode. If the hydraulic fracture encounters a natural fracture nonortho-gonally, 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: where τ β 1 is the combined shear stress and σ β 2 is the combined normal stress, as noted in Figure 3; S 0 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 S 0 − μσ β 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: 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 K 1 is the stress level (i.e., substitution variable) required to reinitiate a fracture. Note that K 1 is the larger solution of the quadratic polynomial equation below: where A = cos 2 β 2 − 3 tan 2 α sin 2 β 2 + 2 tan α tan β 2 cos β + 1, ð11Þ where T 0 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].

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: If fracture crossing does not occur and the sum of the lengths of all connected fractures is smaller than a max , the fracture is reinitiated at the end of the natural fracture using Equation (6).

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 fineto very fine-grained sandstones, muddy sandstones, and sandy mudstones. The matrix porosity and matrix permeability are less than 8% and 1.0E-3 μm 2 , respectively. Table 1 presents the vertical stress σ v , minimum horizontal stress σ h , pore pres-sure 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 (R 2 ) was 0.901, and the root-mean-square error (RMSE) was 1.369.
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 T 0 was 4.05 MPa, and the cohesion S 0 was zero. The model results agreed with all the experimental results.    [50] was used to quantitatively compare fracture propagation behaviours in shale gas formations in Section 4. The BI is defined as follows: where E min and E max 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 E D into the static Young's modulus E using Equation (16) [77] as follows:

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   7 Geofluids analysis to investigate the effects of the geomechanical properties of the reservoir on the contribution of the sliding mode to the mixed mode. Figure 6 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, a max 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 E D , 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, E D 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 T 0 . Here, the values of T 0 were constants of 3, 4, and 6 MPa for the Marcellus, Barnett, and Eagle Ford formations, respectively.

Experimental Setting.
Although this study utilises field data acquired from references (Table 4), 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 a max is given to each case.

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 Table 4. The inclination angle α of 15°yielded a fracture initiation angle θ i of 7.5°. The simulation results with the mixed mode were also 8 Geofluids compared with the simulation results with the opening mode to highlight the effects of the mixed mode on fracture propagation. Figure 7 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.
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 m 3 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 m 3 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.

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

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 fracturecrossing angle, ν, E, and ε H /ε h were selected as sensitivity analysis factors. Figure 9 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.

Effects of Poisson's Ratio.
As shown in Figures 9 and 10, all three shale formations showed a similar trend in terms of the minimum fracturecrossing 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.

Effects of Young's Modulus.
We also analysed the fracture propagation behaviour considering the Young's modulus E, 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 12 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.

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. Figure 13 shows the minimum fracturecrossing angle obtained for various values of the strain ratio ε H /ε h . As the strain ratio increased, the minimum fracturecrossing 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.

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.

Nomenclature a:
Fracture half-length, m a max : Maximum fracture half-length, m E: Young's modulus, GPa E D : Dynamic Young's modulus, GPa Relative difference (%) Figure 14: The effect of tectonic strain ratio on the relative difference.
11 Geofluids E max : Maximum Young's modulus in the interval of interest, GPa E min : Minimum Young's modulus in the interval of interest, GPa K I : Stress intensity factor of opening-mode fractures, MPa·m 1/2 K II : Stress intensity factor of sliding-mode fractures, MPa·m 1/2 K 1 : Stress level of opening-mode P: Hydraulic fracturing pressure, MPa p: Pore pressure, MPa r: Radial distance from the fracture tip, m S 0 : Cohesion of the natural fracture interface, MPa T 0 : Tensile strength of porous media, MPa α: Inclination angle between the direction of hydraulic fracturing and the direction of maximum horizontal stress, degrees β: Intersection angle between hydraulic fracture and natural fracture, degrees γ: Fracture reinitiation angle, degrees δ 1 : Biot's coefficient ε H : Maximum tectonic compressive strain ε h : Minimum tectonic compressive strain θ: Polar angle to the direction of fracture tip, degrees θ i : Fracture initiation angle, degrees μ: Friction coefficient ν: Poisson's ratio ν max : Maximum Poisson's ratio in the interval of interest ν min : Minimum Poisson's ratio in the interval of interest σ H : Maximum horizontal stress, MPa σ h : Minimum horizontal stress, MPa σ r,β 2 : Normal stress on natural fracture interface, MPa σ tip,β 2 : Normal stress at hydraulic fracture tip, MPa σ v : Vertical stress, MPa σ x : Normal stress to the x-direction, MPa σ y : Normal stress to the y-direction, MPa σ β 2 : Combined normal stress, MPa τ r,β 1 : Shear stress on natural fracture interface, MPa τ tip,β 1 : Shear stress at hydraulic fracture tip, MPa τ xy : Shear stress on xy-plane, MPa τ β 1 : Combined shear stress, MPa.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.