Three‐dimensional complex fracture propagation simulation: Implications for rapid decline of production capacity

In this study, a simulation model of fracture network geometry in fractured and porous elastic reservoirs is established based on the globally embedded three‐dimensional cohesive zone model (3D CZM). The effects of stress interference between natural fractures (NFs), horizontal stress difference (HSD), injection rate, and fluid viscosity on fracture geometry were studied and the phenomenon of rapid failure of unconventional reservoir productivity was explained from a coupling fluid flow/geomechanics perspective. The intersection behavior of hydraulic fracture (HF) and NF is simulated; the simulation results are in good agreement with Blanton's criterion, proving the reliability of the 3D CZM. The research results show that in the process of hydraulic fracturing, in addition to the strong stress interference between HFs, there is also strong stress interference between NFs. Stress interference has an important impact on the traction separation of cohesive elements (CEs) and a continuous dynamic impact on the NF opening. For reservoirs with a small HSD, higher injection rate and fluid viscosity produce smaller stimulated reservoir volume (SRV) length, but larger SRV width. Complex fracture networks tend to develop in the vicinity of horizontal wellbores, which leads to rapid failure of productivity. For this kind of reservoir, minimizing the impact of stress interference with a small injection rate is recommended instead of maximizing the injection rate and fluid viscosity. For reservoirs with a large HSD, lower injection rate and fluid viscosity produce larger SRV length and smaller SRV width. The fracture geometries tend to develop simple straight fractures, which leads to low initial productivity. Thus, adopting a relatively large injection rate and fluid viscosity is recommended to make full use of stress interference effects to increase the development of complex fractures near horizontal wellbores. The results of this study can serve as a guide in evaluating fracture complexity, SRV, proppant migration, drainage reservoir volume (DRV), and well completion design; they can also promote the in‐depth understanding of the coupling effect between fluid flow and geomechanics.


| INTRODUCTION
As a mature stimulation measure in the oil industry, hydraulic fracturing has become the leading method to stabilize oil field production, from the earliest conventional vertical well fracturing to the horizontal well multi-cluster fracturing technology in unconventional reservoirs. 1,2 In the development of shale, tight oil, and tight gas, and in the refracturing of old wells, hydraulic fracturing technology is generally adopted for reservoir reconstruction. For unconventional reservoirs, especially low-permeability fractured reservoirs, under certain conditions, hydraulic fracturing can take advantage of developed NFs and the coupling effect of fluid flow and geological stress to form a highly permeable fracture network. [3][4][5][6] Microseismic field data confirmed that some reservoir fracture zones can be formed after hydraulic fracturing. 7 The corresponding capacity analysis showed that the larger the SRV, the larger the initial capacity. However, field production statistics also show that high displacement fracturing wells can achieve high initial productivity followed by a rapid decline in productivity. 8 At least 50% of new wells in China's Changing Shale Development Demonstration Zone face rapid decline. 9,10 In hydraulic fracturing, the fracture extension exhibits a complex behavior, affected by the injection rate of fluid, fluid properties such as viscosity, rock parameters (such as the elastic modulus of the rock, tensile strength, and shear strength), NF parameters (such as NF approaching angle, NF location, length, tensile strength, and shear strength of NF), and the coupling effect between fluid flow and geomechanics. [11][12][13] At present, the main research methods for fracture networks include the discontinuous displacement method (DDM), finite element method (FEM), extended finite element method (XFEM), discrete element method (DEM), phase field method (PFM), and CZM. [14][15][16][17][18] The 3D CZM can be used to simulate the initiation and extension of three-dimensional HFs, and the intersection, crossing, blocking, and capturing behaviors between HFs and NFs. The 3D CZM is a direct extension of the two-dimensional CZM, as shown in Figure 1. Based on the two-dimensional CZM, zero-thickness cohesive units are embedded in the rock to form the potential paths of three-dimensional HF and NF propagation. The mechanical parameters of the NF CE set and the matrix CE set are defined to have corresponding mechanical properties. Based on the damage criterion, when failure occurs, the CE loses its stiffness due to traction separation. The behaviors of HF and NF, such as intersection, crossing, extension, turning, and capture, occur in the CEs. Compared with DDM, XFEM, and DEM, 3D CZM eliminates the singularity at the fracture tip and overcomes the disadvantage of linear elastic fracture mechanics that can diverge, properties that are at the core of its ability to simulate three-dimensional fracture network expansion. [19][20][21] Barenblatt 22 first proposed the concept of a cohesive zone simulating crack propagation in brittle materials. Yao 23 established a 3D hydraulic fracture propagation model using the pore cohesive zone model, simulated hydraulic fractures with K E Y W O R D S fluid flow/geomechanics, fracture network geometry, horizontal stress difference (HSD), hydraulic fracturing, three-dimensional cohesive zone model (3D CZM) F I G U R E 1 Schematic of 3D CZM different rock characteristics, discussed the difference between linear elastic fracture mechanics and cohesive fracture mechanics, and proposed an effective fracture stiffness method. Carrier and Granet 24 developed a two-dimensional hydraulic fracturing model to simulate the process of hydraulic fracture propagation using the zero-thickness cohesive element and the conventional finite element method. The simulation results were compared with analytical asymptotic solutions under a zero-fluid lag assumption in the four limiting propagation regimes, and the advantage of the cohesive zone model was proved. Mokryakov 25 and Wang et al 26 simulated the fracturing process in soft rock and the hydraulic fracturing expansion process, respectively, under the plastic influence of the fracture tip and rock interior using the 2D CZM. Guo et al 27 established the HF and NF intersection model based on the CZM and discussed the influence of HSD and approaching angle on the intersection of HF and NF. The results showed that a smaller approaching angle and HSD increased the likelihood of NF initiation and expansion. Guo et al 28 then used this model to simulate the propagation behavior of HFs in layered shale. The study showed that HFs mainly penetrated the interface and formed branched fractures along the interface. Wang 29 established a fully coupled nonplanar hydraulic fracture extension model based on the EFEM and CZM and studied the important influence of inelastic deformation caused by hydraulic fracture on the expansion pressure and fracture geometry. The ductile formation was simply regarded as soft rock with a low Young's modulus; this assessment caused great error in the simulation results. Dahi-Taleghani et al 30 simulated the interaction between HFs and preset NFs by using the CZM. To ensure the accuracy of CE parameters, they first obtained the cohesive properties of natural fractures and traction-separation law (TSL) characteristics of rock based on experimental methods. The results showed that the intersection angle between fractures and the nature of cementitious materials in natural fractures may produce an initial direction of hydraulic fracture flow that is perpendicular to the minimum horizontal stress. Wang 10 also established a model to simulate the propagation of complex fractures based on the CZM. The research results showed that only complex fractures, rather than a complex fracture network, are formed in the reservoir after fracturing, due to the influence of natural fracture parameters (including the density, size, direction, and location of natural fractures), stress distribution, and hydraulic fracturing construction parameters. Yu et al 31 simulated the propagation of fractures in the rock matrix and the network of intersecting natural fractures based on a coupled fluid flow and deformation finite element model with adaptive insertion of three-dimensional cohesive elements. The effects of different pumping programs on the complexity of fracture networks were studied. Li et al 32 simulated the interaction between HFs and NFs based on the two-dimensional modified pore pressure CZM and analyzed the impact of temporary fracturing on fracture complexity. Liu et al 33 simulated the influence of vertical stress on fracture extension in the interaction between HFs and NFs using the 3D CZM. The results showed that in addition to HSD, cluster spacing, and construction parameters, vertical stress also has a significant impact on the direction of fracture extension.
In this study, the reservoir is regarded as a porous elastic rock formation with natural fractures. The intersecting initiation and extension models of HFs and NFs are established by the 3D CZM. The important influence of coupled fluid flow and geomechanics on the formation of a complex fracture network is discussed, and the influence of changing horizontal stress difference, injection rate, and fluid viscosity on fracture geometry is analyzed. Our main purpose is to study the complex fracture morphology under different conditions, analyze the potential production capacity of the fracture morphology, and provide corresponding construction suggestions to explain the rapid failure of unconventional reservoir production capacity. Our research offers insight into the formation mechanism of underground fracture networks, and guidance in evaluating fracture complexity, SRV, drainage reservoir volume (DRV), and well completion design.

| MATHEMATICAL MODEL
The mathematical model is divided into three parts: the fracture propagation model within the three-dimensional cohesive elements, the traction-separation criteria used to judge the initiation and propagation of fractures, and the coupled model of fluid flow/geomechanics. In the fracturing process, temperature has little influence on fracture propagation in a few hours, so the effect of temperature on fracture is ignored in our model.

| Fracture propagation model
For HF and NF, initiation and propagation can be regarded as the failure of CE and loss of stiffness. The 3D CZM embedded in the porous matrix contains all potential extension paths after initiation of HF and NF. Therefore, all initiation and propagation behaviors occur in the cohesive element zone. The driving force of these behaviors is the injected viscous fracturing fluid. In the cohesive element zone, the physical behavior of fracturing fluid flow can be controlled by the following four equations.

| Material balance equation
According to the principle of material balance, the material balance equation for fractures is 23,24 (1) where q f is the flow rate within the fracture, w is the fracture opening, q l is the filtration rate of the fracturing fluid, Q (t) is the injection rate, and (x, y) is the Dirac delta function.

| Flow equation
Reservoirs with NFs, such as shale, are generally hydraulically fractured with slippery water. The fracturing fluid has Newtonian fluid characteristics and is incompressible; the pressure drop equation in an HF can be expressed as 34 where p f is the fluid pressure in the fracture in MPa, and w is the liquid viscosity in mPa s.

| Fracture width equation
The stress interference between HFs and NFs, and between NFs continuously affects crack expansion and the opening degree. Using FEM to calculate the stress distribution and damage coupling of CE, the crack opening is determined by the displacement of the top and bottom of the fracture wall, and is the result of coupling CE characteristics, fluid filtration, pore pressure, fluid pressure, and stress distribution. [35][36][37] where n represents the unit normal vector on the hydrofracture surface (dimensionless), and u t and u b represent the displacement of the top and bottom wall surfaces of the fracture, respectively, as shown in Figure 1D.

| Fluid filtration
For a complete fluid-solid coupled model, Darcy's law is used to establish a coupling fluid filtration model with leakoff coefficient, fluid pressure, and pore pressure as the variable terms. 38,39 where c t and c b represent the fracturing fluid leak-off coefficients at the top and bottom crack wall, respectively, in m 2 / (s Pa), q t and q b denote the filtration rates of the fracturing fluid on the top and bottom element surfaces, respectively, and p t and p b represent the pore pressures on the top and bottom element surfaces, respectively, in Pa.

| Traction-separation law
During the formation of a fracture network, the judgment criteria of initial damage, initiation, and expansion of fracture are critical. In this model, the traction-separation law is based on the 3D CZM, which defines the constitutive relationship of the bonding performance of the cohesive element at the crack tip. In complex three-dimensional fracture propagation, there are many mixed extension modes due to interactions between HFs and NFs. Therefore, the initiation and extension of the crack are judged using the quadratic nominal stress criterion. [40][41][42] In this criterion, the effects of normal stress and shear stress are considered and determined by three nominal stress components and three nominal stress peaks. This criterion can be expressed as where t n , t s , and t t represent the normal, the first, and the second shear stress components, respectively, t 0 n , t 0 s , and t 0 t refer to the normal and two shear strength of intact cohesive surfaces, respectively. The ⟨⟩ symbols represent Macaulay brackets, indicating that the pure extrusion deformation or stress state does not cause damage. Equation (6) indicates that damage begins when the sum of squares of the total nominal stress ratio reaches 1. The nominal stress component and the elastic parameters are determined by the linear elastic constitutive relation and the material properties of the CE. In addition, the stress component of the traction-separation model is also affected by the damage variable D during the failure of the CE according to where t = t n , t s , t t refers to the nominal stress vector; t = t n , t s , t t represents the stress components predicted by the elastic traction-separation behavior without damage for the current displacement. In this model, a bilinear cohesive law 43,44 is used to describe the relationship between traction force and displacement, as shown in Figure 2. In Equation (7), dimensionless D represents the overall damage of the material without damage cohesion unit during fracturing. D varies with the fracturing time and increases linearly from 0 to 1 after the initiation of damage. In the model, the overall damage variable D is determined by the traction-separation law and is calculated as where 0 m and f m denote the displacements when the traction reaches T max and the fracture is completely damaged. 0 m is a constant related to the properties of the material that can be measured experimentally. m is the effective displacement, defined as where n , s , and t represent the normal, the first, and the second shear displacement components of a CE, respectively. In Equation (8), f m can be defined as where T represents the tensile or shear strength and is a constant calculated form the parameters of the CE material and G c refers to the mixed-mode fracture energy. Once a CE breaks, the damage can be evaluated according to fracture energy theory. The damage evolution of a CE during fracture propagation is determined by the Benzeggagh-Kenane fracture criterion, 30 defined as where G n , G s , and G t represent the work done by the traction and its conjugate relative displacement in the normal, the first, and the second shear directions, respectively. These parameters can be obtained using a traction-separation curve. To characterize the material properties more accurately, the tractionseparation curve is obtained by experimental method. 30 Here, , G c s , and G c t represent the fracture energies in the normal, the first, and the second directions, respectively. 45

| Continuity equation
It is assumed that the control volume of the unit is V, the surface area is S, and the volume in the initial state is V 0 . Free fluid can enter the unit through the surface. If the fluid volume in the control body is V w , the total mass of the fluid in the control body is expressed as where w is the fluid density and n w = dV w .
The change rate of fluid mass per unit time is equal to the mass of the fluid passing through the control body and surface such that 47,48 where J = dV 0 , v w represents the fluid flow velocity, and n denotes the unit normal vector.
By differentiating Equation (12), we obtain where x is coordinate component.

| Flow equation
During fracturing, the fracturing fluid has two main roles: forming fractures through hydraulic driving, and entering the reservoir matrix through the fracture wall surface to cause the pore elastic effect. The flow behavior of the filtrated fracturing fluid is characterized by Darcy's law, where k is the permeability vector, g is the acceleration due to gravity, and ∇ represents the pore pressure gradient.

| Stress balance equation
According to the virtual work principle, the stress balance equation for the elastic model of porous media is 49 where is the imaginary strain, represents the Cauchy stress, v denotes the imaginary displacement, f S is the force per unit area, f and sn w w together form the volume force per unit volume.

| Dispersion of coupling equations
Equation (16) is discretized by the polynomial difference method. The form function N N (representing the displacement form of the element) and N (reflecting the relationship between the virtual displacement and the virtual strain of the node of the element) are introduced. The system of equations after discretization becomes The weak form of the fluid seepage stress equation is where u w represents the variation of the pore fluid, J denotes the change rate of the volume of the porous medium, and v represents the seepage velocity after the equation becomes a weak form.
Through the divergence theorem, Equation (20) becomes where 0 w represents the fluid density in the reference configuration. The Newton iterative method is used to solve Equations (19) and (21).

| MODEL VERIFICATION
In this section, the KGD model and Blanton's criterion are used to validate the proposed global embedded 3D CZM. 36,50 The curves of the net pressure, injection point width, and the injection time obtained based on the KGD model and 3D CZM are shown in Figure 5A. It can be observed that the results of our model are basically consistent with those of the KGD model. During the formation of a fracture network, the intersection criterion determines the extension mode after the intersection of HF and NF, which is the key to studying the fracture network geometry. In this section, we use Blanton's criterion for the intersection of HFs and NFs to verify the proposed globally embedded 3D CZM. Blanton's criterion is based on the stress elastic solution in the interactive region, and it is assumed that the HF is passivated when it meets the NF. According to the criterion, when the pressure in the joint at the interaction is greater than the normal stress on the wall of the NF, the NF opens. The HF will pass through the NF when the fracture initiation pressure is lower than the pressure required for the fracture to open. Under different HSDs and different approaching angles, the extension mode after the intersection of HF and NF is shown as crossing and capturing. Based on the 3D CZM, a 3D intersection model of HF and NF with different HSDs and different intersection angles is established, as shown in Figure 3A. Figure 3B shows the connection of CEs at the intersection point. The red dots represent the fluid pressure nodes in the cohesive unit. At the intersection point of the four CEs (CE1, CE2, CE3, and CE4), the fluid pressure nodes have the same position coordinates and share one fluid pressure node. In the modeling process, to ensure accuracy and save calculation cost, the grid was encrypted near the prefabricated CE and was sparse further from the CE. The main calculation parameters in the model are shown in Table 1. Figure 4 shows the HF and NF intersection results with HSD = 3 MPa and intersection angles of approximately 30° and 60°. In the model, the minimum horizontal principal stress is along the x-axis and the maximum horizontal principal stress is along the y-axis. The simulation results show that when the HSD is 3 MPa and the approaching angle is 30°, the HF will not pass through the NF, but will be blocked and captured by the NF, changing the original propagation direction of the fracture. When the HSD is 3 MPa and the approaching angle is 60°, the HF passes through the NF because the normal stress on the natural crack wall increases with the increase of the approaching angle, making it more difficult for the natural crack to open. To fully demonstrate the correctness of the 3D CZM, we conducted multiple sets of numerical simulations. The comparison between the simulation results and Blanton's criterion is shown in Figure 5B. Based on Blanton's criterion, the area under the curve indicates that with an HSD and approaching angle in this range, HFs are captured by NFs and expand along the direction of NFs. In the upper part of the curve, the HF passes through the NF. The position on the curve indicates that under this condition, HF may pass through or be trapped, depending mainly on stress and the grid. The simulation results are consistent with Blanton's criterion, as shown in Figure 5B.

AND ANALYSIS
Field production data analysis has shown that the productivity of many shale gas wells decreases rapidly after fracturing. To determine the cause and explain the internal mechanism from a coupling fluid flow/geomechanics perspective, the influence of the stress interference effect between HFs and NFs and the stress shadow effect on the fracture network geometry is discussed. The geometry of the fracture network with different HSDs and different injection displacements was studied, and the influence of fracture network geometry on later productivity was analyzed. Corresponding suggestions are presented for different geological conditions.

| Stress interference between NFs
Many studies in the literature have shown that when multicluster fractures are simultaneously fractured in homogeneous reservoirs, stress interference such as repulsion, attraction, and extrusion between fractures exists between HFs. In a reservoir with NFs, the question remains if there is stress interference between NFs. The answer lies in proving the effect of fluid flow coupling and geomechanics on the opening, expansion, and geometry of NFs. Thus, we established two physical models to simulate the stress interference between NFs, as shown in Figure 6. The length, width, and height of the two models are 40 m, 20 m, and 10 m, respectively. In Figure 6, Case I, NF-a and NF-b are pre-embedded according to the 3D CZM, and the injection point is on NF-a. Both NF-a and NF-b have the same characteristics. NF-c is embedded in Case II to verify the effect of NF-c on NF-b, with the same fracture position as in Case I. Parameters used in the simulation are shown in Table 2.
The simulation results of the two cases are shown in Figure 7; the crack opening is magnified 100 times. Figure 8 shows the three-dimensional fracture of the two cases at the end of fracturing using MATLAB program software. opening is smaller at the root (closer to NF-a) and larger in the middle (away from NF-a). As NF-a continues to expand, the NF-b opening in region C, which is inhibited by NF-a, continues to decrease but does not close (the opening does not reach zero, as shown in Figure 9). In Case II, the CEs distributed on the right side of NF-b and NF-c exhibit traction-separation behavior (regions D and E). The CEs on the left side of NF-b and NF-c were not damaged by the extrusion of NF-a. During the traction separation of NF-b and NF-c, the change in opening size is similar to that of NF-b in Case I. Inhibited by NF-a, the openings of NF-b and NF-c in region F gradually decrease but do not close (the opening does not reach zero, as shown in Figures 8 and 9). To more clearly reflect the influence of NF-c on NF-b, the opening changes of four CEs on NF-b in Cases I and II (Figure 8) were recorded for comparison, as shown in Figure 9. As observed in Figure 9, the opening of the CE increases first and then decreases with the increase of injection time. The opening of the CE on NF-b in Case I is always greater than in Case II, proving that there is stress interference between NF-c and NF-b. In addition, the results of both cases show that the closer to NF-a, the smaller the opening, which is consistent with the previous analysis. The results show that coupling between fluid and geomechanics occurs in the process of initiation, extension, and crossing. Influenced by the dynamic change of stress, the NF opening shows a trend of increasing first, then decreasing, and finally reaching a constant size. With the increase of NFs, the flow distributed by NFs will decrease, and the stress interference between NFs will increase, leading to the decrease of fracture opening size and proppant migration ability. Therefore, if the NFs in the reservoir are extremely developed, they are not necessarily beneficial to the hydraulic fracturing operation.

| Effect of HSD on fracture geometry
The HSD has great influence on the fracture network geometry. To study the influence of HSD on a fracture network, we established a double-seam synchronous hydraulic fracturing model, as shown in Figure 10. The reservoir matrix has the characteristics of shale with staggered development of bedding, and the 3D CZM is embedded into the rock to simulate the NF network of nearly orthogonal geometry. The length, width, and height of the model are 20 m, 20 m, and 10 m, respectively. The main input parameters in the simulation are shown in Tables 1 and 2.

F I G U R E 8
Three-dimensional fracture geometry of two cases at the end of fracturing In the simulation, three HSD cases (0 MPa, 1.5 MPa, 3 MPa) are simulated, while the injection rate and injection time are kept constant. To increase the displayed effect, the crack network opening data are reprocessed by MATLAB program software. Figure 11 is a three-dimensional view of the crack network at different times under different HSD conditions; the color band in the figure represents the crack opening size. The results show that a smaller HSD produces a more complex fracture network and a longer total fracture length, as shown in Figure 12. Although the fracture network is formed, the complex fracture network is mostly concentrated near the horizontal wellbore, and the transformation in the direction perpendicular to the horizontal wellbore is not deep enough. However, when the HSD is increased to 3 MPa, complex fractures are formed, but not a complex fracture network (The complex fracture network reflects the existence of branching and connecting fractures between HFs and HFs. Complex fractures represent the existence of branching cracks, but there are no connecting branching cracks between HFs and HFs). Although no fracture network forms near the wellbore, more fracturing fluid drives the fracture deeper into the reservoir. The reason for this is stress interference between cracks. When the initial HSD is small, due to the effect of stress interference (the stress increases more quickly in the direction of the minimum horizontal main stress than in the direction of the maximum horizontal main stress), the in situ stress field changes and the maximum and minimum horizontal main stresses reverse locally, and fracture network spreading along the horizontal wellbore occurs (when the HSD is 0 MPa and 1.5 MPa, as shown in Figure 11). However, when the initial HSD is large, even with stress interference, stress inversion does not occur, and the crack mainly expands along the direction of the original maximum horizontal principal stress. From a production perspective, for reservoirs with a F I G U R E 9 The opening curve of four cohesive elements on NF-b in Cases I and II small HSD, fracture networks tend to form near the wellbore after fracturing. This type of well has a high initial productivity, but the capacity decreases quickly because the fracture network control volume is not deep enough. In addition, because the fracture network extends along the wellbore, the injected proppant tends to accumulate near the wellbore rather than migrate deep into the reservoir. For reservoirs with a large HSD, complex fractures form more easily than complex fracture networks. After fracturing, the initial productivity increase of such wells is limited and the original purpose of volumetric fracturing is not achieved.

| Effect of injection rate on fracture network geometry
The ultimate goal of hydraulic fracturing is to improve the recovery of oil and gas resources, rather than pursuing the complexity of the fracture network. For the field fracturing designer, the most effective, fastest, and cost-effective means of affecting the fracture network geometry is to change the injection rate of the fracturing fluid. Based on the physical model established in Figure 10, we simulated the fracture network with different injection rates with HSDs of 0.5 MPa and 3 MPa. Figure 13A,B,C represent the fracture geometry when the HSD is 0.5 MPa and the injection rate is 0.005 m 3 /s, 0.01 m 3 /s, and 0.015 m 3 /s, respectively. It is observed that with a low HSD, a complex fracture network can be formed after fracturing, and injection rate has a great influence on the fracture network. When the injection rate is lower, the stress interference effect and the stress reversal effect are smaller. In this case, the maximum length of the fracture network (SRV) is larger and the width is relatively small, as shown in Figure 13A,D. When the injection rate is high, the stress interference is strong. Fractures are more likely to deflect and extend along the horizontal wellbore. Thus, the maximum width of the SRV is larger and the length is smaller. From the production point of view, when the reservoir HSD is small, it is not advisable to increase the displacement because the fracture network that develops around the horizontal wellbore results in rapid decrease of productivity after fracturing. Figure 14A,B,C show the fracture geometry with different injection rates and an HSD of 3 MPa. The results show that when the HSD increases, a complex fracture network cannot be formed at any injection rate. When the injection rate is low, only complex fractures can be formed, as shown in Figure 14A (for the production of an unconventional reservoir, such complex fractures cannot greatly improve productivity). When the HSD is 3 MPa, with the increase of injection rate, it is observed that the fracture turns more easily and expands the SRV width; the SRV range is not limited to the vicinity of the wellbore. Therefore, from a production perspective, for reservoirs with a high HSD, it is beneficial to increase the injection rate appropriately. Originally, with a high HSD, fractures tend to be simple and straight, characterized by a small increase in production capacity in the early stage and stable production in the later stages. After increasing the injection rate, the complex fracture network can be constructed using strong stress interference and horizontal principal stress reversal, thus greatly improving the initial production capacity. Figure 15A,B,C show the fracture geometry with a stress difference of 3 MPa and fluid viscosity of 1 mPa s, 20 mPa s, and 40 mPa s. Figure 15D shows the relationship between the length and width of SRV impact range and the fluid viscosity. It is observed in Figure 15A,B,C that as the fluid viscosity increases, the geometry of fractures near the wellbore becomes more complex, and the fractures are more likely to turn and expand along the wellbore. The results in Figure 15D indicate that increasing the fluid viscosity and increasing the displacement have similar effects on the SRV. As the fluid viscosity increases, the maximum length along the y-axis decreases, while the maximum width along the x-axis increases. From a production perspective, complex geometric fractures near the wellbore are easily formed by hydraulic fracturing with high viscosity and can greatly improve initial productivity. Low viscosity fracturing fluids extend fractures deep into the reservoir, providing a more long-term, stable source of oil, and gas production. For reservoirs with a low HSD, the fracturing fluid viscosity can be appropriately reduced to form complex fractures with long-term productivity. For reservoirs with a high HSD, the complex fracture network can be formed near the wellbore, taking full advantage of the strong stress shadow effect if the fracturing fluid viscosity is appropriately increased.

| CONCLUSIONS
Based on the 3D CZM, a hydraulic fracturing model in a fracture pore elastic reservoir is established in this study. From a coupled fluid flow/geomechanics perspective, the stress interference effect between NFs, the effect of HSD, injection rate, and fluid viscosity on fracture geometry were studied. Based on the research results, some suggestions are provided for construction strategy. The intersection behavior of the HF and the NF was simulated. The simulation results agree with Blanton's criterion, proving the reliability of the 3D CZM. The results show that: 1. In the process of hydraulic fracturing, in addition to the strong stress interference between HFs, there is also strong stress interference between NFs. This stress interference causes the NF openings to dynamically increase and then decrease. With an increasing number of NFs, the stress interference is stronger. 2. For a given injection rate, a complex fracture network of a near horizontal wellbore forms more easily with a low HSD; the complex fracture network of near straight fracture forms more easily with a high HSD. 3. For reservoirs with a low HSD, higher injection rate and fluid viscosity produce a smaller maximum SRV length and a larger maximum SRV width. The fracture network tends to develop near the wellbore, which can cause rapid production failure and proppant accumulation. Therefore, for such reservoirs, maximum injection rate and fluid viscosity are not recommended. With the sand carrying capacity of fracturing fluid, the influence of stress interference should be minimized. 4. For reservoirs with a high HSD, lower injection rate and fluid viscosity produce a larger maximum SRV length and a smaller maximum SRV width. The fracture geometry tends to develop simple straight fractures, leading to a smaller initial production capacity. For such reservoirs, a higher injection rate and fluid viscosity are recommended to make full use of the stress interference effect to increase the development of the complex fracture network near the horizontal wellbore.
Research results presented in this study can be applied to hydraulic fracturing production practice. These results can guide hydraulic fracturing design, productivity prediction, and proppant migration analysis, and promote the understanding of fluid flow/geomechanics coupling.