The propagation behavior of hydraulic fracture network in a reservoir with cemented natural fractures

There are still some problems in the study of hydraulic fracture (HF) network evolution in cemented naturally fractured reservoirs, such as microseismic mapping showing exaggerated stimulated reservoir volume in some cases. In addition, the dominant role of natural fracture (NF) cementation strength, injection rate, in situ stress difference, NF distribution, and fracture initiation sequence of perforations in synthetically influencing fracture network formation needs to be further studied. For this purpose, a three‐dimensional matrix hexahedral element global coupled 0‐thickness cohesive element hydraulic fracturing model was developed. Results show that each interaction between HF and NF causes HF diameter shrinkage, which increases the propagation pressure of HF. When the cementation strength of the NF is low, the HF tends to deviate toward the tip of the NF to form a complex fracture network. Increasing the injection rate and the number of NFs can significantly enhance the complexity of the HF network, but does not change the HF and NF interaction pattern. The in situ stress differences dominate the morphology of the HF network when the cementation strength of NFs is constant. The stress interference of multiple fractures under segmented fracturing may form “S”‐shaped HFs, and the HFs are difficult to maintain a symmetrical morphology in the direction of the well axis. In addition, some NFs in inactivated damaged zones have developed a certain width geometrically due to the induced effect of HF, but they are still isolated by the low permeability matrix and might only generate some microseismic events.

(NFs), with oil and gas production exceeding the expected production of the shale reservoirs. 4 Therefore, the traditional bi-wing planar fracture model could not satisfy the demands of industry and complex hydraulic fracture (HF) network calculation.
Although NFs exist as a random distribution, most classical numerical simulations only simulated the interaction of NFs with HFs and ignored the effect of NF properties on the formation of the HF networks. These results showed that pre-existing NFs have a significant influence on the propagation of HFs and the formation of HF networks. 5 Both physical experiments and numerical simulations have observed three possibilities of interaction events between HF and pre-existing NF: opening, crossing, and arrested. 6,7 The viscosity and injection rate of the fracturing fluid also had a strong impact on the opening width of the NF and the geometry of the HF. 8 When multiple fractures overlap, the borehole injection pressure of neighboring wells will increase because of the compressive stress around the propagating fracture that squeezes the fluid within the fracture in the neighboring well. 9 In addition, the smaller interaction angle and the coefficient of friction on the NF surface allow the HF to be easily transferred to the NF. 10 Finally, the high injection pressure tends to drive the HF to cross the NF near the wellbore. 11 However, these studies only introduced stress interference between multiple HFs, and still do not fully understand the interaction between HF and cemented NF groups in a naturally fractured reservoir and the factors affected by the formation of the complex fracture network, for example, the impacts of cementation strength of NFs, in situ stress difference, injection rate, distribution characteristics of NFs, and stress interference with multiple fractures on the interaction performance when HF encounters multiple NFs, as illustrated in this paper.
Simulating HF propagation in the formation is a challenge, with at least three physical processes to be coupled simultaneously: (a) fracture initiation and propagation 12 ; (b) fluid flow within the fracture 13 ; (c) rock deformation due to fluid pressure on the fracture surface 14 ; with the upgrading of computer hardware, many new methods and models have been proposed to solve the numerical solution of complex coupled equations and thus capture the dynamic interaction process between HF and NF. An extended finite element method (XFEM)-based flow formulation was proposed by Shi and colleagues, [15][16][17][18] and the computer programs they prepared were able to simulate the interaction between HF and NF. Ghaderi et al. 19 combined the XFEM and the discrete element method (DEM) to investigate the HF propagation pattern in porous media with NF. They evaluate the production flow under different pore sizes and lengths of HF. Escobar et al. 20 used the XFEM to simulate the hydraulic fracturing process, considered the propagation of several vertically nonplanar fluid-driven fractures, and investigated the effect of stress shading on the pressure required for fracture propagation and the effect of fracture geometry. Tan et al. 21 established a numerical model for transverse isotropic laminated shale with a transition zone using the XFEM based on the cohesive zone model (CZM). The well-known DEM has also been used for modeling HF networks. Krzaczek et al. 22 combined the DEM with computational fluid dynamics to develop a numerical model to characterize the behavior of fluid-driven fractures in rocks. Yang et al. 23 developed a discrete element model (DEM) for synchronous hydraulic fracturing of two holes. The results showed that the azimuth and spacing affect the fracture propagation and fracture pressure, and the stress shadow effect between the two holes is easy to cause HF migration. The CZM is often used to construct models of complex fracture networks as well. Li 24 presented a twodimensional (2D) pore pressure CZM to simulate HF propagation in naturally fractured formation. Li et al. 25 and Zhai et al. 26 used the CZM to build a reservoir model with randomly distributed weak planes. In recent years, the phase-field method (PFM) has become increasingly popular in the fracturing community. Zhang et al. 27 combined XFEM and the PFM to solve discontinuous and continuous hydraulic fracturing formulations, respectively, for simulating hydraulic fracturing in a shale formation with frictional and cemented NFs. Xu et al. 28 developed a modified PFM to describe crack branching behavior with a wide range of fluid viscosities and injection rates. The branching behavior of the crack is quantitatively described by a phase diagram in which tension and shear failure modes alternately dominate it, followed by the periodic fluctuation of the crack tip velocity and equivalent driving term.
An effective simulation of the interaction between HF and NF to capture the dynamic evolution of the HF network in a cemented naturally fractured reservoir is particularly important to measure the reservoir stimulation effect. This requires that the numerical model should contain both HF and NF, and assign different rock mechanical parameters for them. Moreover, the numerical solution of the HF propagation path should depend on the fracture criterion, damage evolution criterion, and boundary conditions completely, rather than along a specified path.
In this paper, a method based on an explicit integration scheme with global coupled 0-thickness cohesive elements in rock matrix elements of 3D hexahedra is proposed. Through this method we achieved stochastic propagation of HF and interaction behavior with NF, it also simulated the evolution of complex fracture network evolution for staged fracturing of naturally fractured shale reservoirs, and analyzed the key factors influencing the complexity of the HF network. The structure of this paper is as follows. Section 2 describes the numerical simulation methodology, including the governing equations and cohesive model. Next, we focus on describing the way coupled cohesive elements are implemented in Section 3. The impacts of cementation strength of NFs, injection rate, in situ stress difference, distribution characteristics of NFs, and stress interference under segmented fracturing on fracture network evolution are described in Section 4. A summary and some insights based on the numerical results are given in Section 5.

| SIMULATION METHODOLOGY
In this study, the propagation process of HF and the NF was simulated by coupling fluid seepage in a porous medium with the elastic deformation of the matrix based on the theory of porous medium. The porous medium theory uses effective stress to describe the transfer of stress and is considered as a more reliable model for hydraulic fracturing simulation calculations at this stage.

| The flow-structure interaction model of the porous medium
The rock matrix is subjected to external forces such as intra-fracture fluid pressure and pore pressure during hydraulic fracturing. Assuming that the fracture expansion process is force balanced, the stress balance equation of the current system can be expressed as Equation (1) 29 for any moment t in the domain based on the principle of virtual work.
where σ is the true (Cauchy) stress, δε is the virtual rate of strain, δν is the virtual velocity, t is the traction force per unit area, fˆis the body force per unit volume, dV is the volume element, and dS is the area element.
Define the weight of the wetting liquid as f sn n ρ g = ( + ) , where s is the saturation, n is the porosity, n t is the liquid volume absorbed per unit matrix, (sn + n t ) indicates the total water cut, ρ w is the density of the wetted liquid, and g is the gravitational acceleration. For simplicity, an explicit form is written for the gravitational loading of the wetting liquid, so that the gravitational terms in fˆwill be related to the nonwetting porous medium only. Thus, we can rewrite the virtual work Equation (1) as follows: sn n ρ gδvdV where f are all body forces except the weight of the wetting liquid. We consider that a certain volume of rock (this volume occupies space V with surface S) contains a fixed amount of solid material in the matrix, then the total mass of wetting liquid in this volume of rock is calculated as where V w is the wetting liquid volume flowing through V at any one time, V t is the trapped liquid volume, n w is the wetting liquid volume dV w divided by the total volume dV.
Using the wetting fluid pore pressure as the nodal variable and interpolating on the element, the continuity equations for the porous medium can be obtained in weak form as Equating the change rate of the total mass of fluid in the control volume V to the mass velocity of the fluid through the surface S gives the fluid mass continuity equation 30 as shown below where n is the outward normal direction to the surface S and dt is the time increment step, v w is the seepage velocity.
The differential form of the fluid continuity equation is given by the divergence theorem, as shown below ZHANG ET AL.

| 1645
The pore pressure boundary condition is considered to be p p . Constrain the boundary node displacements in x, y, and z directions. The coupled equation matrix of the stress balance equation and fluid continuity equation is calculated by the direct method so that distribution characteristics of stress, strain, displacement, and other related parameters in the simulated area can be obtained.

| Cohesive model
To avoid the defects of the Griffith fracture theory, Dugdale 31 and Barenblatt 32 first proposed the CZM. Based on the theory of cohesion damage, the propagation process of HF is the result of the separation of materials in the cohesion zone by overcoming the cohesion effect. The fluid flow within a HF and fluid seepage in a porous medium are described by hydrodynamic equations. The cohesion zone existing at the fracture tip is simultaneously subjected to cohesive force and fluid tension force. The stress singularity at the tip is compensated. When the traction force on the cohesive zone reaches the energy damage threshold, the softening stage will be entered. At this time, the overall stiffness weakening coefficient scalar stiffness degradation variable (SDEG) is used to describe the overall scalar damage variable D. 33 Figure 1 illustrates the fracture propagation driven by liquid based on the cohesive element method. When the fluid completely fills the fracture zone, the cohesive elements are completely damaged and the SDEG = 1. 34 If the fluid does not completely fill the fracture zone, the cohesive elements from the tip of the material to the tip of the cohesive elements will be damaged, but the degree of damage gradually decreases and at this time 0 < SDEG < 1.
The pore pressure cohesive element is defined as linear elastic, and the strain is increasing as the stress on the element increases. Before damage occurs on the pore pressure cohesive element, the constitutive relation as shown in Equation (8).
where t is the stress vector of the cohesive element, t n is the normal traction force of cohesive element, t s is the first tangential traction force, t t is the second tangential traction force, which does not exist in the 2D case, ε n , ε s , ε t are the strain components in the corresponding direction, and K is the stiffness matrix of the cohesive element.
In the CZM, fracture initiation and propagation are governed by the traction-separation law. 35 In this paper, the fracture propagation form using a mixed mode of normal and tangential propagation. The mixed-mode damage law is shown in Figure 2. The figure shows the traction force on the vertical axis and the magnitudes of the normal and the shear separations along the two horizontal axes. The unshaded triangles in the two vertical coordinate planes indicate the response under normal and shear deformation, respectively. The shaded is the projection of the damage response (yellow plane) in the normal and the shear direction under the mixed mode. When the normal displacement between the upper and lower surface of the cohesive element is less than the damage initiation displacement of δ m 0 , the material behaves linearly elastic behaves as linear elastic. When the normal stress of the cohesive element reaches the material tensile strength of t n 0 , the cohesive element is F I G U R E 1 Description diagrams of the cohesive zone with fluid drive (cohesive zone demonstration for either hydraulic fracture or natural fracture).
The bilinear cohesive traction-separation law.
damaged, and the damage softening of the cohesive element occurs subsequently. Finally, when the displacement of the upper and lower surface reaches δ m f , the cohesive element will fail completely. For the judgment of fracture initiation, common criteria include the quadratic nominal stress criterion, the quadratic nominal strain criterion, the maximum normal stress criterion, and the maximum normal strain criterion. 8 Through trial and error as well as validation, we found that the simulation results converge stably and maintain high accuracy with the maximum nominal stress criterion. The maximum normal stress criterion can be expressed as follows: where t n 0 , t s 0 , and t t 0 represent the normal stress and the two tangential stresses when the cohesive element is damaged. The symbol  indicates that the cohesive element can only withstand tensile stresses and will not be damaged when it is subjected to compressive stresses. The cohesive element stiffness attenuation is used to describe fracture propagation. The equation is described as Equation (10).
where t ̅ n , t ̅ s , and t ̅ t represent the stress components in the normal and tangential directions of the cohesive element before damage under the current strain conditions.
The damage variable D is calculated by where δ m 0 and δ m f represent the separation displacements at the beginning and complete damage of the element. When the cohesive elements produce a mixed-mode deformation, the effective displacement 36 is introduced to describe the damage evolution occurring in the cohesive elements under the combined action of normal and tangential deformation, as shown below where δ n , δ s , δ t refer to the normal, the first, and the second shear displacement components. The BK criterion 37 is chosen to describe the tensorshear composite fracture propagation, and the first, as well as second tangential energy release rates, are assumed to be equal. In this case, the critical energy release rate G c of the cohesive element in the mixed mode can be calculated by where G n C and G s C refer to the critical energy release rate in the normal, first tangential direction of the cohesive element, G n , G s , and G t refer to the energy release rate in the current normal, first, and second tangential directions, and η is a material parameter. The usual tangential flow behavior of fracturing fluid within a HF includes Newtonian fluid and power-law fluid, and so forth. In this study, it is assumed that the fluid in the fracture is an incompressible Newtonian fluid, 38 and the tangential flow equation can be described as where q is the tangential flow in the fracture, w is the fracture width, μ is the fluid viscosity, and ∇p is the fluid pressure gradient along the fracture propagation direction.
In addition to tangential flow, the fluid inside the fracture can also seep into the rock along the upper and lower surface of the fracture, the equation can be described as follows: where q t and q b refer to the volume flow rate per unit time on the upper and lower surface of the fracture, respectively; c t and c b refer to the leak coefficient on the upper and lower surface of the HF, respectively; p t and p b refer to the pore pressure on the upper and lower surface parts of the fracture, respectively; and p i is the fluid pressure inside the fracture.
The mass conservation equation can be expressed as (16) where Q inj is the injection rate. Substituting Equations (14) and (15) into Equation

| NUMERICAL IMPLEMENTATION
ABAQUS requires the pre-setting cohesive elements to serve as potential propagation paths for HF, which means that the HF can only propagate along the "presetting" path and cannot propagate randomly in the rock matrix. To achieve random propagation of HF and NF, there are three key issues that need to be solved before this numerical model work. First, the pore pressure at adjacent overlapping nodes should be equal. Secondly, all the adjacent cohesive elements should share a single seepage node. Finally, the nodes' order distribution of cohesive elements must ensure that the integration points are in the correct direction for each element. Taking some 3D hexahedral elements as an example, we illustrate the process of embedding cohesive elements globally in detail, and the steps are shown below.
Step 1: Mesh the rock matrix. A structured meshing method is used to create a finite element mesh of the rock matrix, and an 8-node trilinear displacement and pore pressure element (C3D8P type) is selected as the mesh element. As shown in Figure 3A, a C3D8P has eight displacement pore pressure nodes and six faces, respectively. The nodes are numbered in counterclockwise order, and each face has four integration points. When multiple C3D8P elements are united, the common edges of adjacent elements shared the same node. The boundary conditions and different material properties could be defined for different targeted zones within the matrix region where the mesh has been divided. In this way, we provide the basis for describing the heterogeneity of the rock matrix and the random propagation process of fractures.
Step 2: Split the mesh nodes of the rock matrix. Before starting the analysis, we defined the nodes of the matrix elements as the parent nodes. We can see from Figure 3A that there is a face overlap between two adjacent matrix elements, which will result in some edges of the element overlapping two or four times. All we need to do is to build the set of edges that overlap more than two times and record their parent node numbers. Then, we split the parent nodes by the number of times that the common edges are repeated. For example, the parent nodes numbered 2, 5, 14, 17, 9, 12, 10, and 7 are repeated two times, which will create two virtual nodes with different numbers at the same position after splitting one time. The difference is that the parent nodes 8 and 11 are repeated four times. After splitting three times, they will create four virtual nodes with different numbers at the same position.
Step 3: Update the rock matrix and cohesive zone node numbers. After the above two steps, we have inserted the 8-node cohesive element (COH3D8 type) in the mesh of matrix elements. Create a set of cohesive elements directly by editing the INP file. Next, we need to renumber the nodes of the matrix elements and cohesive elements, which is a very critical step. To record the relationship between the split virtual nodes and the parent nodes for python scripting purposes, we specified that the digit of the new number was the maximum digit of the parent node number plus one, and the first digit of the new number is the number of times the parent has split, while the other numbers are the same as the parent. As shown in Figure 3B, the parent nodes retain their original numbering. The embedded COH3D8 cohesive elements are numbered 5, 6, 7, 8 sequentially, and each COH3D8 element shares eight pore pressure nodes. These pore pressure nodes do not need to be specifically defined, ABAQUS automatically gives the nodes calculation capabilities according to the type of element they belong to. The numbering of the new nodes must satisfy that the first edge of the cohesive element formed by the counterclockwise connection has length and the second edge has a thickness of 0. This rule is very important for the correct output of the integration points.
Step 4: Modify the element type and generate the interface. Following the above steps, the default generated COH3D8 elements only support non-fluidsolid coupling cases in fracture mechanics, and do not support numerical simulations of hydraulic fracturing with fluidsolid coupling. Therefore, we must modify COH3D8 cohesive elements to COH3D8P cohesive elements. Figure 3C shows that a COH3D8P cohesive element has one more pore pressure node at the thickness edge compared to a COH3D8 cohesive element, so it shares a total of 12 pore pressure nodes. We write the dictionary set for the newly generated COH3D8P cohesive element, which the key as the element number and the value as the node number. In this way, a new interface is generated between the discrete meshes, and the parent nodes' numbers are critical to identifying the interface between two elements.
Step 5: Merge the intersecting cohesive element nodes. Up to now, we have embedded some COH3D8P cohesive elements in the matrix mesh, but they still cannot work because the pore pressure at the intersection point of the intersecting cohesive elements is unequal to achieving continuous transfer of fluid pressure, which results in the failure of fracture propagation. We judge the spatial location of the intersecting nodes by the coordinates of the nodes and merge the four pore pressure nodes of the cohesive | 1649 element located at the thickness edge under the same coordinates to make special nodes with pore pressure freedom, thus ensuring the continuity of the stress and fluid fields at the intersecting nodes. As shown in Figure 3C, in fact, the location of the parent node numbered 8 contains a total of one parent node and seven new nodes. As can be seen from Figure 3D, when the nodes were merged, the same position contains a total of one parent node and four new nodes with nodes numbered 8, 108, 208, 308, and 321.
Finally, the rock matrix elements are connected to their adjacent cohesive elements by sharing two common nodes, and the two adjacent cohesive elements are connected to each other by sharing one node. The cohesive elements will experience loading, damage, stiffness degradation, and fracturing under internal/external traction loads.

| NUMERICAL SIMULATION
The fracturing sequence of a single well was changed from "commuter-frac" fracturing to "Texas two-step" fracturing in the early stage to enhance stress interference between multiple fractures so that the fracture complexity could be increased. 1,39 To study the evolution process of this complex fracture network, many scholars have numerically analyzed various factors impacting the intersection of NF and HF. [40][41][42][43] In this paper, we developed a new improved 3D HF model with a global embedded cohesive zone and analyzed the impacts of cementation strength of NFs, injection rate, in situ stress difference, distribution characteristics of NFs, and stress interference on the HF network.

| Model setup
Detournay 44 believes that two types of energy dissipation are present in the HF propagation process, viscosity-dominated and toughness-dominated. Viscosity-dominated means that hydraulic energy is consumed primarily in the form of fracture fluid flow friction, and toughness-dominated means that hydraulic energy is consumed in the form of overcoming rock fracture energy. Studies 45 have shown that the HF propagation process is mainly viscosity-dominated even when slickwater is used as a fracturing fluid. Table 1  Based on the method described in Section 3, we wrote the ABAQUS plug-in by using Python language to achieve bulk embedding of cohesive elements, and the final model was composed of rock matrix elements and cohesive elements. Since there are a large number of natural cracks in the model, the use of cohesive elements to represent the initiation and propagation of fracture may cause the computation not to converge easily and affect the calculation accuracy. For this reason, the UEL subroutine was written to improve the convergence of the embedded cohesive elements in the calculation process. Figure 4 shows the schematic diagram of the plug-in writing process.
To verify the accuracy of the model, we build a model with dimensions of 100 m × 60 m × 1 m, which makes the simulation domain much larger than the opening and length of fracture ( Figure 5). The solid element type of the model is C3D8P type, and COH2D4P elements are embedded globally. σ H , σ h , and σ z are applied as the maximum principal stress, minimum principal stress, and vertical stress in the x, y, and z directions of the model, respectively. Figure 6 shows the analytical solution 46,47 compared with the simulation results. The results show that numerical results agree well with the analytical solution except for the net pressure at the beginning of the injection stage. This is due to the small cohesive zone in front of the fracture tip in the numerical model, which causes the rock to show linear elasticity, however, the analytical solution does not consider the pressure accumulation period before the fracture initiation. Figure 7 shows the schematic diagram of this numerical simulation work, it can be seen from the figure that the size of the model is 100 m × 60 m × 1 m, the edge length of meshes is 1 m, the displacement on the boundaries of model x, y, and z directions are constrained, and the pore pressure of the boundaries are set to 20 MPa. The horizontal wellbore direction is defined as the direction of the minimum principal stress (x-direction). The mesh of the numerical model is composed of three parts: the rock matrix elements, the  | 1651 cohesive elements that support random propagation of HF, and the cohesive elements that represent NFs. The parameters of these three parts are listed in Table 1.

| Numerical results
Numerous studies have shown the presence of NFs in shale reservoirs like Marcellus Shale 48 and Barnett Shale. 49 The potential role of NFs has been confirmed by production analyses. For the case of Barnett Shale, NFs are mostly filled with calcite and quartz from the underlying Ellenburger limestone formation, but the NFs that are sealed still have weak cementation strength. 50 Cemented fracture zones are often acted as preferential weak paths for inducing fracture deflection to enhance the stimulation effect of hydraulic fracturing, expecting their interaction with HFs to form a more complex fracture network. In fact, hydraulic fracturing in the field has to face more complex geological environments and engineering conditions. Therefore, this subsection focuses on the impacts of geological and engineering factors on the HF network evolution in a cemented naturally fractured reservoir. A series of numerical modeling schemes were designed by using orthogonal tests. The benchmark scheme has the smallest input parameters, and the comparison schemes have sequentially increasing parameters. In this paper, we keep the values of other parameters as standard values when adjusting the target parameter to study the impacts of sensitive parameters on fracture propagation. Figure 8 shows the propagation of HF as the cementation strength of NFs (C NF ) varies from 1 to 4 MPa. When the cementation strength of the NFs is 1 MPa, the HF can completely activate the NFs near the injection point and use the tip of the NF as the new initiation location ( Figure 8A). Because the properties of the rock matrix and the NFs differ significantly, the NF tip will gather a high fluid pressure before being breached, which will cause other adjacent NFs to be induced to be destroyed but not activated by the injected fluid. We observed two damage events, one is that the pore pressure at the contact surface between the NF and the rock matrix will be increased because of the percolation of fluid in the porous medium, and the side close to the rock matrix has different degrees of damage. The other one is the stress interference in the branch fractures close to each other after the HF activated NFs, and the superposed stress field caused different degrees of damage to the rock matrix. Increasing the cementation strength of the NFs to 1.5 MPa, the tip of the NF is no longer used as the initiation location after the HF activates the NF, and the new initiation location tends to be shifted toward the interaction point ( Figure 8B). This indicates that the degree of NF cementation has an influence on the new initiation point of the intersecting HF, and the tip of the NF with stronger cementation no longer serves as the new initiation location of the HF, but gradually shifts to the interaction point. In this situation, the branching HFs are close to each other and stronger stress disturbances may lead to unequal lengths of branching fractures propagation. When the cementation strength of NFs is 2 MPa, HF tends to activate NFs while breaking through the interaction point, and after that, the length of other activated NFs in the propagation path of HF is reduced ( Figure 8C). The HF will not be able to activate the NFs completely when the cementation strength of the NFs continues to increase ( Figure 8D). At this time, the propagation length of HF is increased greatly. When the cementation strength of NFs is 4 MPa, the HF will lose the ability to activate the NFs, and finally, propagate in the form of a conventional symmetrical double-wing fracture ( Figure 8E). As shown in Figure 9A, the cementation strength of NFs has a greater effect on the opening width of HFs, each intersection between HF and NFs will reduce the opening width of HF, and this diameter shrinkage phenomenon is improved gradually as the cementation strength of natural cracks increases. The fluid within the HF needs to overcome the fracture energy of the NFs to activate the NFs. As we can see in Figure 9B, we can see that improving the cementation strength of NFs can reduce the interaction behavior of HFs with NFs and increase the propagation distance of HFs. In addition, the propagation pressure of HF decreases significantly as the cementation strength of NF increases. This is because a fractured reservoir presents strong heterogeneity when the cementation strength of NFs is low. As the cementation strength of NFs increases, the reservoir inhomogeneity is improved and the pressure required for HF propagation decreases.

| Impact of the injection rate
Injection rate (Q inj ) is an important controllable parameter of stimulation treatment. Increasing the injection rate will facilitate the complexity of the HF network. 51,52 In this section, four cases of injection rate, 0.010, 0.014, 0.018, and 0.022 m 3 /s are investigated. Comparing Figures 8A and 10, it can be seen that increasing the injection rate can improve the complexity of HF networks and increase the SRV. We observed that increasing the injection rate does not seem to change the propagation behavior of the HF after interacting with the NF. After HFs activated the NF, the location of new initiation still tends to the tip of the NF (Figure 10A-D). When the injection rate is 0.018 m 3 /s, the HF even completely followed the NFs propagation in the horizontal direction and no longer propagated along the direction of the maximum principal stress ( Figure 10C). When the injection rate is 0.022 m 3 /s, the HF not only activates the NFs around the path but also increases the propagation length to 80 m ( Figure 10D). This indicates that the injection rate dominates the complexity of the HF network with other conditions remaining constant.
In addition, increasing the injection rate may cause more NFs to be shear-damaged. These shear-damaged NFs have formed a certain width geometrically, but these NFs and HFs are still separated by the rock matrix. The lower rock matrix permeability does not allow fluid flow to these induced open NFs, so they do not have the ability to increase the well production either. As shown in Figure 11A, when the injection rate is 0.010 and 0.014 m 3 /s, the HF propagation length and opening width are increased, which is beneficial to expand the reservoir drainage area. Continuing to increase the injection rate to 0.18 m 3 /s, the HF activated the NF in the horizontal direction continuously, which caused the propagation length to decrease although the opening width of the HF increased. At this time, the main component of the fracture network is NFs. As we can see in Figure 11B, we can see that the initiation pressure of HF increases from 32.4 to 34.1 MPa as the injection rate increases. Moreover, there is an increased chance of HF interacting with NFs as the injection rate increases. When a HF intersects with a NF, the HF breaks through the tip of the NF to produce a branch fracture, and there is a significant rise in injection point pressure at this time. This is because changing the direction of HF propagation requires greater injection pressure support. The fluctuation of the injection pressure can be clearly seen in the interaction zone (the area within the dashed box in the figure) from Figure 11B. In this zone, each lift of injection pressure indicates that the HF has interacted with the NF and created a branching fracture. When the injection rate increases to 0.022 m 3 /s, the diameter shrinkage phenomenon after the intersection of HFs and NFs decreases significantly, which makes the injection pressure fluctuation in the interaction zone become flat. Figure 12 illustrates the case of HF propagation with minimum principal stress of 8 MPa and maximum principal stresses of 10, 12, 14, and 16 MPa, respectively. When the in situ stress difference is 2 MPa (Figure 12A), the HF tends to propagate along the NF tip direction. The NFs with lower fracture energy at low in-situ stress difference provide a dominant channel for the propagation of HF and lead the propagation direction of HF. In this situation, the two branching fractures after the HF activation of the NF no longer seem to propagate simultaneously, and one of the branching fractures has a limited propagation length, which will diminish the complexity of the HF network. Increasing the in-situ stress difference to 4 MPa ( Figure 12B), the HF is able to activate the NFs completely and then propagate along the maximum principal stress direction after crossing the interaction point. When the in-situ stress difference continues to increase ( Figure 12C), the HF can only completely activate the NF near the injection point, so the NFs far from the injection point are activated less. When the in-situ stress difference is 8 MPa (Figure 12D), the NFs with weak cementation strength and near the injection point can still be partially activated by the HF, but the NFs which are far away from the wellbore will not be affected by the HF. This is because the higher level of in situ stress difference increases the normal stress acting on the NF surface, thus increasing the resistance to opening the NFs. The net pressure on the fracture interface will have to rise to exceed the stress anisotropy, and the number of activated NFs in the fracture network system is regulated by the stress anisotropy. Therefore, when the in situ stress difference is low, NFs as weak surfaces are still the main factor controlling the formation of a complex fracture network. After the in-situ stress difference is higher than 4 MPa, the number of activated NFs decreases significantly, and the stress dominates the propagation direction of HFs at this time.

| Impact of the in situ stress difference
In this case, HF steering and branching are reduced, which may result in a less complex fracture network. As can be seen from Figure 13A, the width of HFs decreases with increasing in-situ stress difference, but the opening width increases significantly. Moreover, Moreover, the increase of the in-situ stress difference can significantly improve the diameter shrinkage of the HF after activating the NFs. Figure 13B shows that increasing the in situ stress difference is unlikely to increase the injection pressure. This is because the HF initiation is mainly to overcome the sum of the minimum horizontal principal stress and the rock tensile strength. In this simulation, the minimum horizontal principal stress is a constant value, so the variation of HF initiation pressure is small. Meanwhile, with the increase of the insitu stress difference, the opportunity for interaction between HF and NF decreases. The propagation pressure of the HF decreases from 24.8 to 20.9 MPa, and the curve shape becomes smooth. In this condition, the HF propagates along the maximum principal stress and is less influenced by the cemented NF.

| Impact of the spatial distribution of NFs
Simulate the evolution of HF networks in a cemented naturally fractured reservoir with the random orthogonal characteristic distribution by increasing the number of NFs in the vertical wellbore direction. It can be seen from Figure 14 that the distribution characteristics of NFs in the reservoir significantly affect the morphology of the HF network. The variety of NF distribution can change the path of HF propagation, and make the HF potentially deviate from the direction of maximum principal stress. As the number of NFs increases, the distance of HFs "slipping" in the NFs increases significantly, which is beneficial for the HF to activate the NFs far away from the injection point as far as possible under the lower injection pressure and results in a higher percentage of NF length in the total length of the fracture network. In addition, the distribution of NFs with different angles helps to offset the loss of hydraulic pressure head caused by the HF activating the NF of high-angle. After a HF opens the NF of a high angle, the possibility of continuing to activate the NFs of low angle increases, which will enhance the complexity of fracture networks. However, we need to be aware that reservoir heterogeneity is enhanced as the number of NFs increases. As a result, the HF may activate the NFs distributed near the injection point, and the NFs far from the injection point may not be activated. The results of this simulation illustrate that when other conditions hold unchanged, the distribution of cemented NFs dominates the propagation direction of the HF and determines the complexity of the fracture network. From Figure 15A, we can see that increasing the number of orthogonal NFs at the same injection rate causes the opening width of HFs to decrease but the length to increase. We observe that more NFs are induced to open. This is because an increase in NFs will enhance the reservoir heterogeneity, which results in higher concentrated stress at the tip of the HF during propagation, and the nearby NFs have an increased possibility to be induced to open. Besides, the HF activating more low-angle NFs is beneficial to improve the diameter shrinkage phenomenon of HFs. Figure 15B shows that the NFs with several angles in the reservoir do not change the fracture initiation pressure of the HF, but have a greater impact on the propagation pressure of the HF.

| Impact of the multiple HFs
We increased the number of perforations to three as a way to simulate synchronous hydraulic fracturing, commuter hydraulic fracturing, and Texas two-step hydraulic fracturing. As can be seen from Figure 16, the stress interference effect among the HFs enhances the complexity of fracture networks as the number of perforations is increased. As a result of the stress interference effect, the interaction pattern of HFs and NFs becomes complicated, the event of NFs being activated is difficult to predict, and the HFs show asymmetry in the propagation process. Figure 16A shows that during synchronous hydraulic fracturing, the three fractures are interfered with by the stress of each other at the same time. Moreover, the adjacent parts of the intersection are open to forming a "cross"-type fracture network, which will provide a high-speed pathway for the transport of oil and gas. Figure 16B illustrates that the HFs in the later initiation step are likely to be in the stress shadow area of the earlier initiated HFs. The earlier initiated HFs tend to form an L-shaped fracture network after activating the NFs. The continuous propagation and reorientation to the maximum principal stress direction of the subsequent HFs may lead to forming of an S-shaped fracture, and restrict the growth of the HF length. As shown in Figure 16C, the Texas two-step hydraulic fracturing can ensure enough stress interference among multiple fractures to form a complex U-shaped fracture network, permitting tighter shot spacing without the concern of generating so much induced stress that the intermediate fracture propagation is limited. The numerical results also support the understanding in the fracturing community that the Texas two-step hydraulic fracturing is the preferred method for hydraulic fracturing in a naturally fractured reservoir.

| CONCLUSIONS
In this paper, a 3D, fully coupled hydraulic fracturing model has been established using the 3D hexahedral elements to describe the rock matrix and globally embedded 0-thickness cohesive elements. In the proposed model, the propagation of fractures and fluid flow within the HFs are governed by the CZM. The coupled equations are discretized and solved by using an explicit temporal integration method, and a subroutine is used to improve convergence during the simulation of the model. Finally, the model was used to investigate the geological and engineering factors that influence the evolution of the HF network in a cemented naturally fractured reservoir. The following conclusions can be obtained from the simulation results.
1. The cemented NFs increase the tortuosity of the HF propagation path, and each interaction of HFs with NFs caused the diameter shrinkage of the HFs. When the cementation strength of NFs is small, the HFs tend to offset toward the tip of NFs and then form a complex HF network. However, when the cementation strength of NFs is larger, the HFs prefer to pass through the NFs directly from the interaction point. In this situation, the propagation path of HFs is more regular, but the complex fracture network is difficult to form. 2. Increasing the in situ stress difference reduces the possibility of NFs being activated, the path of HF propagation becomes longer, and a fracture network is difficult to form. In this situation, the dominant factor in the morphology of the HF network is stress. The change in injection rate and the number of NFs do not modify the interaction behavior of HFs and NFs, which work together with the cementation strength of NFs to determine the morphology of the HF network. 3. Changes in the in situ stress difference, the cementation strength of natural cracks, and the number of NFs do not seem to influence the HF initiation pressure but have a significant impact on the pressure during the propagation of the HF. The interaction of a HF with NFs will require higher injection pressure than HF propagation alone. This is because hydraulic fracturing needs a higher fluid pressure to overcome the energy loss caused by the propagation direction shifting after reorientation. 4. Increasing the injection rate is beneficial to improve the possibility of interaction between HFs and NFs, and cause a significant increase in injection pressure. The HF activating the NF is more likely to branch and intersect. The HF appears to create more zones of inactivated damage that can be recorded by microseismic monitoring solution when interacting with NFs as the injection rate increases, which may cause an incorrect assessment of the SRV. 5. Although the pre-existing NFs create a favorable geological condition for segmented fracturing to form a HF network, the stress interference caused by multiple fractures may alter this advantage. This is because some HFs may form S-shaped fractures with short lengths after reorientation and propagation under the impact of stress interference, which is not conducive to increasing the reservoir drainage area.