Simulation of the interactions between multiple hydraulic fractures and natural fracture network based on Discrete Element Method numerical modeling

Obtaining the most Stimulated Reservoir Volume (SRV) by forming a complex fracture network is the prerequisite for shale gas hydraulic fracturing. In this paper, a coupled DFN–DEM model is presented based on a realistic discrete fracture network from fracturing treatment, to simulate the complex interactions between a created hydraulic fracture and a natural fracture system in shale. The sensitivity analysis of numerical simulation was performed to investigate the effects of hydraulic fracture length, multi‐fracture spacing, and internal friction angle on the stimulation reservoir area. The simulation results revealed that regardless of the DFN distribution and the length of hydraulic fracture, shear failure always tracked on the existing natural fracture plane. The amount and depth of the natural fracture shear into formation increased as the hydraulic fracture grew in length. Moreover, the actual dynamic morphology of fracture network was the radial or dendritic. A drastic change of permeability along the distance from injection may stimulate the fracture networks especially near the wellbore. For the multi‐fracture, the opening of the second fracture was able to stabilize the sheared fractures at the tip of the initial fracture. The shear failure area around the second fracture tended to increase with the growth of the length of the second fracture. Fracture network connectivity in the system is better for the optimal fracture spacing, and then fluid penetration during injection will influence the permeability evolution during stimulation. The stimulation reservoir area of the fracture network and the natural fracture friction angle had a quadratic relationship. Finally, the results by the data from field measurement in Sichuan, China verified that the proposed coupled DFN–DEM model had an accuracy up to 91.7%. The research results provide a reference for predicting the development of fracture network far away for multistage fracturing treatment.


| INTRODUCTION
The exploration and development of unconventional resources are greatly accelerating as the conventional natural resources are running out. Efficient exploitation of shale gas and oil has already obtained gigantic achievements in North America, such as Marcellus Shale, Barnett Shale, Lewis Shale, Fayetteville Shale, and Woodford Shale. [1][2][3][4][5][6][7] Obviously, shale gas has received much attention from the world's energy industry.
Artificial hydraulic fracturing, as one of the most frequently used technique in the petroleum industry for stimulating various unconventional applications, has been widely used to increase the permeability of a shale reservoir by the injection of a pressurized fluid into a formation. [8][9][10] The formation of complex fracture networks is essential for the cost-effective exploitation and high-efficient production of shale gas. Therefore, efforts should be made to address the nonlinear behaviors of naturally fractured reservoirs due to the increasingly intensive development of unconventional reservoirs. Exploring the main controlling factors of fractures formed in the process of shale reservoir fracturing should consider the discontinuity characteristics of shale. In most naturally fractured reservoirs, fractures are the primary flow channels, and rock matrix only provides the primary storage capacity. The existence of the pre-existing natural fractures should be considered in both stimulation and production processes, which is the key to shale development. [11][12][13] The fundamental point for the modeling and design of shale gas hydraulic fracture stimulations is that, in the process of fracturing, how a created hydraulic fracture can induce natural fractures to open or propagate, realizing a better communication between an artificial hydraulic fracture and a naturally fractured formation, to form a fracture network system with high fracture flow conductivity, and to increase the stimulation reservoir area of fracture network to the utmost, as illustrated in Figure 1.
The key problem of the complex interaction between a created hydraulic fracture and a natural fracture system in shale is of great interest for the energy resource industry because the formation of a fracture network system can significantly influence the effectiveness of hydraulic fracturing.
Even though the formation of complex fracture network has been extensively investigated, some problems require further research. Firstly, natural fracture and artificial fracture were studied separately, with their coupling effects ignored. [14][15][16] Secondly, the natural fractures in shale reservoirs are characterized by disordered distribution. Normally, the formation of fracture network is related with the disorderly development of the natural fractures. [17][18][19] Finally, micro-fractures are an effective tool to open the migration pathways of oil and gas, thereby greatly enhancing the productivity. Actually, the fractures identified by naked eyes are limited in quantity for shales, and that, more large cracks or fractures in the development may lead to lower gas production rate, which indicates that a large crack appears not to be conducive to the preservation of shale gas, but micro-fractures. 20,21 As the micro-fractures of shales are extremely developed, the closed natural fractures due to cementation will be induced open through hydraulic fracturing. Numerous papers have reported the research on hydraulic fracturing in shale. So far, nonlinear and discontinuous numerical methods have become a significant tool in modeling geotechnical engineering failure processes. 22 A variety of methods have been proposed to simulate the hydraulic fracturing propagation and interaction with the pre-existing natural fractures. One of the foundational work on the naturally fractured reservoirs was conducted in 1962 by Warren and Root. 23 Moreover, due to the huge increase in unconventional reservoir development, Cundall, and Strack proposed the Discrete Element Method (DEM) to study the mechanical characteristics of rocks and soils, which was more versatile and computationally intensive than dual porosity models. 24,25 Safariand Ben et al, employed the displacement discontinuity method (DDM) techniques to model the hydraulic fracturing in naturally fractured formations. 26,27 Considering the effects of nonlinear normal deformation and the shear dilation of fractures, Ki-Bok Min F I G U R E 1 Illustration of the communication between an artificial hydraulic fracture and a naturally fractured formation to form a fracture network system 2924 | LI conducted a series of numerical simulations to calculate the permeability of simulated fractured rock masses using a 2D DEM. 28 Qiao proposed a method for fluid-solid coupling using the varying joint stiffness method (VJSM) based on the characteristic that the compressibility of a fracture decreases with the stress increasing. 29 To predict the conductivity during fracture closing, Zhu employed a new method in channel fracturing by implementing analytical solution in DEM. 30 Nagel employed both 2D and 3D discrete element block models to evaluate the simulation results of the natural fracture shear in the hydraulic fracturing of naturally fractured formations. 31,32 The natural fractures were prone to shear failure when there was an interference between the hydraulic fracture and the natural fracture. However, the simulation results by Nagel were rather discrete, and the internal laws between the basic parameters of rock mechanics and shear failure area or SRV have not been systematically demonstrated.
Considering the existence of a large number of random natural fractures and joints in shale, mathematical analytical models cannot be effectively applied to modeling. Also, finite element numerical method is not suitable for the modeling case where lots of natural fractures exist, due to the unsatisfactory calculation speed and convergence performance. By contrast, discrete element numerical simulation methods have certain advantages in dealing with discontinuous media-related problem, as they can form a random discrete fracture network (DFN) to directly simulate the interaction between hydraulic fractures and natural fractures. In this study, a coupled DFN-DEM model is presented based on a realistic discrete fracture network from fracturing treatment, to simulate the complex interactions between a created hydraulic fracture and a natural fracture system in shale. Simultaneously, this research investigates the influences of hydraulic fracture length, multi-fracture spacing, and internal friction angle on the SRV, fracture paths, and patterns due to the creation of a single fracture alone and parallel fracture via sensitivity analysis. The numerical simulation results provide effective reference for revealing the disturbance and expansion behaviors of natural fractures during hydraulic fracturing in the far wellbore zone. The research results provide novel insights for the design of shale gas fracturing.

| Joint deformation model
A rock is divided into multiple block elements by joints and fissures, and each block element which is quasi-rigid maintains its shape and size during the computing process. Also, each block element can contact with neighboring nodes in the process of movement due to the control of joints and other discontinuities, also can be apart by removing the tangential and normal forces of contact point. For those fractures with two rough surfaces to contact each other in the fracture network, Figure 2 illustrates a fracture subject to normal and shear stress between them.
A fracture normal deformation equation induced from a large body of experimental data is given by (Bandis et al) 33 where ′ n is the effective normal stress, D n is the normal closure of fracture, aa = 1/K ni , aa/b = D n max .
So, Equation (1) can be rewritten as.
where K n is the normal stiffness, derived as a function of D n max or ′ n or For the i-th fracture, the effective normal displacement ΔD n and the effective normal stress Δ n of the fracture satisfy where K n is related to normal stress, a n and e n are the model parameters.
Similarly, there is a linear relationship between the change in shear stress and the change in shear displacement before yielding. Shear displacement is also related to shear stiffness K S , as expressed by An approximate linear relation is used for the aperture increase ΔD n−dilation due to shear movement.
where d is the dilation angle.

| The permeability of fractured network
The rock mass cut by natural joints or fractures is a dual medium, the rock is a porous medium, and the network system composed of fractures is a discrete medium. For most of rocks, the permeability coefficient of a complete block is extremely weak. In the process of hydraulic fracturing, fracturing fluid is injected at a certain flow rate, and then flow along the fracture network interface only. In this process, the fracture conductivity and mechanical deformation interact with each other.
The flow analysis of an idealized fracture intersection in a fracture network is based on the domain structure of fractures, as illustrated in Figure 3. The intersections in the figure are the locations where two or more fractures meet, and the fracture between two adjacent intersections is called a segment. Then, the sets of fracture segments form a fracture network for conducting flow. 34 The numerical implementation for fluid flow effects is illustrated in Figure 4. For a closely packed network system of domains, each domain is assumed to be filled with fluid at uniform pressure and communicates with its neighbors through contacts. The features of fluid flow in the fractures through the domain are determined by the pressure difference between the adjacent fluid domain, and the fluid pressure depends on the value at the center of the domain. Referring to Figure 2, domains 1, 3, and 4 represent joints, domain 2 is located at the intersection of two joints, and domain 5 is a void space. The letters A to F represent the contact points where the forces of mechanical interaction between blocks are applied.
The cubic law for flow in a planar fracture (Witherspoon et al 1980) can be used. 35 The fluid balance equation is expressed as.
(6) K n = a n e n n

F I G U R E 3 An idealized fracture intersection in two dimensions
for flow

F I G U R E 4 Flow in joints modeled as flow between domains
LI where q f is the flow rate in the fracture per unit, q int is the interface flow rate per fracture length per unit, ρ f is the fluid density, ΔL is the length of fracture segment and q s is the production rate per unit. The fluid flow in natural joints consists of tangential flow and normal flow. The volume flow rate of tangential flow is defined by the pressure gradient along the fracture surface. 36 In this case, Darcy's law is still applicable. The flow rate in the fracture can be derived by where k f is the fracture permeability and can be obtained from the Cubic law, and w f is the fracture aperture where c f is the fluid compressibility. Combining Equations (10) and (11), the net fracture flow rate is expressed by The fracture volume is directly related to the change of fracture aperture Equation (13) can be rewritten as Substituting Equation (11) into Equation (15) yields Combining Equations (9)-(16), the fluid diffusion equation in the fracture network is obtained

Method (DEM)
The DEM is usually used in dealing with the noncontinuum medium (eg, joints and fractures in rock mass), and it is based on the Lagrangian algorithm that is well-suited to model the large movements and deformations of a blocky system. Early in 1971, the American scholar P. Cundall first raised the DEM law. 25 At present, the Universal Distinct Element Code (UDEC), a 2D numerical calculation program based on the DEM for discontinuous modeling, has been widely used to simulate the response of a jointed rock mass. 37 The essence of the method is the step-by-step integration for the critical damping vibration equation, which can turn a nonlinear static problem into a dynamic problem. In order to guarantee the accuracy, the quality and stiffness damping are used to absorb the kinetic energy of the system. In a 2-D fracture, it is assumed that the fracture geometry is planar, and the intersections of fractures are assigned as nodes.
The major features of the DEM include: 1. Large deformation simulation (slip and extend) along the structural surface in discrete media; 2. Discrete media are expressed by rounded block aggregates, and discrete blocks are treated as deformed or rigid; 3. Stable solutions for the static and dynamic analysis of jointed rock masses; 4. Explicitly recognize new contacts automatically or interfaces between the discrete bodies, and allow for finite displacements, rotations of discrete bodies, fracture matric deformation, and failure including fracture propagation.

DFN-DEM
To explore whether there is interference between hydraulic fractures and natural fractures and whether natural fractures will connect to each other, it is necessary to observe the fracture propagation tendency and fracture morphology of the intersection area during the extending process of hydraulic fracturing. There are three main features compared with previous research: (a) the simulation considers joint/fracture hydro-mechanical coupling effects; (b) The discontinuity sets of coupled DFN-DEM model was developed based on the shale outcrops and the micro-seismic events from field measurement together. The fluid balance equation in the fracture includes the flow from the connected fractures and the interface flow from the connected matrices. The hydro-mechanical coupled mechanism of the interactions between hydraulic fractures LI and natural fracture network should be explained according to stress conditions and fracture geometries. Figure 5 shows the modeling procedures of such interactions in a fluid-mechanically coupled system using DFN-DEM model. Three components should be coupled simultaneously in hydro-mechanical modeling the hydraulic fracture propagation in naturally fractured formation: (a) the fracture deformation induced by fluid pressure; (b) the permeability of fractured network; and (c) the active fracture network. The geometrical basis for fracture network model with DEM is a DFN generated in a square region. The DFN-DEM model is constructed by the new FISH code to determine the stimulation reservoir area. As shown in Figure 5, the workflow consists of combining the dip direction, spacing, trace length, gap length, fracture density, fisher constant, and other fracture parameters used for DFN generation.

| Discontinuity sets and model boundary
The determination of fracture or joint occurrence is a prerequisite for the establishment of a DEM model, and different combinations of natural fractures determines different fracture network morphologies. According to some statistics, jointed rock faces occur in groups or individuals, and also appear in form of multiple intersections. Rough joints often have a significant status in a group or multiple sets of joints. Therefore, rough joints can be considered as the major failure surface in rock mass, and the rest of joints can be regarded as the secondary failure surface.
The key issue for modeling is to construct a fracture network system as real as possible. According to the observation of shale outcrops, natural joint inclination arrangement has a certain regularity, as shown in Figure 6. Through the field measurement in the Wei-201-HX reservoir, the probability distribution of natural fractures was obtained. Based on the direction and inclination of the seismic source breaking plane, the discrete statistics of natural joints or fracture directions in shale were obtained to get the modeling basis for the DFN model based on statistical results. The statistical data of the micro-seismic fractures in the Wei-201-HX reservoir are taken as an example for interpretation, as shown in Figure 7. Combining the distribution of natural joint inclination from the shale outcrops and the statistical data of the micro-seismic fractures from field measurement, the natural fractures with different combinations of dip angles were preset in this model, and the randomly generated fracture network is composed of F I G U R E 5 The modeling procedures for a fluid-mechanically coupled system in DFN-UDEC model 2928 | LI two groups of fractures with case A (0°/70°) and case B (45°/135°). As the effective flow channels are formed by the connected fractures, the disconnected fractures are ignored for efficient calculation. Adjacent natural fractures are interconnected in closed-loop way.
This research ignores the presence of irregular fractures on a small scale in the shale reservoir, the surface roughness, and the permeability of rock matrix, and only considers the major natural fractures. The Monte Carlo simulation is a statistical method for stochastic simulation. 38 The random distribution function of natural fractures generated by Monte Carlo simulations is embedded into the fracture generator, and a DFN model is established in which artificial fractures and natural fractures interfere with each other, as shown in Figure 8. The probability density function of the random variables that generate the sample sets is estimated mainly by giving statistical sample sets. Fisher distribution is the most commonly used probability density function for discontinuity orientations. According to the distribution characteristics of the natural fractures in shale reservoir, each group of natural fracture sets is given a reasonable assignment to the dip, length, spacing length and trace length, gap length, fracture density, and fracture sets orientation dispersion (Fisher constant K), respectively, as illustrated in Table 1.
It is assumed that the DFN 2D model with fracture network is 600 m × 600 m in this study. For the constructed model, the horizontal wellbore is along the X-axis (minimum horizontal principal stress σ h ), the bottom line of the model represents the horizontal wellbore wall, and the hydraulic fracture extends along the Y-axis perpendicular to the wellbore (maximum horizontal principal stress direction σ H ). Obviously, the distance between the two fractures represents the fracture spacing, and the hydraulic fracture length can be adjusted optionally along the Y-axis in the scope of the model. In the process of model building, the amplitude of the load is divided into

| 2929
LI 20 steps. The fracturing fluid enters the first hydraulic fracture after it is injected from the wellbore. The plane load of the fracturing fluid and the pore pressure on the internal surface of the first hydraulic fracture is applied to simulate the dynamic consolidation process step by step, and then followed by the addition of the second fracture. As the liquid pressure in the fracture increases gradually, the induced stress generated during the process of hydraulic fracture propagation makes the shear stress of the surrounding natural fractures exceed the shear strength and thereby cause failure. As a result, the hydraulic fractures and natural fractures interpenetrate to form a fracture network.

| Mechanical parameters
The distribution of large number of joints in rocks in a random way changes the mechanical properties of rock mass. Generally, the mechanical properties of rock mass mainly depend on the mechanical properties of intact rock and joints which have important influences on the stress, strain and strength of rock mass. Due to be large number of complicated natural fractures, shale reservoir exhibits strong heterogeneity characteristics. We set the mechanical parameters of the rock mass with reference to the model established by Nagel, 29 as listed in Table 2.
After establishing the fracture network numerical model, the main parameters of the reservoir are set as follows. The density of the formation fluid is 1.0 g/cm 3 , the vertical stress is 94.5 MPa, the maximum horizontal principal stress is 90.5 MPa, the minimum horizontal principal stress is 86.2 MPa, the initial pore pressure of the formation is 77.91 MPa, and the net pressure is 4.3 MPa.

| Influence of fracture length
The numerical simulation result shows that the stress concentration zone is distributed at the hydraulic fracture tip, which is favorable to make the fracture network grow along the trend of the natural fracture and extend to the far well zone, as shown in Figure 9. Moreover, the dynamic change in the hydraulic fracture length can be one of the effective criteria of the interference between the hydraulic fracture and the natural fracture. In this case, the half lengths of the preset hydraulic fractures in the numerical model are 50 m, 100 m, 200 m, and 300 m, respectively. Considering the natural fracture dip combinations of case A (0°/70°) and case B (45°/135°), the effects of fracture length on the expansion of the stimulation reservoir area are examined, as presented in Figures 10 and 11. The results reveal that the amount of sheared natural fracture increases as the hydraulic fracture grows in length. From a half-length fracture of 200-300 m, significant change in size is observed in the stimulation area, while the shear area is 150-250 m wide (perpendicular to the hydraulic fracture) into the formation. Moreover, two obvious behaviors can be observed, that is, under the induced role of hydraulic fracture, the shear failure of the fracture network would occur along the direction of the natural fracture plane, and the depth into the formation would increase as the hydraulic fracture grows in length. The significant conclusions from the result of single fracture DEM simulation indicate that the actual dynamic morphology of hydraulic fracturing is not the traditional plane double fracture, but the radial or dendritic fracture.
For case A (0°/70°) and case B (45°/135°), the stimulation reservoir area under different hydraulic fracture lengths is shown in Figure 12. The analysis clearly reveals that when the hydraulic fracture lengths are of the same size, the interconnected area of the natural fracture network in case B (45°/135°) is larger than that in case A. With the extension of hydraulic fractures, the difference in stimulation reservoir area between the two cases gradually increases. This phenomenon can be explained as follows. For larger natural fracture dip angle, the relative change in the induced stress in the direction of the minimum horizontal principal stress (perpendicular to hydraulic fracture) increases accordingly, and then the disturbance intensity increases, which leads to a larger probability of complex fracture network formation. Initial natural

Items Value
Intact rock fracture permeability K 0 = 0.0001 md, assuming fractured reservoir extension permeability is K. The evolution of permeability along the distance from injection with a halflength fracture of 50-300 m is shown in Figure 13. We observe that increasing the fracture half-length will lead to an increase in the spatial distribution of reservoir permeability enhancement within a certain range of distance from injection. At the beginning, the drastic changes of permeability along the distance from injection may stimulate the fracture networks especially near the wellbore. As fractures expand, permeability of reservoir gradually weakens and tends to reach a steady state due to the change in effective stress driven by fluid and influence on reservoir shear failure.

| Influence of multi-fracture spacing
Only considering the interference effect of a single hydraulic fracture on the surrounding reservoir stress field is not sufficient to describe the forming condition of the complex fracture network. The stress field produced in the extending process of another hydraulic fracture may affect the stress field around the original hydraulic fracture. With the adjacent hydraulic fracture opening and extending to further zones, the interference between the adjacent hydraulic fractures in the fracturing target zones is definitely crucial to the formation of a complex fracture network system. Therefore, a multi-fracture interference model was constructed to investigate the morphology after the hydraulic fracture interferes in the further zones for case A (0°/70°). Figure 14 shows the natural fracture shear results for an initial half-length hydraulic fracture 250 m, then for the addition of 100 m, 200 m and 300 m fracture. Figure 15 shows Take a combination of half-length fracture 200 m-250 m as an example, the optimal fracture spacing is investigated. Figure 16 shows the top view of the fracture propagation after the first and the second fractures interfering in further zones, with fracture spacings of 20 m, 25 m, 30 m, 60 m, and 80 m, respectively. The evolution of permeability of different fracture spacing for case A (0°/70°) is shown in Figure 17. The numerical results show the relationships between multi-fracture spacing and shear failure area. Specifically, as the fracture spacing decreases, the shear failure area diminishes. In addition, when the fracture spacing is 30 m, the natural fractures between the first and the second hydraulic fractures just intersect with each other to form a fracture network. When the fracture spacing is more than 30 m, the overall shear failure area further becomes larger. Moreover, the phenomenon in the case with fracture spacing shorter than 30 m, might be caused by the shear stress distribution change induced by the interaction between fractures while the second fracture propagates. In addition, we observe that there is a peak in the evolution of permeability when the fracture spacing is 30 m. Such behavior illustrates that when fracture network connectivity in the system is better for the optimal fracture spacing, and then fluid penetration during

| Influence of internal friction angle
The shear failure area of natural fracture is the connected area of fracture network which can be calculated by selecting the node coordinates of the failure boundary for fracture propagation area. The preset hydraulic fracture half-length in the numerical model is 200 m. We calculated the connected area of the natural fracture network under different internal friction angles (φ = 10°, 20°, 30°, 40°) for case A (0°/70°) and case B (45°/135°). The results are shown in Figure 18.
A large number of numerical calculation results are used to fit the function relation between the connected area and the internal friction angle for case A (0°/70°) and case B (45°/135°). The Influence of internal friction angle on the fracture network extension is illustrated in Figure 19. As can be seen from the figure, the connected area and the internal friction angle have a quadratic relation. The square of the correlation coefficient of the natural fracture angle for case A (0°/70°) and case B (45°/135°) is R 2 = 0.9928 and R 2 = 0.9945, respectively. It can be seen that the fitting results have high consistency with the actual results. Taking the internal friction angle of 20° as an example, the connected areas for case A (0°/70°) and case B (45°/135°) are 18 970 m 2 and 24 270 m 2 , respectively. Besides, when the internal friction angle is the same, the connected area for case B (45°/135°) is larger than that for case A (0°/70°).
The reasons for the above phenomena are as follows. As internal friction angle increases, the relative change in the induced stress along the direction of the minimum horizontal principal stress (perpendicular to hydraulic fracture) increases accordingly, and then the disturbance intensity increases, and the probability of forming a complex fracture network is relatively large. In addition, the calculation results show that the smaller the internal friction angle is, the larger the connected area is. This is because smaller internal friction angle results in the weakening of the shear strength of the rock, so that the shear fracture easily occurs, and the reservoir reformation is easy to form a complex fracture network.

| CASE ANALYSIS
Generally speaking, SRV is defined as the rock volume with increased permeability after hydraulic fracturing. Although SRV is not a measurement of production, it can be used to calibrate the DFN model to improve historical consistency and predict the Estimated Ultimate Recovery (EUR).Taking the field measurement data of the Wei 201-HX well in Sichuan, China as an example, the results of the permeability tensor of the fractured reservoir are presented in Figure 20. Yellow/orange/red/magenta indicates a permeability greater than 3000 Darcy, and green/blue indicates a permeability tensor below 3000 Darcy. If a fracture is contained in a grid cell in the volume, the grid cell is counted, even if the value of calculated permeability tensor is extremely small. All geological honeycomb units with permeability tensors are summed to obtain the total SRV of the stimulation well. The total SRV depends on the size of the modeling grid and can be modified according to the known reservoir flow characteristics. The calculations were also performed for the DEM-DFN model. The hydraulic fracture length is 200 m, hydraulic fracture average height is 120 m, the number of fracturing segments for a multi-stage hydraulic fracturing treatment in the Wei 201-HX well is 12. Figure 21 compares the field measurement and calculation results from DEM-DFN model. It is obvious that the numerical method in this study agrees well with the SRV from the field measurement data, with an average accuracy of 91.7%. So, it can be known that the proposed model can be used as the basis of choosing the stimulation reservoir area, and even the SRV.

| CONCLUSION
In this paper, based on coupled DFN-DEM numerical modeling, a series of numerical simulations have been conducted to evaluate the mutual interactions between a created hydraulic fracture and a natural fracture system in the shale formation. This paper draws the following conclusions: 1. The shear failure of fracture network stimulation occurred along the direction of natural fracture plane, and the stimulation reservoir area of the natural fracture shear into the formation would increase as the hydraulic fracture grows in length. The dynamic morphology of hydraulic fracture is not the traditional plane double fracture, but the radial or dendritic fracture. For large natural fracture angle, the disturbance intensity increases, which leads to a larger probability of forming a complex fracture network. Increasing the fracture half-length will lead to an increase in the spatial distribution of reservoir permeability enhancement within a certain range of distance from injection. 2. For the multi-fracture case, regardless of the DFN and hydraulic fracture length, the opening of the second fracture is able to stabilize the sheared fractures at the tip of the initial fracture. The connected area around the second fracture tends to increase with the growth of the length of the second hydraulic fracture. The natural fractures between the first and the second hydraulic fractures just intersect with each other to form a fracture network when the optimal fracture spacing is 30 m. Fracture network connectivity in the system is better for the optimal fracture spacing,