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

Coal Mining and Designing Branch, China Coal Research Institute, Beijing 100013, China Engineering Geology and Resource Geotechnics Group, Simon Fraser University, Burnaby, BC, V5A 1S6, Canada State Key Laboratory of Coal Mining and Clean Utilization, Beijing 100013, China Coal Mining and Designing Department, Tiandi Science and Technology Co., Ltd., Beijing 100013, China CCTEG Coal Mining Research Institute, Beijing 100013, China


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.
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][23][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. 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.
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]. Figure 4, the lattice is composed of numerous quasirandomly distributed nodes connected by springs. Joints are overlaid on the lattice using the SJM methodology.

Mechanical Formulation. As shown in
The central difference method is employed to compute the transitional degrees of freedom [32]: where _ u where ∑M ðtÞ i 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]: where N represents "normal," S represents "shear," F represents the spring force. k N and k S 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] If F n − pA < 0 then F n = 0, F s i = 0, where F n represents the normal force, F i S 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: if F n − pA + σ c A < 0 or jF s i j > τ c A, the bond fails in tension or shear (where σ c is the bond tensile strength, τ c is the bond shear strength), else, F s i ← F s i , the bond remains intact. Figure 5, 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.

Flow Formulation. As shown in
Lubrication theory is used to calculate the flow rate from node A to node B along a pipe: where a represents the aperture; μ f represents the fluid viscosity; p A and p B represents the hydraulic pressures at nodes "A" and "B," respectively; z A and z B 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, k r , is a function of saturation, s: The pressure increment, ΔP, during the flow timestep, Δt f , is calculated as 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 K f is the apparent fluid element bulk modulus.

Hydromechanical Coupling.
The fluid flow and mechanical process are fully coupled (see Figure 6). 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.
where E represents the Young's modulus. If K I < K IC (where K IC is the rock toughness); then, the spring tensile strength is utilized to detect spring failure. Otherwise, K I is compared to K IC to detect spring failure.   Table 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 = 5 MPa. 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 m 3 /s and a fluid viscosity of 1 mPa·s for 1 s. As shown in Figure 8, 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  5 Geofluids 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-stressdifference 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 9 shows the variations in the maximum vertical growth (L z ) and maximum lateral growth (L x ) of the HFN for the five models. With an increase in the stress difference, L x decreased from 4.62 to 2.49 m (by 2.13 m). Conversely, L z increased with an increase in the stress difference (except for the discrete points, the case of Δσ = 7:5 MPa was examined). For the case of Δσ = 0 MPa, L x was greater than L z . When the stress difference exceeded zero, L z was greater than L x . Moreover, the difference between L z and L x increased as the stress difference increased from 2.5 to 10 MPa. Figure 10 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.

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 m 3 /s. The volume of the injection fluid was 0.005 m 3 for all the models. The other parameters were as follows: σ x = σ y = 5 MPa, σ z = 7:5 MPa, and μ = 1 mPa·s.
As shown in Figure 11, 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 m 3 /s), and each branch was longer than the corresponding branch in the highest injection rate case (Q = 0:007 m 3 /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 m 3 /s).
As shown in Figure 12, the L z was greater than the L x for all five models, and among the models, the difference between the L z and the L x was the smallest for Q = 0:007 m 3 /s. As shown in Figure 13, α increased with the injection rate. There was no significant variation in α as the injection rate increased from 0.002 to 0.004 m 3 /s. However, there was a significant increase (53°) in α as the injection rate increased from 0.004 to 0.007 m 3 /s.

Geofluids
As shown in Figure 14, 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 higherviscosity 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.
As shown in Figure 15, L z was greater than L x (except for the case of μ = 4 mPa·s), and the difference between L z and L x 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.
As shown in Figure 16, α 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.
As shown in Figure 17, 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°.
As shown in Figure 18, as the dip angle increased from 0°t o 60°, L z decreased from 3.06 to 0.84 m (by 2.22 m). Conversely, L x increased from 1.66 to 4.72 m. In the case of θ = 0°, L x was greater than L z , whereas L z exceeded L x for θ > 15°. The maximum difference between L z and L x occurred for θ = 60°.
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.
The displacement field is shown in Figure 20. 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.

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 (Figure 21). 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 9 m × 9 m × 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 = 13 MPa, σ x = σ y = 5. Simulations were conducted under the injection rate of Q = 0:003 m 3 /s for 12 s.
As shown in Figure 21(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 Figure 22(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  The HF propagation for type C differed significantly from the two aforementioned cases. As shown in Figure 23(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.

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.
As shown in Figure 25, 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].
As shown in Figure 26, 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.
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 Figure 27, 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.
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.   Figure 28, 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 mixedmechanism 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.

Analysis of Stimulation Mechanism in NFN. As shown in
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 Figure 8). Additionally, the NF is more likely to experience

Geofluids
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 Figure 29, 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.
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 Figure 19, 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 Z X Y (a)  13 Geofluids 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.

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][23][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.

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

Data Availability
All relevant data used to support the findings of this study are included within the article.   Figure 29: Schematic of the stimulated (dilated) zone around the sand-propped zone [41].