# Numerical Investigation of Hydraulic Fracture Propagation in Naturally Fractured Reservoirs Based on Lattice Spring Model

Authors: Kaikai Zhao; Pengfei Jiang; Yanjun Feng; Xiaodong Sun; Lixing Cheng; Jianwei Zheng
Journal: Geofluids
 (2020)

---

## Abstract

Hydraulic fracturing has been extensively employed for permeability enhancement in low-permeability reservoirs. The geometry of the hydraulic fracture network (HFN) may have implications for the optimization of hydraulic fracturing operations. Various parameters, including the in situ stress, treatment parameters (injection rate and fluid viscosity), and orientation of natural fractures (NFs), can significantly affect the interactions between hydraulic fracture (HF) and NFs and the final HFN. In this study, a lattice-spring code was employed to determine the impact of various parameters on the geometry of the HFN. The modelling results indicated that with a large stress difference, the global orientation of the fracture propagation was restricted to the direction of maximum principal stress, and the number of fracture branches was reduced. The geometry of the HFN changed from circular to elliptical. In contrast, with an increase in the fluid viscosity/injection rate, the evolution of the geometry of the HFN exhibited the opposite trend. The global orientation of HF propagation tended to remain parallel to the direction of maximum principal stress, regardless of the branching and tortuosity of the fracture. The variations in the ratio of tensile fracture (HF) to shear fracture (shear slip on NF) can be significant, depending on the stress state, treatment parameters, and preexisting NF network, which determine the dominant stimulation mechanism. This study provides insight into the HF propagation in naturally fractured reservoirs.


---

## Body

## 1. Introduction

Hydraulic fracturing treatment has been widely applied in the shale gas reservoir [1] as well as coal seam gas reservoir [2]. It is increasingly being used for preconditioning of the orebody in cave mining [3]. In the naturally fractured formation, a hydraulic fracture (HF) may encounter natural fractures (NFs) of various scales, such as joints, bedding planes, and faults. Several types of interactions (e.g., diversion, offsetting, and crossing) can occur when the HF encounters the NF. Therefore, hydraulic fracturing treatment in a naturally fractured reservoir may give rise to a complex hydraulic fracture network (HFN) instead of a symmetric, planar, biwing HF [4]. By predicting the HFN geometry, the accuracy of the hydrofracture simulation for fractured reservoirs can be improved [5].Various parameters, including the in situ stress, treatment parameters (injection rate and fluid viscosity), and geometric and mechanical properties of the NF, can significantly affect the HF–NF interactions and the final HFN. Numerous experiments have been performed to investigate the interactions between HFs and NFs. Zhou et al. [6, 7] argued that the stress difference, shear strength of the NF, and approach angle (intersection angle between the HF and NF) are crucial factors determining the HF propagation behavior in the fractured formation. HFs tend to cross preexisting NFs under a large stress difference and approach angle, whereas they undergo diversion/deflection due to the NF under a small stress difference and approach angle. In laboratory experiments performed by Beugelsdijk et al. [8] With a small value of the product of the injection rate Q and the fluid viscosity μ (Qμ), fluid tended to leak into the NFs, resulting in tortuous HF propagation paths following the NFs. With a high Qμ value, the HF tended to cross most NFs, and the overall propagation path was relatively straight. Zou et al. [9] performed a series of experiments to investigate HF propagation using computed tomography scanning technology. The results indicated that the NF network (NFN) was activated for a small horizontal stress difference of <6 MPa, and a simple transverse fracture pattern was observed for a large horizontal stress difference of >9 MPa. Additionally, a dominant HF was observed for treatment with a high injection rate, whereas the NFN was activated to a large extent under a low injection rate. In the field, the HFN has been characterized by combining microseismic analysis with surface and downhole tilt fracture mapping [10]. As shown in Figure 1, field observations revealed varying degrees of complexity, ranging from a simple, relatively planar fracture to a complex fracture network. Mayerhofer et al. [11] proposed the concept of the stimulated reservoir volume (three-dimensional (3D) volume of the microseismic event cloud) as a correlation parameter for good performance. The stimulated reservoir volume can approximate the size of the created HFN. Because of the limitations regarding the size of the rock sample and the precision of the measurement device, it is challenging to perform a sensitivity analysis or quantitively evaluate the effects of various parameters on the interactions between an HF and multiple NFs, as well as the final HFN.Figure 1
Schematics of the levels of complexity observed in HFs [4].Comprehensive numerical models have been proposed to investigate the HF propagation in naturally fractured formations, which can be categorized according to their numerical methods: the finite element method (FEM), including the extended finite element method (XFEM) [12, 13] and cohesive zone method [14, 15], boundary element method (BEM) [16], displacement discontinuity method (DDM) [17, 18], distinct element method (DEM) [19, 20], and lattice method [21]. The simulation methods for hydraulic fracturing have recently been comprehensively reviewed [22–24]. Taleghani and Olson [25] presented an XFEM model considering the interactions between HFs and NFs. The modelling results indicated that the fracture-pattern complexity is significantly affected by the stress anisotropy, rock toughness, and NF strength, as well as the orientation of the NF. Abbas et al. [26] adopted an XFEM model to examine the effects of different combinations of parameters (i.e., formation moduli, far-field stresses, and injection rates) on the HF height and the size of the HF opening. Ghaderi et al. [27] used the XFEM method to simulate the deformation of the NF during the approaching stage of the HF. The results indicated that the tensile and shear debonding of the NF change with respect to the angle and distance from an NF. Zhang et al. [16] investigated the HF deflection behaviors at bedding interfaces using a two-dimensional BEM model. The HF deflection and fluid invasion depend on various parameters, including the elastic-modulus contrasts, in situ stresses, interfacial frictional coefficients, and fluid viscosities. Olson [28] presented a complex fracture network model to simulate HF propagation and the interaction between an HF and an NF using a pseudo-3D DDM. Kresse and Weng [18] developed an unconventional fracture model for simulating HF propagation, rock deformation, and fluid flow in a complex fracture network.The DEM has been widely adopted in various rock engineering projects [29, 30]. The typical DEM can effectively reproduce the open/slip of an NF and the interaction between blocks and an HF. However, no new HF propagation beyond the prebuilt trajectories can be reproduced. The synthetic rock mass (SRM) approach compensates for the deficiency of the prebuilt trajectory in conventional DEM models [31]. The SRM scheme has been incorporated in the lattice scheme code XSite [32]. Bakhshi et al. [33] adopted XSite to investigate the intersection of an HF with an NF with consideration of the effects of the intersection angle and the mechanical properties of the NF. Zhao et al. [34] employed XSite to simulate the 3D interaction between an HF and an NF, with consideration of the effects of the stress difference, treatment parameters, and NF properties. Liu et al. [35] adopted XSite to study the stress interference between multiple HFs in a horizontal well. Wan et al. [36] used XSite to investigate the effects of the rock properties and in situ stresses on HF containment. In most of the foregoing studies, the wellbore/open-hole was treated as an injection point or a predefined fracture path, and the effects of stress concentration around the wellbore (and potentially on the fracture tortuosity and branching near the wellbore) were neglected.This study focused on the numerical modelling of the HF propagation in a naturally fractured formation and the evolution of the HFN geometry. A series of XSite simulation was conducted to investigate the effects of the stress difference (∆σ=σ1–σ3), the treatment parameters (fluid viscosity and injection rate), and the orientation of the NF. Additionally, the HF propagation in several typical NFNs was analyzed.


## 2. Modelling Methodology

The code XSite, based on the lattice method (see Figure2) and SRM approach (see Figure 3), is a more computationally efficient version of the 3D particle flow code. In SRM, the bonded particle model (BPM) was employed to represent intact material, and the smooth joint model (SJM) was used to describe the joints behaviors [32]. Lattice simulation is a simplification of the BPM in which the particles and contacts are replaced by nodes and springs.Figure 2
Schematic diagram of lattice simulation method.Figure 3
Schematic diagram of the synthetic rock mass (SRM) approach [21].There are two methods used to generate the springs that connect the nodes: regular and Voronoi. Voronoi lattice tessellation is utilized in the presented simulations, with the springs placed based on Voronoi tessellation in 3D space, where the springs are created at common faces of the discretization domains. The lattice is created by multiplication of the periodic brick (p-brick) in three orthogonal directions. The p-brick is a quasirandom arrangement of nodes within a cube of unit edge length. The final model geometry is achieved by trimming of the “excess” lattice extending outside the analyzed domain [32].


### 2.1. Mechanical Formulation

As shown in Figure4, the lattice is composed of numerous quasirandomly distributed nodes connected by springs. Joints are overlaid on the lattice using the SJM methodology.Figure 4
Schematic of a lattice array [37].The central difference method is employed to compute the transitional degrees of freedom [32]:
(1)u̇it+Δt/2=u̇it−Δt/2+∑FitΔtm,(2)uit+Δt=uit+u̇it+Δt/2Δt,where u̇itand uitrepresent the velocity and position of component i (i=1,3) at time t. ∑Fit represent the sum of all the force components i acting on a node of mass m, with time step Δt.The angular velocities of componenti at time t are calculated using
(3)ωit+Δt/2=ωit−Δt/2+∑MitIΔt,where ∑Mit is the sum of all moment-components acting on the node of moment of inertia I.The force change in the spring is determined by the displacements of the node [32]:
(4)FN←FN+u̇NkNΔt,(5)FiS←FiS+u̇iSkSΔt,where N represents “normal,” S represents “shear,” F represents the spring force. kN and kS represent the spring normal and shear stiffness, respectively. If the force exceeds the spring strength, the spring breaks and a microcrack is formed.Joint slip and opening follow the relationship [37]
(6)IfFn−pA<0thenFn=0,Fis=0,elseFis⇐Fis∣Fis∣minFn−pAtanϕ,∣Fis∣,where Fn represents the normal force, FiS represents the shear force vector, p represents the pressure, A represents the apparent area, and ϕ represents the friction angle.The following relationship determines the bonded joint status: ifFn−pA+σcA<0 or Fis>τcA, the bond fails in tension or shear (where σc is the bond tensile strength, τc is the bond shear strength), else,Fis←Fis, the bond remains intact.


### 2.2. Flow Formulation

As shown in Figure5, the flow in HF is simulated using fluid elements linked by pipes. The fluid elements act as microcracks which are positioned at the centers of the broken springs or springs overlapped by the joint.Figure 5
Schematic of the pipe network [32].Lubrication theory is used to calculate the flow rate from node A to node B along a pipe:(7)q=βkra312μfpA−pB+ρwgzA−zB,where a represents the aperture; μf represents the fluid viscosity; pA and pB represents the hydraulic pressures at nodes “A” and “B,” respectively; zA and zB represents the elevations of nodes “A” and “B,” respectively; ρw represents the fluid density; g represents the gravitational acceleration; and β is a calibration parameter that reflects the conductivity.The relative permeability,kr, is a function of saturation, s:
(8)kr=s23−2s,The pressure increment,ΔP, during the flow timestep, Δtf, is calculated as
(9)ΔP=QVKfΔtf,where Q is the sum of flow rate from the pipes connected to the fluid element, V is the volume of the fluid element, and Kf is the apparent fluid element bulk modulus.


### 2.3. Hydromechanical Coupling

The fluid flow and mechanical process are fully coupled (see Figure6). The mechanical deformation and damage are computed based on the variation in fluid pressure. In contrast, the variation in fluid pressure depends on the mechanical deformation. The HF permeability is decided by the HF aperture and mechanical deformation.Figure 6
Couplings between mechanical and hydraulic processes.


### 2.4. Fracture Propagation Criteria

The criteria for HF propagation are based on anJ-integral formulation. The stress intensity factor, KI, can be calculated as
(10)KI=JE,where E represents the Young’s modulus. If KI<KIC (where KIC is the rock toughness); then, the spring tensile strength is utilized to detect spring failure. Otherwise, KI is compared to KIC to detect spring failure.


### 2.5. Model Setup

Figure7 shows the typical model setup. The dimensions of the rock mass were 5m×5m×0.3m. NFs with dimensions of 0.25m×0.25 m were evenly distributed in the rock mass. The initial aperture of the NFs was 0.1 mm. A wellbore was placed at the center of the model, parallel to the y-axis. The wellbore had a radius of 0.075 m and a length of 0.3 m. The horizontal stress was set as σx=σy=5MPa for all the models. The mechanical and hydraulic parameters for the typical model are presented in Table 1.Figure 7
Setup of the base model.Table 1
Mechanical and hydraulic model input parameters.
CategoriesVariablesValuesIn situ stressStress difference,∆σ (MPa)0, 2.5, 5, 7.5, 10Treatment parametersInjection rate,Q (m3/s)0.002, 0.003, 0.004, 0.005, 0.007Fluid viscosity,μ (mPa·s)1, 2, 3, 4, 5Natural fractureDip angle,θ (°)0, 15, 30, 45, 60°Intact rockYoung’s modulus,E (GPa)11.74Poisson’s ratio,ν0.221Tensile strength,σt (MPa)7.5UCS,σc (MPa)75


## 2.1. Mechanical Formulation

As shown in Figure4, the lattice is composed of numerous quasirandomly distributed nodes connected by springs. Joints are overlaid on the lattice using the SJM methodology.Figure 4
Schematic of a lattice array [37].The central difference method is employed to compute the transitional degrees of freedom [32]:
(1)u̇it+Δt/2=u̇it−Δt/2+∑FitΔtm,(2)uit+Δt=uit+u̇it+Δt/2Δt,where u̇itand uitrepresent the velocity and position of component i (i=1,3) at time t. ∑Fit represent the sum of all the force components i acting on a node of mass m, with time step Δt.The angular velocities of componenti at time t are calculated using
(3)ωit+Δt/2=ωit−Δt/2+∑MitIΔt,where ∑Mit is the sum of all moment-components acting on the node of moment of inertia I.The force change in the spring is determined by the displacements of the node [32]:
(4)FN←FN+u̇NkNΔt,(5)FiS←FiS+u̇iSkSΔt,where N represents “normal,” S represents “shear,” F represents the spring force. kN and kS represent the spring normal and shear stiffness, respectively. If the force exceeds the spring strength, the spring breaks and a microcrack is formed.Joint slip and opening follow the relationship [37]
(6)IfFn−pA<0thenFn=0,Fis=0,elseFis⇐Fis∣Fis∣minFn−pAtanϕ,∣Fis∣,where Fn represents the normal force, FiS represents the shear force vector, p represents the pressure, A represents the apparent area, and ϕ represents the friction angle.The following relationship determines the bonded joint status: ifFn−pA+σcA<0 or Fis>τcA, the bond fails in tension or shear (where σc is the bond tensile strength, τc is the bond shear strength), else,Fis←Fis, the bond remains intact.


## 2.2. Flow Formulation

As shown in Figure5, the flow in HF is simulated using fluid elements linked by pipes. The fluid elements act as microcracks which are positioned at the centers of the broken springs or springs overlapped by the joint.Figure 5
Schematic of the pipe network [32].Lubrication theory is used to calculate the flow rate from node A to node B along a pipe:(7)q=βkra312μfpA−pB+ρwgzA−zB,where a represents the aperture; μf represents the fluid viscosity; pA and pB represents the hydraulic pressures at nodes “A” and “B,” respectively; zA and zB represents the elevations of nodes “A” and “B,” respectively; ρw represents the fluid density; g represents the gravitational acceleration; and β is a calibration parameter that reflects the conductivity.The relative permeability,kr, is a function of saturation, s:
(8)kr=s23−2s,The pressure increment,ΔP, during the flow timestep, Δtf, is calculated as
(9)ΔP=QVKfΔtf,where Q is the sum of flow rate from the pipes connected to the fluid element, V is the volume of the fluid element, and Kf is the apparent fluid element bulk modulus.


## 2.3. Hydromechanical Coupling

The fluid flow and mechanical process are fully coupled (see Figure6). The mechanical deformation and damage are computed based on the variation in fluid pressure. In contrast, the variation in fluid pressure depends on the mechanical deformation. The HF permeability is decided by the HF aperture and mechanical deformation.Figure 6
Couplings between mechanical and hydraulic processes.


## 2.4. Fracture Propagation Criteria

The criteria for HF propagation are based on anJ-integral formulation. The stress intensity factor, KI, can be calculated as
(10)KI=JE,where E represents the Young’s modulus. If KI<KIC (where KIC is the rock toughness); then, the spring tensile strength is utilized to detect spring failure. Otherwise, KI is compared to KIC to detect spring failure.


## 2.5. Model Setup

Figure7 shows the typical model setup. The dimensions of the rock mass were 5m×5m×0.3m. NFs with dimensions of 0.25m×0.25 m were evenly distributed in the rock mass. The initial aperture of the NFs was 0.1 mm. A wellbore was placed at the center of the model, parallel to the y-axis. The wellbore had a radius of 0.075 m and a length of 0.3 m. The horizontal stress was set as σx=σy=5MPa for all the models. The mechanical and hydraulic parameters for the typical model are presented in Table 1.Figure 7
Setup of the base model.Table 1
Mechanical and hydraulic model input parameters.
CategoriesVariablesValuesIn situ stressStress difference,∆σ (MPa)0, 2.5, 5, 7.5, 10Treatment parametersInjection rate,Q (m3/s)0.002, 0.003, 0.004, 0.005, 0.007Fluid viscosity,μ (mPa·s)1, 2, 3, 4, 5Natural fractureDip angle,θ (°)0, 15, 30, 45, 60°Intact rockYoung’s modulus,E (GPa)11.74Poisson’s ratio,ν0.221Tensile strength,σt (MPa)7.5UCS,σc (MPa)75


## 3. Results

### 3.1. Effect of Stress Difference

For the five models, the vertical stress,σz, was set as 5, 7.5, 10, 12.5, and 15 MPa, and σx=σy=5MPa. The corresponding stress differences were 0, 2.5, 5, 7.5, and 10 MPa, respectively. Simulations were conducted under a constant injection rate of 0.005 m3/s and a fluid viscosity of 1 mPa·s for 1 s.As shown in Figure8, multiple radial branches from the wellbore were simulated, and the tortuosity of the propagation pathway due to the interactions between the NF and the HF was determined. The NF locally altered the HF propagation pathway through the diversion of the HF or induced branching of the HF. As the stress difference changed from 0 to 10 MPa, the dominant propagation direction of the HFN became more evident, tending to follow the direction of maximum principal stress. The geometry of the stimulated region (represented by red dotted lines) changed from circular to elliptical, and the long axes tended to be along the direction of maximum principal stress. More branches were observed in small-stress-difference cases than in large-stress-difference cases. The branches were inhibited under a large stress difference, and the dominant propagation direction was restricted to the direction of maximum principal stress.Figure 8
HF propagation in the NFN with varying stress differences: (a) 0 MPa, (b) 2.5 MPa, (c) 5 MPa, (d) 7.5 MPa, and (e) 10 MPa.Figure9 shows the variations in the maximum vertical growth (Lz) and maximum lateral growth (Lx) of the HFN for the five models. With an increase in the stress difference, Lx decreased from 4.62 to 2.49 m (by 2.13 m). Conversely, Lz increased with an increase in the stress difference (except for the discrete points, the case of ∆σ=7.5MPa was examined). For the case of ∆σ=0MPa, Lx was greater than Lz. When the stress difference exceeded zero, Lz was greater than Lx. Moreover, the difference between Lz and Lx increased as the stress difference increased from 2.5 to 10 MPa.Figure 9
Variations of the maximum vertical growth (Lz) and lateral growth (Lx) under varying stress differences.Figure10 shows the variation of the maximum deviation angle α (maximum angle between the branching and the direction of σ1) for the five models. In general, α decreased with the increasing stress difference. α decreased significantly (by 31°) as the stress difference increased from 0 to 5 MPa. α decreased slightly (by 13°) as the stress difference increased from 5 to 10 MPa.Figure 10
Variation of the maximum deviation angle (α) for varying stress differences.


### 3.2. Effect of Fluid Injection Rate

Five models were used for the simulations, with assumed injection rates of 0.002, 0.003, 0.004, 0.005, and 0.007 m3/s. The volume of the injection fluid was 0.005 m3 for all the models. The other parameters were as follows: σx=σy=5MPa, σz=7.5MPa, and μ=1 mPa·s.As shown in Figure11, in general, the geometry of the stimulated region changed from elliptical to circular with an increase in the injection rate. Additionally, there were fewer primary branches at lower injection rates (see Figures 8(a) and 8(b)). Four primary branches were observed in the lowest injection rate case (Q=0.002 m3/s), and each branch was longer than the corresponding branch in the highest injection rate case (Q=0.007 m3/s). The HF mainly propagated along the direction of σ1 for low injection rates. However, no clear dominant propagation direction was observed at the highest injection rate (Q=0.007 m3/s).Figure 11
HF propagation in the NFN with injection rates of (a) 0.002 m3/s, (b) 0.003 m3/s, (c) 0.004 m3/s, (d) 0.005 m3/s, and (e) 0.007 m3/s.As shown in Figure12, the Lz was greater than the Lx for all five models, and among the models, the difference between the Lz and the Lx was the smallest for Q=0.007 m3/s.Figure 12
Variations of the maximum vertical and lateral growth for different injection rates.As shown in Figure13, α increased with the injection rate. There was no significant variation in α as the injection rate increased from 0.002 to 0.004 m3/s. However, there was a significant increase (53°) in α as the injection rate increased from 0.004 to 0.007 m3/s.Figure 13
Variation of the maximum deviation angle for different injection rates.


### 3.3. Effect of Fluid Viscosity

Five models were used for the simulations, with hydraulic fracturing fluid viscosities of 1, 2, 3, 4, and 5 mPa·s. The stress state was set as follows:σx=σy=5MPa, σz=7.5MPa. The simulations were conducted under a constant injection rate of 0.002 m3/s for 2.5 s.As shown in Figure14, the geometry of the stimulated region changed from elliptical to circular as the fluid viscosity increased from 1 to 5 mPa·s. The size of the stimulated region decreased with the increasing fluid viscosity. Additionally, more branches were induced under treatment with higher-viscosity fluids. The HF mainly propagated along the direction of σ1 for cases with a relatively low-viscosity fluid. However, no dominant propagation direction was observed for the high-viscosity cases of μ=4 and 5 mPa·s.Figure 14
HF propagation in the NFN with different fluid viscosities: (a) 1 mPa·s, (b) 2 mPa·s, (c) 3 mPa·s, (d) 4 mPa·s, and (e) 5 mPa·s.As shown in Figure15, Lz was greater than Lx (except for the case of μ=4 mPa·s), and the difference between Lz and Lx for the high-fluid viscosity cases of μ=3, 4, and 5 mPa·s was smaller than that for the lower-fluid viscosity cases of μ=1 and 2 mPa·s.Figure 15
Variations of the maximum vertical and lateral growth for different fluid viscosities.As shown in Figure16, α tended to increase with the fluid viscosity. A significant increase (53°) in α was observed when the fluid viscosity increased from 2 to 4 mPa·s. The deviation angle was maximized (approximately 90°) in the high-fluid viscosity cases of μ=4 and 5 mPa·s.Figure 16
Variation of the maximum deviation angle for different fluid viscosities.


### 3.4. Effect of NF Orientation

Five models were used for the simulations, with dip angles of 0°, 15°, 30°, 45°, and 60°. The angles between the NFs and the direction ofσ1 were 90°, 75°, 60°, 45°, and 30°, respectively. The stress state was as follows: σx=σy=5 MPa, σz=5.5 MPa. Simulations were performed under a constant injection rate of 0.0005 m3/s for 2 s.As shown in Figure17, the HFN geometry varied significantly with changes in the NF dip angle. The direction of the long axis of the stimulated region changed from subhorizontal to subvertical as the NF dip angle increased from 0° to 60°. As shown in Figure 14(a), the two HF branches were initiated at the wellbore and then encountered the two NFs closest to the wellbore. Fluid invasion into the two NFs occurred, causing the open and shear slip of the NFs and then extending from the edges of the NFs. This process occurred again when the HF encountered the next NF. Thus, a step-like HFN was formed in the case of θ=0°. With an increase in the dip angle (and a corresponding reduction in the angle between the NF and the direction of σ1), the number of primary branches decreased, and the main propagation direction of the HF network became closer to the direction of σ1. There were four primary branches in the case of θ=0°, whereas only two branches were observed for θ=60°.Figure 17
HF propagation in the NFN with different NF dip angles: (a) 0°, (b) 15°, (c) 30°, (d) 45°, and (e) 60°.As shown in Figure18, as the dip angle increased from 0° to 60°, Lz decreased from 3.06 to 0.84 m (by 2.22 m). Conversely, Lx increased from 1.66 to 4.72 m. In the case of θ=0°, Lx was greater than Lz, whereas Lz exceeded Lx for θ>15°. The maximum difference between Lz and Lx occurred for θ=60°.Figure 18
Variation of the maximum deviation angle for different NF dip angles.Taking caseθ=45° as an example for detailed investigation, Figure 19 shows the microcracks (tensile failure of intact rock) and slip event (shear failure in NF). The HFN comprised the tensile failure in the intact rock and the shear failure in the connected NFs. The connected NFs appeared to be fully activated by shear slippage. Additionally, isolated NFs were subjected to shear failure, and slip events were observed on isolated NFs.Figure 19
Microcracks in intact rock and shear slip on NFs forθ=45°.The displacement field is shown in Figure20. Three primary branches divided the intact block into three separate blocks, which were characterized by different displacement fields. In general, the displacement decreased with the increasing distance from the wellbore. Most of the block experienced displacement of >0.2 mm. The largest displacement was observed for block 3.Figure 20
Discontinuous displacement field forθ=45°.


### 3.5. HF Propagation in Different NFNs

We investigated the HF propagation behavior in different NFNs and considered three types of simplified NFNs. The type A NFN comprised several large, parallel-distributed NFs (Figure21). In the type B NFN, there were vertical NFs between adjacent horizontal NFs, and all the NFs were isolated from each other (Figure 22). In type C NFN, the vertical NFs were connected with the horizontal NFs (Figure 23). The rock mass volume was 9m×9m×0.6 m. A starter crack with a radius of 0.15 m was located at the center of the model, normal to the x-axis. The stress state was as follows:σz=13MPa, σx=σy=5. Simulations were conducted under the injection rate of Q=0.003 m3/s for 12 s.Figure 21
Type A NFN: (a) model setup, (b) fracture aperture, (c) fracture fluid pressure, and (d) shear slip on the NF.
(a)(b)(c)(d)Figure 22
Type B NFN: (a) model setup, (b) fracture aperture, (c) fracture fluid pressure, and (d) shear slip on the NF.
(a)(b)(c)(d)Figure 23
Type C NFN: (a) model setup, (b) fracture aperture, (c) fracture fluid pressure, and (d) shear slip on the NF.
(a)(b)(c)(d)As shown in Figure21(b), a vertical HF was induced and crossed all the horizontal NFs, yielding a fishbone-like HFN. In general, the aperture of the horizontal NF decreased with the increasing distance from the starter crack (the injection point). For instance, the two horizontal NFs closest to the injection point had greater apertures (>1 mm) than the other horizontal NFs. Note that the aperture is varied in along the horizontal NF plane. Figure 21(c) shows the approximately radial distribution of the fracture fluid pressure, which decreased with increasing distance from the injection point. Figure 21(d) shows the shear slip on the horizontal NFs. Among the NFs, the largest shear slip was observed on the two horizontal NFs closest to the injection point.As shown in Figure22(b), the vertical flow pathway contained the newly induced HF and several vertical NFs. The flow network did not connect most of the vertical NFs. No significant difference in the fracture aperture was observed in comparison with type A. The approximately radial distribution pattern of the fracture fluid pressure is presented in Figure 22(c). The fluid pressure varied slightly compared with type A. The induced shear slip on the horizontal NFs was similar between types A and B. Additionally, the vertical NFs connected by the HF appeared to be fully activated by shear failure. Most of the other vertical NFs also underwent a slight slip. In general, the HF propagation was similar between types A and B.The HF propagation for type C differed significantly from the two aforementioned cases. As shown in Figure23(b), the fluid was forced to pass through the connected NFN, and only a few new HFs were induced. The HF extended from the tip of the vertical NF, which was located on the edge of the NFN. The dominant flow path was still along the preferential fracture plane (the direction of maximum principal stress). As shown in Figure 23(c), the distribution pattern of the fracture fluid pressure was approximately radial. However, the fluid pressure gradient was reduced in comparison with those for types A and B, owing to the high connectivity of the NFN. Specifically, the size of the region with high fluid pressure (>20 MPa) was smaller than those for types A and B. The fluid pressure on the edge of the network exceeded approximately 8 MPa. The number of vertical NFs that experienced shear slip was significantly larger than that for type B (Figure 23(d)). A large shear zone along the vertical direction was simulated, in which the vertical NFs appeared to be fully activated by shear slippage.


## 3.1. Effect of Stress Difference

For the five models, the vertical stress,σz, was set as 5, 7.5, 10, 12.5, and 15 MPa, and σx=σy=5MPa. The corresponding stress differences were 0, 2.5, 5, 7.5, and 10 MPa, respectively. Simulations were conducted under a constant injection rate of 0.005 m3/s and a fluid viscosity of 1 mPa·s for 1 s.As shown in Figure8, multiple radial branches from the wellbore were simulated, and the tortuosity of the propagation pathway due to the interactions between the NF and the HF was determined. The NF locally altered the HF propagation pathway through the diversion of the HF or induced branching of the HF. As the stress difference changed from 0 to 10 MPa, the dominant propagation direction of the HFN became more evident, tending to follow the direction of maximum principal stress. The geometry of the stimulated region (represented by red dotted lines) changed from circular to elliptical, and the long axes tended to be along the direction of maximum principal stress. More branches were observed in small-stress-difference cases than in large-stress-difference cases. The branches were inhibited under a large stress difference, and the dominant propagation direction was restricted to the direction of maximum principal stress.Figure 8
HF propagation in the NFN with varying stress differences: (a) 0 MPa, (b) 2.5 MPa, (c) 5 MPa, (d) 7.5 MPa, and (e) 10 MPa.Figure9 shows the variations in the maximum vertical growth (Lz) and maximum lateral growth (Lx) of the HFN for the five models. With an increase in the stress difference, Lx decreased from 4.62 to 2.49 m (by 2.13 m). Conversely, Lz increased with an increase in the stress difference (except for the discrete points, the case of ∆σ=7.5MPa was examined). For the case of ∆σ=0MPa, Lx was greater than Lz. When the stress difference exceeded zero, Lz was greater than Lx. Moreover, the difference between Lz and Lx increased as the stress difference increased from 2.5 to 10 MPa.Figure 9
Variations of the maximum vertical growth (Lz) and lateral growth (Lx) under varying stress differences.Figure10 shows the variation of the maximum deviation angle α (maximum angle between the branching and the direction of σ1) for the five models. In general, α decreased with the increasing stress difference. α decreased significantly (by 31°) as the stress difference increased from 0 to 5 MPa. α decreased slightly (by 13°) as the stress difference increased from 5 to 10 MPa.Figure 10
Variation of the maximum deviation angle (α) for varying stress differences.


## 3.2. Effect of Fluid Injection Rate

Five models were used for the simulations, with assumed injection rates of 0.002, 0.003, 0.004, 0.005, and 0.007 m3/s. The volume of the injection fluid was 0.005 m3 for all the models. The other parameters were as follows: σx=σy=5MPa, σz=7.5MPa, and μ=1 mPa·s.As shown in Figure11, in general, the geometry of the stimulated region changed from elliptical to circular with an increase in the injection rate. Additionally, there were fewer primary branches at lower injection rates (see Figures 8(a) and 8(b)). Four primary branches were observed in the lowest injection rate case (Q=0.002 m3/s), and each branch was longer than the corresponding branch in the highest injection rate case (Q=0.007 m3/s). The HF mainly propagated along the direction of σ1 for low injection rates. However, no clear dominant propagation direction was observed at the highest injection rate (Q=0.007 m3/s).Figure 11
HF propagation in the NFN with injection rates of (a) 0.002 m3/s, (b) 0.003 m3/s, (c) 0.004 m3/s, (d) 0.005 m3/s, and (e) 0.007 m3/s.As shown in Figure12, the Lz was greater than the Lx for all five models, and among the models, the difference between the Lz and the Lx was the smallest for Q=0.007 m3/s.Figure 12
Variations of the maximum vertical and lateral growth for different injection rates.As shown in Figure13, α increased with the injection rate. There was no significant variation in α as the injection rate increased from 0.002 to 0.004 m3/s. However, there was a significant increase (53°) in α as the injection rate increased from 0.004 to 0.007 m3/s.Figure 13
Variation of the maximum deviation angle for different injection rates.


## 3.3. Effect of Fluid Viscosity

Five models were used for the simulations, with hydraulic fracturing fluid viscosities of 1, 2, 3, 4, and 5 mPa·s. The stress state was set as follows:σx=σy=5MPa, σz=7.5MPa. The simulations were conducted under a constant injection rate of 0.002 m3/s for 2.5 s.As shown in Figure14, the geometry of the stimulated region changed from elliptical to circular as the fluid viscosity increased from 1 to 5 mPa·s. The size of the stimulated region decreased with the increasing fluid viscosity. Additionally, more branches were induced under treatment with higher-viscosity fluids. The HF mainly propagated along the direction of σ1 for cases with a relatively low-viscosity fluid. However, no dominant propagation direction was observed for the high-viscosity cases of μ=4 and 5 mPa·s.Figure 14
HF propagation in the NFN with different fluid viscosities: (a) 1 mPa·s, (b) 2 mPa·s, (c) 3 mPa·s, (d) 4 mPa·s, and (e) 5 mPa·s.As shown in Figure15, Lz was greater than Lx (except for the case of μ=4 mPa·s), and the difference between Lz and Lx for the high-fluid viscosity cases of μ=3, 4, and 5 mPa·s was smaller than that for the lower-fluid viscosity cases of μ=1 and 2 mPa·s.Figure 15
Variations of the maximum vertical and lateral growth for different fluid viscosities.As shown in Figure16, α tended to increase with the fluid viscosity. A significant increase (53°) in α was observed when the fluid viscosity increased from 2 to 4 mPa·s. The deviation angle was maximized (approximately 90°) in the high-fluid viscosity cases of μ=4 and 5 mPa·s.Figure 16
Variation of the maximum deviation angle for different fluid viscosities.


## 3.4. Effect of NF Orientation

Five models were used for the simulations, with dip angles of 0°, 15°, 30°, 45°, and 60°. The angles between the NFs and the direction ofσ1 were 90°, 75°, 60°, 45°, and 30°, respectively. The stress state was as follows: σx=σy=5 MPa, σz=5.5 MPa. Simulations were performed under a constant injection rate of 0.0005 m3/s for 2 s.As shown in Figure17, the HFN geometry varied significantly with changes in the NF dip angle. The direction of the long axis of the stimulated region changed from subhorizontal to subvertical as the NF dip angle increased from 0° to 60°. As shown in Figure 14(a), the two HF branches were initiated at the wellbore and then encountered the two NFs closest to the wellbore. Fluid invasion into the two NFs occurred, causing the open and shear slip of the NFs and then extending from the edges of the NFs. This process occurred again when the HF encountered the next NF. Thus, a step-like HFN was formed in the case of θ=0°. With an increase in the dip angle (and a corresponding reduction in the angle between the NF and the direction of σ1), the number of primary branches decreased, and the main propagation direction of the HF network became closer to the direction of σ1. There were four primary branches in the case of θ=0°, whereas only two branches were observed for θ=60°.Figure 17
HF propagation in the NFN with different NF dip angles: (a) 0°, (b) 15°, (c) 30°, (d) 45°, and (e) 60°.As shown in Figure18, as the dip angle increased from 0° to 60°, Lz decreased from 3.06 to 0.84 m (by 2.22 m). Conversely, Lx increased from 1.66 to 4.72 m. In the case of θ=0°, Lx was greater than Lz, whereas Lz exceeded Lx for θ>15°. The maximum difference between Lz and Lx occurred for θ=60°.Figure 18
Variation of the maximum deviation angle for different NF dip angles.Taking caseθ=45° as an example for detailed investigation, Figure 19 shows the microcracks (tensile failure of intact rock) and slip event (shear failure in NF). The HFN comprised the tensile failure in the intact rock and the shear failure in the connected NFs. The connected NFs appeared to be fully activated by shear slippage. Additionally, isolated NFs were subjected to shear failure, and slip events were observed on isolated NFs.Figure 19
Microcracks in intact rock and shear slip on NFs forθ=45°.The displacement field is shown in Figure20. Three primary branches divided the intact block into three separate blocks, which were characterized by different displacement fields. In general, the displacement decreased with the increasing distance from the wellbore. Most of the block experienced displacement of >0.2 mm. The largest displacement was observed for block 3.Figure 20
Discontinuous displacement field forθ=45°.


## 3.5. HF Propagation in Different NFNs

We investigated the HF propagation behavior in different NFNs and considered three types of simplified NFNs. The type A NFN comprised several large, parallel-distributed NFs (Figure21). In the type B NFN, there were vertical NFs between adjacent horizontal NFs, and all the NFs were isolated from each other (Figure 22). In type C NFN, the vertical NFs were connected with the horizontal NFs (Figure 23). The rock mass volume was 9m×9m×0.6 m. A starter crack with a radius of 0.15 m was located at the center of the model, normal to the x-axis. The stress state was as follows:σz=13MPa, σx=σy=5. Simulations were conducted under the injection rate of Q=0.003 m3/s for 12 s.Figure 21
Type A NFN: (a) model setup, (b) fracture aperture, (c) fracture fluid pressure, and (d) shear slip on the NF.
(a)(b)(c)(d)Figure 22
Type B NFN: (a) model setup, (b) fracture aperture, (c) fracture fluid pressure, and (d) shear slip on the NF.
(a)(b)(c)(d)Figure 23
Type C NFN: (a) model setup, (b) fracture aperture, (c) fracture fluid pressure, and (d) shear slip on the NF.
(a)(b)(c)(d)As shown in Figure21(b), a vertical HF was induced and crossed all the horizontal NFs, yielding a fishbone-like HFN. In general, the aperture of the horizontal NF decreased with the increasing distance from the starter crack (the injection point). For instance, the two horizontal NFs closest to the injection point had greater apertures (>1 mm) than the other horizontal NFs. Note that the aperture is varied in along the horizontal NF plane. Figure 21(c) shows the approximately radial distribution of the fracture fluid pressure, which decreased with increasing distance from the injection point. Figure 21(d) shows the shear slip on the horizontal NFs. Among the NFs, the largest shear slip was observed on the two horizontal NFs closest to the injection point.As shown in Figure22(b), the vertical flow pathway contained the newly induced HF and several vertical NFs. The flow network did not connect most of the vertical NFs. No significant difference in the fracture aperture was observed in comparison with type A. The approximately radial distribution pattern of the fracture fluid pressure is presented in Figure 22(c). The fluid pressure varied slightly compared with type A. The induced shear slip on the horizontal NFs was similar between types A and B. Additionally, the vertical NFs connected by the HF appeared to be fully activated by shear failure. Most of the other vertical NFs also underwent a slight slip. In general, the HF propagation was similar between types A and B.The HF propagation for type C differed significantly from the two aforementioned cases. As shown in Figure23(b), the fluid was forced to pass through the connected NFN, and only a few new HFs were induced. The HF extended from the tip of the vertical NF, which was located on the edge of the NFN. The dominant flow path was still along the preferential fracture plane (the direction of maximum principal stress). As shown in Figure 23(c), the distribution pattern of the fracture fluid pressure was approximately radial. However, the fluid pressure gradient was reduced in comparison with those for types A and B, owing to the high connectivity of the NFN. Specifically, the size of the region with high fluid pressure (>20 MPa) was smaller than those for types A and B. The fluid pressure on the edge of the network exceeded approximately 8 MPa. The number of vertical NFs that experienced shear slip was significantly larger than that for type B (Figure 23(d)). A large shear zone along the vertical direction was simulated, in which the vertical NFs appeared to be fully activated by shear slippage.


## 4. Discussion

### 4.1. Analysis of Geometry of HFN

The objective of this study was to examine the effects of various parameters on the HFN geometry and the complexity of the stimulation mechanism in naturally fractured formations. It has been well documented that the HF propagation pathway is dominated by in situ stress. According to the fracture mechanism and HF theory, it was argued that a larger stress difference leads to a shorter distance between the reoriented HF and the location of the maximum principal stress [38]. Liu et al. [39] concluded that the HF always propagates along the path of least resistance, regardless of the direction of HF initiation. Additionally, an NF is opened or dilated when the HF propagates to the intersection point, resulting in a complex HFN with an elliptical stimulated region. Figure 24 shows a schematic of several types of behavior that have been observed in minebacks or laboratory tests. Despite the tortuosity of propagation pathway, HF always propagates along the path of least resistance [40]. Those interaction behaviors (e.g., crossing, diversion, and arrest) have been well reproduced in XSite simulations.Figure 24
Schematic of types of observed fracture behavior crossing interfaces [40].As shown in Figure25, at the local scale, fluid tends to follow the NF rather than induce a new HF in the intact rock, as it must minimize the local work. However, on a large scale, the global orientation of the HF tends to remain parallel to σ1 owing to the global work minimization requirements [41].Figure 25
Schematic of the global orientation of fracture propagation in naturally fractured formation [41].As shown in Figure26, in a large size triaxial experiment, Chen et al. [42] found that the fracture network evolution pattern is dependent on the horizontal principal stress difference. Under small stress difference, a radial fracture network would be induced. In contrast, under large stress difference, a dominate fracture with some small multibranch fractures would occur.Figure 26
Experimental results of HF propagation pattern under different horizontal stress difference: (a) 5 MPa and (b) 10 MPa [42].
(a)(b)In our models, multiple radial HFs can be initiated at the wellbore, and an HFN with complex geometry can be induced owing to the HF–NF interactions (e.g., crossing, diversion, and arrest). In general, the global orientation of HF propagation tended to remain parallel to the direction of maximum principal stress. As shown in Figure27, the HFN under a small stress difference was characterized by multiple radial branches, which are evenly distributed. With the increasing stress difference, the number of HF branches decreases, and the global orientation of the HF propagation is restricted to the direction of σ1. However, with an increase in the fluid viscosity or injection rate, the evolution of the HFN geometry exhibits the opposite trend.Figure 27
Evolution of the HFN geometry.Treatment with a sufficiently high-viscosity fluid prevents the fluid from leaking into the surrounding rock, which results in a significant stress concentration. As the injection rate increases, the leakage into the surrounding rock becomes less important [43]. Moreover, poroelastic stress changes can locally modify the given tectonic stress regime [44, 45]. Consequently, with a high-viscosity fluid/high injection rate, multiple radial HFs can be initiated from the wellbore, and the HF pathway becomes relatively independent of the far-field stress field. Variations in the NFN, e.g., in the NF dip angle or the connectivity, can significantly alter the geometry of the corresponding HFN. The HF–NF interaction behaviors, such as crossing, diversion, and arrest, lead to high complexity in the analysis of the HFN.The in situ stress dominates the HF propagation pathway on a large scale, whereas the treatment parameters (fluid viscosity/injection rate) and the NFs can alter the HF path on a local scale. The in situ stress, treatment parameters, and the NFs may act together to describe the propagation of HFN. The dominant stimulation mechanism in the NFN is determined by the contributions of different factors, which can vary significantly.


### 4.2. Analysis of Stimulation Mechanism in NFN

As shown in Figure28, McClure and Horne [46] presented four conceptual models for stimulation mechanisms: pure-opening model (POM), pure-shear stimulation (PSS), primary fracturing with shear stimulation leak off (PFSSL), and mixed-mechanism stimulation (MMS). The POM model assumes that no shear slippage occurs on the NF. The PSS model assumes that stimulation occurs through shear slippage on the NF and is hardly affected by the propagation of the new HF. The PFSSL model assumes that the continuous HF grows away from the wellbore, and fluid leaks into the connected NFs, resulting in shear slippage. The MMS model assumes that the HF can be terminated against the NF. This inhibits the development of the continuous HF, forcing the fluid to pass through a network consisting of the new HF and NF.Figure 28
Four proposed stimulation mechanisms. The new HFs and the NFs are indicated by red and blue lines, respectively. The mechanisms are POM, PSS, PFSSL, and MMS [44].POM and PSS are extreme cases that are unlikely to occur for naturally fractured formations. PFSSL and MMS are more probable for field treatment. The variations in the ratio of the tensile fracture (new HF) to the shear fracture (shear slip on NF) can be significant. They depend on the stress state, treatment parameters, and NFN, which determine the dominant stimulation mechanism. For instance, under a small stress difference, multiple radial HFs grow from the wellbore, connecting more NFs and forming a larger stimulated region compared with the large-stress difference case (see Figure8). Additionally, the NF is more likely to experience shear slip under smaller stress difference. Therefore, the shear slippage may play a significant role in the stimulation mechanism under a small stress difference. Conversely, the new induced HF may play a major role under a large stress difference. With a larger stress difference, the stimulated region is smaller, and the HF propagation is restricted to the direction of σ1. Thus, fewer NFs can be reached by the HF, and the NFs are more stable. The dominant stimulation mechanism can be altered (at least to some extent) by changing the treatment parameters (injection rate and fluid viscosity) to modify the HFN geometry (see Figures 11 and 14). The type of NFN also affects the simulation mechanism. For a rock mass that contains well-connected NFs, the fluid tends to follow the NFN, and only a few new HFs are induced. In this case, shear slippage plays a dominant role in the stimulation of the NFN (see Figure 23). In contrast, for rocks containing isolated NFs, the formation of the continuous HF is more dominant.Quantitative evaluation of the stimulated region remains difficult. As shown in Figure29, Chen et al. [41] reported a discrepancy between the stimulated (dilated) zone and the sand zone. The shear dilation effect can radiate outward by 200–300 m in some cases, resulting in a stimulated zone around the sand zone that has a larger volume than the sand zone.Figure 29
Schematic of the stimulated (dilated) zone around the sand-propped zone [41].Neither the analysis of McClure and Horne nor that of Dusseault considered the isolated NF, which may experience shear slip. In the present study, there was a dilated zone involving unconnected NFs, which experienced shear slip. As shown in Figure19, many unconnected NFs also underwent shear slip. The slip on unconnected NFs may be attributed to poroelastic stress changes in the rock mass, which is supported by microseismicity interpretation [47] and theoretical analysis [48]. Even though activation of unconnected NFs cannot enhance the conductivity of the HFN, the stress acting on the NF plane and the elastic strain energy can be reduced, which may be beneficial for mitigating rockburst or destressing underground excavations [49]. As shown in Figure 20, the variation in the displacement in a large area around the HF was observed. The displacement field can reflect the redistribution of stress. The stimulated region may contain the unconnected NFs experienced shear slip, modified displacement field, and redistribution of stress; it is not limited to the connected HFN. The stimulated region is larger than typically acknowledged, and a highly complex stimulation mechanism can be expected in naturally fractured formations.


### 4.3. Modelling Considerations and Future Work

The assumed NFN geometry might not precisely represent the actual complex NFN, which may consist of irregularly oriented NFs with different mechanical properties. The assumed NFN is a reasonable simplification because it properly accounts for the interaction behaviors between the HF and multiple NFs and can thus be used to predict an HFN with a simplified geometry. A discrete fracture network (DFN), which is based on geological mapping, stochastic generation, and geomechanical simulation, is recommended for a more realistic representation of the NFN [50, 51]. A DFN can explicitly represent the geometric properties of individual NFs (e.g., the size, position, orientation, shape, and aperture), as well as the topological relationships between individual NFs and NF set [52]. Numerical simulations have contributed significantly to our understanding of HF propagation in fractured formation. However, numerical methods, such as the XFEM, BEM, and DEM, have their own merits and limitations [22–24]. Hydraulic fracturing is a nonlinear and multiscale process that involves mechanical deformation, fluid flow, fracture propagation, and their interaction. Additionally, the mechanical uncertainty and spatial variability of naturally fractured formations present considerable challenges for sensitivity and risk analyses [53, 54]. In addition to numerical simulations, support from methodologies, analysis, and experiments are required for a clearer understanding of the formation of HFN in naturally fractured formations.


## 4.1. Analysis of Geometry of HFN

The objective of this study was to examine the effects of various parameters on the HFN geometry and the complexity of the stimulation mechanism in naturally fractured formations. It has been well documented that the HF propagation pathway is dominated by in situ stress. According to the fracture mechanism and HF theory, it was argued that a larger stress difference leads to a shorter distance between the reoriented HF and the location of the maximum principal stress [38]. Liu et al. [39] concluded that the HF always propagates along the path of least resistance, regardless of the direction of HF initiation. Additionally, an NF is opened or dilated when the HF propagates to the intersection point, resulting in a complex HFN with an elliptical stimulated region. Figure 24 shows a schematic of several types of behavior that have been observed in minebacks or laboratory tests. Despite the tortuosity of propagation pathway, HF always propagates along the path of least resistance [40]. Those interaction behaviors (e.g., crossing, diversion, and arrest) have been well reproduced in XSite simulations.Figure 24
Schematic of types of observed fracture behavior crossing interfaces [40].As shown in Figure25, at the local scale, fluid tends to follow the NF rather than induce a new HF in the intact rock, as it must minimize the local work. However, on a large scale, the global orientation of the HF tends to remain parallel to σ1 owing to the global work minimization requirements [41].Figure 25
Schematic of the global orientation of fracture propagation in naturally fractured formation [41].As shown in Figure26, in a large size triaxial experiment, Chen et al. [42] found that the fracture network evolution pattern is dependent on the horizontal principal stress difference. Under small stress difference, a radial fracture network would be induced. In contrast, under large stress difference, a dominate fracture with some small multibranch fractures would occur.Figure 26
Experimental results of HF propagation pattern under different horizontal stress difference: (a) 5 MPa and (b) 10 MPa [42].
(a)(b)In our models, multiple radial HFs can be initiated at the wellbore, and an HFN with complex geometry can be induced owing to the HF–NF interactions (e.g., crossing, diversion, and arrest). In general, the global orientation of HF propagation tended to remain parallel to the direction of maximum principal stress. As shown in Figure27, the HFN under a small stress difference was characterized by multiple radial branches, which are evenly distributed. With the increasing stress difference, the number of HF branches decreases, and the global orientation of the HF propagation is restricted to the direction of σ1. However, with an increase in the fluid viscosity or injection rate, the evolution of the HFN geometry exhibits the opposite trend.Figure 27
Evolution of the HFN geometry.Treatment with a sufficiently high-viscosity fluid prevents the fluid from leaking into the surrounding rock, which results in a significant stress concentration. As the injection rate increases, the leakage into the surrounding rock becomes less important [43]. Moreover, poroelastic stress changes can locally modify the given tectonic stress regime [44, 45]. Consequently, with a high-viscosity fluid/high injection rate, multiple radial HFs can be initiated from the wellbore, and the HF pathway becomes relatively independent of the far-field stress field. Variations in the NFN, e.g., in the NF dip angle or the connectivity, can significantly alter the geometry of the corresponding HFN. The HF–NF interaction behaviors, such as crossing, diversion, and arrest, lead to high complexity in the analysis of the HFN.The in situ stress dominates the HF propagation pathway on a large scale, whereas the treatment parameters (fluid viscosity/injection rate) and the NFs can alter the HF path on a local scale. The in situ stress, treatment parameters, and the NFs may act together to describe the propagation of HFN. The dominant stimulation mechanism in the NFN is determined by the contributions of different factors, which can vary significantly.


## 4.2. Analysis of Stimulation Mechanism in NFN

As shown in Figure28, McClure and Horne [46] presented four conceptual models for stimulation mechanisms: pure-opening model (POM), pure-shear stimulation (PSS), primary fracturing with shear stimulation leak off (PFSSL), and mixed-mechanism stimulation (MMS). The POM model assumes that no shear slippage occurs on the NF. The PSS model assumes that stimulation occurs through shear slippage on the NF and is hardly affected by the propagation of the new HF. The PFSSL model assumes that the continuous HF grows away from the wellbore, and fluid leaks into the connected NFs, resulting in shear slippage. The MMS model assumes that the HF can be terminated against the NF. This inhibits the development of the continuous HF, forcing the fluid to pass through a network consisting of the new HF and NF.Figure 28
Four proposed stimulation mechanisms. The new HFs and the NFs are indicated by red and blue lines, respectively. The mechanisms are POM, PSS, PFSSL, and MMS [44].POM and PSS are extreme cases that are unlikely to occur for naturally fractured formations. PFSSL and MMS are more probable for field treatment. The variations in the ratio of the tensile fracture (new HF) to the shear fracture (shear slip on NF) can be significant. They depend on the stress state, treatment parameters, and NFN, which determine the dominant stimulation mechanism. For instance, under a small stress difference, multiple radial HFs grow from the wellbore, connecting more NFs and forming a larger stimulated region compared with the large-stress difference case (see Figure8). Additionally, the NF is more likely to experience shear slip under smaller stress difference. Therefore, the shear slippage may play a significant role in the stimulation mechanism under a small stress difference. Conversely, the new induced HF may play a major role under a large stress difference. With a larger stress difference, the stimulated region is smaller, and the HF propagation is restricted to the direction of σ1. Thus, fewer NFs can be reached by the HF, and the NFs are more stable. The dominant stimulation mechanism can be altered (at least to some extent) by changing the treatment parameters (injection rate and fluid viscosity) to modify the HFN geometry (see Figures 11 and 14). The type of NFN also affects the simulation mechanism. For a rock mass that contains well-connected NFs, the fluid tends to follow the NFN, and only a few new HFs are induced. In this case, shear slippage plays a dominant role in the stimulation of the NFN (see Figure 23). In contrast, for rocks containing isolated NFs, the formation of the continuous HF is more dominant.Quantitative evaluation of the stimulated region remains difficult. As shown in Figure29, Chen et al. [41] reported a discrepancy between the stimulated (dilated) zone and the sand zone. The shear dilation effect can radiate outward by 200–300 m in some cases, resulting in a stimulated zone around the sand zone that has a larger volume than the sand zone.Figure 29
Schematic of the stimulated (dilated) zone around the sand-propped zone [41].Neither the analysis of McClure and Horne nor that of Dusseault considered the isolated NF, which may experience shear slip. In the present study, there was a dilated zone involving unconnected NFs, which experienced shear slip. As shown in Figure19, many unconnected NFs also underwent shear slip. The slip on unconnected NFs may be attributed to poroelastic stress changes in the rock mass, which is supported by microseismicity interpretation [47] and theoretical analysis [48]. Even though activation of unconnected NFs cannot enhance the conductivity of the HFN, the stress acting on the NF plane and the elastic strain energy can be reduced, which may be beneficial for mitigating rockburst or destressing underground excavations [49]. As shown in Figure 20, the variation in the displacement in a large area around the HF was observed. The displacement field can reflect the redistribution of stress. The stimulated region may contain the unconnected NFs experienced shear slip, modified displacement field, and redistribution of stress; it is not limited to the connected HFN. The stimulated region is larger than typically acknowledged, and a highly complex stimulation mechanism can be expected in naturally fractured formations.


## 4.3. Modelling Considerations and Future Work

The assumed NFN geometry might not precisely represent the actual complex NFN, which may consist of irregularly oriented NFs with different mechanical properties. The assumed NFN is a reasonable simplification because it properly accounts for the interaction behaviors between the HF and multiple NFs and can thus be used to predict an HFN with a simplified geometry. A discrete fracture network (DFN), which is based on geological mapping, stochastic generation, and geomechanical simulation, is recommended for a more realistic representation of the NFN [50, 51]. A DFN can explicitly represent the geometric properties of individual NFs (e.g., the size, position, orientation, shape, and aperture), as well as the topological relationships between individual NFs and NF set [52]. Numerical simulations have contributed significantly to our understanding of HF propagation in fractured formation. However, numerical methods, such as the XFEM, BEM, and DEM, have their own merits and limitations [22–24]. Hydraulic fracturing is a nonlinear and multiscale process that involves mechanical deformation, fluid flow, fracture propagation, and their interaction. Additionally, the mechanical uncertainty and spatial variability of naturally fractured formations present considerable challenges for sensitivity and risk analyses [53, 54]. In addition to numerical simulations, support from methodologies, analysis, and experiments are required for a clearer understanding of the formation of HFN in naturally fractured formations.


## 5. Conclusions

The lattice-spring code XSite was employed to determine the effects of various parameters on the geometry of the HFN. Sensitivity analyses were performed to investigate several controlling factors, such as the stress difference, injection rate, fluid viscosity, and NF orientation. The HF propagation in three types of HFNs was analyzed. According to the results, the following conclusions are drawn:(1)
HF propagation tended to remain parallel to the direction of maximum principal stress. The in situ stress significantly affects the global orientation of the fracture propagation on a large scale. In contrast, NFs can alter the fracture pathway on a local scale owing to HF crossing, diversion, or arrest(2)
With a large stress difference, the global orientation of fracture propagation is restricted to the direction of maximum principal stress, changing the geometry of the HFN from circular to elliptical. With an increase in the fluid viscosity or injection rate, the evolution of the HFN geometry exhibits the opposite trend(3)
The growth of multiple branches and the complexity of the HFN are reduced under a large stress difference. Conversely, a high injection rate and fluid viscosity contribute to the growth of multiple branches from the borehole and the complexity of the HFN. With a reduction in the angle between the NF and the maximum principal stress, the fracture branches and the complexity of the HFN are reduced(4)
An NFN with higher connectivity tends to induce a larger shear zone and smaller fluid pressure gradient. The variations in the ratio of the tensile fracture (HF) to the shear fracture (shear slip on NF) can be significant. They depend on the stress state, treatment parameters, and NFN, which determine the dominant stimulation mechanism



---

*Source: 8845990-2020-10-24.xml*