Study on Propagation Behaviors of Hydraulic Fracture Network in Tight Sandstone Formation with Closed Cemented Natural Fractures

Natural fractures in tight sandstone formation play a significant role in fracture network generation during hydraulic fracturing. This work presents an experimental model of tight sandstone with closed cemented preexisting fractures. The influence of closed cemented fractures’ (CCF) directions on the propagation behavior of hydraulic fracture (HF) is studied based on the hydraulic fracturing experiment. A field-scaled numerical model used to simulate the propagation of HF is established based on the flow-stress-damage (FSD) coupled method. This model contains the discrete fracture network (DFN) generated by the Monte-Carlo method and is used to investigate the effects of CCFs’ distribution, CCFs’ strength, and in-situ stress anisotropy, injection rate, and fluid viscosity on the propagation behavior of fracture network. The results show that the distribution direction of CCFs is critical for the formation of complex HFs. When the angle between the horizontal maximum principal stress direction and the CCFs is in the range of 30° to 60°, the HF network is the most complex. There are many kinds of compound fracture propagation patterns, such as crossing, branching, and deflection. The increase of CCFs’ strength is not conducive to the generation of branched and deflected fractures. When the in-situ stress difference ranges from 3MPa to 6MPa, the HF network’s complexity and propagation range can be guaranteed simultaneously. The increase in the injection rate will promote the formation of the complex HF network. The proper increase of fracturing fluid viscosity can promote HF’s propagation. However, when the viscosity is too high, the complex HFs only appear around the wellbore. The research results can provide new insights for the hydraulic fracturing optimization design of naturally fractured tight sandstone formation.


Introduction
Hydraulic fracturing is the key technology to improve the oil and gas development of the unconventional reservoir. The hydraulic fracture (HF) can effectively increase the migration channel of oil and gas and improve production and recovery [1][2][3][4]. The existence of the natural fracture (NF) system in the tight sandstone formation will change the HF's propagation path, which makes the HF generate many propagation patterns [5][6][7][8]. It is of considerable significance to understand the propagation mechanism of HF under the influence of NFs and reveal the formation condition of the complex HF network for the optimization of hydraulic fracturing in the tight sandstone formation.
Due to the lack of accurate and effective monitoring technology, the field cannot directly observe the HF's propagation process in the underground reservoir. At present, the propagation law of HF in unconventional reservoirs is not clear. Many scholars have studied the interaction between HF and NF by laboratory-scale fracturing experiments and numerical simulation. Hou et al., Zou et al., and Xie et al. [9][10][11][12] carried out a series of triaxial hydraulic fracturing experiments with natural shale. They investigated the influence of bedding planes on the HF's propagation behavior through acoustic emission (AE) monitoring technology. Tan et al. [13] classified the propagation pattern of HF in shale formation into four types: simple fracture, fishbone-like fracture, fishbone-like fracture with fissure opening, and multilateral fishbone-like fracture network. Zou et al. [14] carried out triaxial fracturing experiments on layered tight sandstone and studied the effect of natural fractures on CO 2 induced fracture propagation behavior. However, the fractures' distribution in the outcrop of natural rock is random, which is disadvantageous to the quantitative analysis of the fracture propagation mechanism. Therefore, some researchers tried to explore the interaction between HF and NF through the artificial fracturing specimen. In some studies, sandstone slice, gypsum slice, glass slice, or printing paper were placed in the cement mortar to simulate the NFs [15][16][17][18][19][20]. However, solid materials cannot reflect the friction effect and mechanical characteristics of closed NF. Some studies proposed to use the interface between artificial core blocks to simulate NF [21][22][23][24][25][26], but this kind of experiment model can only study the influence of a single preexisting fracture on the HF propagation path.
Due to the high cost and limited scale of laboratory experiments, many scholars have developed different numerical simulation methods to explore the HF propagation mechanism. Guo [27][28][29][30][31][32][33] used the extended finite element method (XFEM) to study the propagation process of complex HFs in fractured reservoirs. Wang proposed a global cohesive zone model and used it to study the formation conditions of the fracture network in fractured reservoir [34]. Zhou [35][36][37][38] used two-dimensional particle flow code (PFC2D) to simulate the effect of NF characteristics on the HF propagation behavior. Amir et al. combined the extended finite element method (XFEM) and discrete element method (DEM) to simulate the HF propagation in porous media with natural fractures [39]. Tang proposed a fully coupled flowstress-damage (FSD) model to simulate the damage and failure process of heterogeneous materials [40]. Many studies used the FSD model to study the propagation process of HF [41][42][43][44]. Li et al. [45] used the FSD model to simulate the initiation, propagation, and associated stress evolution during the fracturing process in the glutenite formation and studied the influence of the confining stress ratio, gravel sizes, and gravel volume content on the hydraulic fracturing pattern in a conglomerate specimen. Wang et al. [46] used the FSD model to study the influence of the injection rate on the propagation behavior of HF in the formation containing a discrete fracture network (DFN). Men et al. [47] simulated the influence of bedding planes' direction and mechanical characteristics on the propagation path of HF in the shale formation. Li et al. [48] used this method to simulate the multistage hydraulic fracturing process and studied the influence of fracture spacing and stress anisotropy on the propagation and reorientation of three HFs.
To explore the HF propagation behavior in tight sandstone formation containing multiple groups of CCFs, a new laboratory-scale fracturing experimental model is designed in this paper. The model not only has a closer mechanical property to tight sandstone but also contains multiple groups of closed cemented preexisting fractures. The hydraulic fracturing experiments with this experimental model are carried out to study the influence of fracture network direction on the propagation behavior of HF. Then a large-scale numerical model used to simulate the hydraulic fracturing process in naturally fractured tight sandstone formation is established by the coupled FSD method. The discrete fracture network (DFN) generated in the model is based on the Monte-Carlo method. Through the numerical simulation, the effects of CCFs' distribution, CCFs' strength, in-situ stress anisotropy, injection rate, and fluid viscosity on the propagation behavior of the HF network are investigated. A new artificial experimental model is used in this study to quantitatively analyze the influence of CCF network direction on HF's propagation. Compared with the existing experimental model, the new model not only contains several groups of CCFs with different directions but also has mechanical properties closer to that of natural tight sandstone. The preparation process of this model includes two key steps:

Hydraulic
Step (a): determination the mortar formula of artificial tight sandstone The materials for making artificial tight sandstone specimens include natural quartz sand, clay, epoxy resin, and densifier. The screening steps of the mortar formula are shown in Figure 1. According to different formulas, the mortar is mixed and poured into the mold. The mortar specimens are solidified under the compaction pressure for 48 hours and are maintained in the incubator for 24 hours. The artificial tight sandstone specimens with the closest physical and mechanical properties to the natural tight sandstone are selected through the mechanical property test and permeability test of each group of mortar specimens.
Step (b): simulation of CCFs in the hydraulic fracturing specimen The structure of the artificial tight sandstone fracturing specimen is shown in Figure 1. The specimen size is 300 mm × 300 mm × 300 mm. The specimen contains three layers, and each layer is 100 mm thick. L-A and L-C represent the upper and lower barriers, respectively. L-B is the hydraulic fracturing target layer, which contains closed cemented preexisting fractures in different directions. The preparation procedures are as follows: (I) the mortar is poured into the mold and solidifies under 40 MPa compaction pressure. After 48 hours, the upper barrier (L-A) is formed. (II) According to 2 Geofluids the formula determined in step (a), the mortar is prepared and poured into the mold which contains L-A. Then, the steel sheets (0.5 mm thick) are inserted into the mortar according to the design direction and position. (III) All steel sheets are taken out from the mortar slowly after L-B has solidified for one hour. The mold containing L-A and L-B is compacted at a pressure of 30 MPa for 48 hours. Then the preexisting fractures in L-B are closed and cemented under compaction pressure. (IV) The lower barrier is prepared according to step (I).
The angle between the CCF and the direction of horizontal maximum principal stress are defined as θ in this paper. Considering that the natural fractures in the formation cannot be in the same direction, the angle θ of CCFs in the three fracturing specimens are set to 0°± 15°, 45°± 15°, and 90°± 15°. The comparison of basic mechanical parameters of the artificial experimental model and natural tight sandstone is shown in Table 1.

Experimental
System and Procedure. The triaxial hydraulic fracturing physical simulation experiment system of Northeast Petroleum University is used in our experiment ( Figure 2). The injection rate is set to 5 ml/min, the horizontal maximum principal stress and the horizontal minimum principal stress are set to 20 MPa and 16 MPa, and the vertical stress is set to 25 MPa. Slick water mixed with red tracer is used as the fracturing fluid. The acoustic emission (AE) system is used to monitor the initiation and propagation of hydraulic fractures in the specimens during the experiment. Eight AE sensors were installed in the groove of the loading plates. The frequency-domain of the AE system was 125-750 kHz, and the gain of the channel preamplifier was 40 dB.

Analysis of Experimental Results
2.3.1. Analysis of Hydraulic Fracture Geometry. Figure 3(a) shows the HFs' propagation patterns observed at the interface between L-A and L-B. The yellow lines and red lines represent the CCFs that are not affected and communicated by hydraulic fractures, respectively. It can be seen that the direction of CCFs has a significant impact on the propagation behavior of HF. In specimen #1 (θ = 0°± 15°), all the four CCFs on the HF's propagation path are communicated by HF. Due to the small angle between the CCFs and the direction of the horizontal maximum principal stress, a simple crack with symmetrical wings is finally formed. In specimen #2 (θ = 45°± 15°), the HF interacts with eight CCFs and has obvious deflection behaviors. There are two possible propagation behaviors of HF after encountering CCFs: deflection along with the CCF and penetration through the CCF. In specimen #3 (θ = 90°± 15°), the HF penetrates the two CCFs near the wellbore with a reinitiation after a short extension distance along with the CCF. Subsequently, the HF crosses through the remaining CCFs and propagates to the direction of horizontal maximum principal stress. This shows that the geometry and propagation pattern of the HFs are the most complex when the angle between the horizontal maximum principal stress direction and CCPF is between 30°and 60°.

Analysis of Injection Pressure-Time Curves.
The characteristics of the injection pressure-time curve can be used to analyze the propagation behaviors of HF. As shown in Figure 4, the initiation and propagation process of HF can be divided into four stages: Stage (I): pressurization and fracture initiation stage. The injection pressure increases sharply with the continuous injection of fracturing fluid. When the injection pressure exceeds the sum of the rock tensile strength and the horizontal minimum principal stress, the rock surrounding the wellbore breaks and the initial HF appears Stage (II): main fracture propagation stage. The main HF propagates in the rock matrix. The injection pressure-time curve is relatively stable and has no significant fluctuation Stage (III): fracture interaction stage. The HF interacts with the preexisting fractures significantly, which makes the injection pressure fluctuate obviously Step (I) Step (II) Step (III) Step (     5 Geofluids decreased twice and then increased again. When the HF intersects with the CCFs, the CCFs dilate and open, and the fracturing fluid leaks into the opened fractures, which makes the pressure drop greatly. When the opened fracture is filled with fracturing fluid, the injection pressure increases again to meet the mechanical conditions of HF propagating in the rock matrix again. At the stage (III) of specimen #2, the injection pressure curve fluctuates violently and has many peaks. At this stage, the propagation patterns of HF include deflection and penetration. Under the joint action of the reinitiation of the rock matrix and the opening of preexisting fractures, the injection pressure appears violent fluctuation. At the stage (III) of specimen #3, when the HF interacts with the first set of preexisting fractures, the injection pressure only fluctuates slightly and lasts for a short time. In the subsequent fracturing process, the change of injection pressure is small. This shows that the HF mainly propagates in the rock matrix and has no significant interaction with the CCFs.

Field-Scale Hydraulic Fracturing Numerical Simulation
Due to the limitation of the simulation scale of laboratory experiments, only three groups of hydraulic fracturing experiments are carried out to study fracture direction's influence on HF's propagation. In this paper, the influence of other key factors on the complex fracture propagation law in fractured tight sandstone formation is studied through a twodimensional numerical model based on the FSD method.

The Coupled Flow-Stress-Damage (FSD) Numerical
Method. The FSD model considers the coupling effects of flow, stress, and damage on the permeability changes caused by the propagation of fractures and the evolution of rock damage.
The assumptions of the model are as follows: (1) The waterrock coupling relationship is based on the Terzaghi principle, the rock seepage conforms to Darcy's law, and the rock deformation process of fluid-solid coupling conforms to the Biot consolidation theory. (2) The material of the model element is assumed to be elastic brittle material with certain residual strength. (3) The damage of the material element follows the maximum tensile strength criterion and Mohr-Coulomb failure criterion. (4) Rock is heterogeneous material, and its physical and mechanical parameters follow Weibull distribution. The element in the FSD coupling model is defined in five states: elastic deformation, tensile failure, compression failure, tensile separation, and compression contact, as shown in Figure 5(a). The current state of each element is determined by its damage history, current stress-strain level, and material properties. Failure paths occur when elements damage, contact, and separate. When the stress state of an element satisfies any strength criterion, it begins to accumulate damage. In the elastic damage mechanics, the elastic modulus of the element decreases with the damage. Based on the assumption of strain equivalence, the elastic modulus of damage element is given by the following formula: where D is the damage variable. There is only one failure mode for the elastic brittle element with residual strength. The compressive stress should be positive, and the stress-strain relationship is shown in Figure 5(b). When the minimum principal strain of the element exceeds the uniaxial tensile strain or the minimum In the multiaxial stress state, the damage variable of the element in the tensile mode can be expressed as: In this case, after the element is damaged, the permeability coefficient of the element can be expressed as follows: When the shear stress of the element satisfies the Mohr-Coulomb failure criterion, it is considered as failure in shear mode.
In the multiaxial stress state, the damage variable of the element in the shear mode can be expressed as: In this case, after the element is damaged, the permeability coefficient of the element can be expressed as follows: where ε is the equivalent principal strain of the member; ε t0 is the elastic ultimate strain; ε t1 is the ultimate tensile strain when the element is completely damaged; ξ is the permeability damage factor; σ 1 and σ 3 are the maximum and minimum effective principal stress; σ t0 and σ c0 are the tensile failure strength and compressive failure strength of the element; σ rt is the residual tensile strength of the element; K is permeability coefficient. Compared with other hydraulic fracturing numerical simulation models, the FSD coupling model is based on the finite element and statistical damage theory. It considers the heterogeneity of material properties and can simulate the dynamic fracturing process of the heterogeneous reservoir.  Table 2.
3.2.1. Analysis of Hydraulic Fracture Geometry. It can be seen from Figure 6 that the fracture propagation patterns of numerical simulation are in good agreement with the experiment. The interaction stage (stage III) of the three numerical models is significantly different. The numerical model of specimen #1 finally appears a simple HF along the direction of horizontal maximum principal stress, which is consistent with the experimental results. It is worth noting that before the HF intersects with the CCF, some damage elements have appeared on the CCF. With the increase of water pressure, the HF encounters the CCF, and the number of damaged elements on the CCF increases rapidly. The HF connects the whole CCF in about 5 time steps. In the numerical model of specimen #2, the HF interacts with the first set of CCFs and appears the composed propagation pattern including penetration and deflection. After the branching fractures are formed in the numerical model, the propagation direction of some HFs changes, and the geometry of HFs becomes complex. In the numerical model of specimen #3, the HF deflects slightly at the first set of CCFs and then propagates into the rock matrix with a short deflection distance. Then, the HF penetrates through the remaining CCFs. In this process, only a few damage elements appear near the interaction point, and the CCFs remain closed. Figure 7, the initiation pressure of specimen #1 is 23.4 MPa. The injection pressure decreases sharply and increases again at step 40 and step 70, respectively. These are the moments when the HF interacts with the CCFs. Due to the normal stress on the surface of CCF in specimen #1 is small, the connected CCFs are completely damaged. Thus, the filtration of the fracturing fluid into the preexisting fracture is significant, and the injection pressure curve is greatly reduced.

Analysis of Injection Pressure-Time Curves. In
The initiation pressure of specimen #2 is 24.4 MPa. The injection pressure fluctuates sharply from step 20 to 114. During this stage, the complex branching propagation fractures appear because of many kinds of interaction patterns between HF and CCFs. The violent interaction is accompanied by frequent filtration and pressure holding, which results in the violent fluctuation of the injection pressure curve. At step 114, the injection pressure has a sharp decline, which is the moment when the HF intersects with the second set of CCFs and deflects along with the opened fractures. After step 120, although the injection pressure fluctuates, the overall pressure is lower than that of step 20-114. During this stage, the HF not only deflects along the opened CCFs but also continue to expand in the rock matrix.
Compared with the other two groups of experiments, the initiation pressure (26.6 MPa) and propagation pressure of specimen #3 are both larger. Besides, the whole fracturing process of specimen #3 lasts a long time (276 steps), which is due to the HF mainly propagates in the rock matrix. At step 78, the injection pressure decreases slightly and rises rapidly. At this time, the HF intersects with the first set of CCFs, and the injection pressure decreases under the effect of filtration. However, due to the large normal stress on the CCFs in specimen #3, the resistance of CCF to open is large. Thus, the HF repenetrates the two CCFs after a short extension distance along with the CCF. After step 90, the injection pressure is relatively stable, which shows that HF penetrates through the remaining CCFs and propagates in the rock matrix.

Numerical Model of Hydraulic Fracturing with Discrete Fracture Network (DFN).
To study the influence of other key factors on the propagation behavior of HFs in tight sandstone formation with CCF network, a field-scale numerical model is established based on the FSD coupling method, as shown in Figure 8. The model has a size of 400 m × 600 m and includes 240000 elements. A two-dimensional DFN generation program is generated based on the Monte-Carlo method. According to the fracture geometry parameters obtained from the field statistics, the DFN is inserted in the numerical model. The fracture parameters of the stochastic modeling method mainly include length, azimuth, opening, and density. This paper focuses on the influence of CCFs on HFs' propagation, so the influence of fracture opening is not considered. The width of CCFs is equivalent to the width of an element. There are two groups of CCFs in the numerical model, and the relevant parameters are in Table 3. The azimuth angle θ is the angle between the direction of CCFs and the direction of horizontal maximum principal stress. In the DFN model, the lengths of CCFs satisfy the normal distribution, and the azimuth angles of CCFs satisfy the logarithmic normal distribution. The physical and mechanical parameters of the numerical model are shown in Table 2.
In the middle of the model, there is a perforation parallel to the direction of the horizontal maximum principal stress with a length of 5 m.

Influence of CCFs Distribution Direction.
In the numerical model, the azimuth angle (θ) of the two groups of natural fractures are set to ±15°, ±30°, ±45°, ±60°, ±75°, ±90°, and the standard deviation is 0. The injection rate is 10m 3 /min, the As shown in Figure 9, when the azimuth angle θ is less than 30°, the HFs mostly propagate along the CCFs with only a small deflection. When θ is small, the normal stress acting Step 12 Step 40 Step 70 Step 134 Specimen #1 (a) Specimen #1 Step 16 Step 50 Step 114 Step 208 Specimen #2

(b) Specimen #2
Step 30 Step 78 Step 134 Step 276 Specimen #3 (c) Specimen #3   9 Geofluids on the CCFs is small. After the HF intersects with the CCFs, the CCFs are easy to open. Under the induction of opened CCFs, the HF propagates along the direction of the fracture development direction. When the azimuth angle θ is ±45°, the deflection amplitude of HF is larger when it encounters the CCFs. Besides, there are some branching fractures in the HFs' propagation path, and the fracture geometry after fracturing is complex. When the azimuth angle θ is ±60°, the HFs deflect sharply along the CCFs. There are many kinds of compound fracture propagation patterns, such as penetration, branching, deflection, and offsetting. When the azimuth angle θ is ±60°, the fracture geometry after hydraulic fracturing is the most complex. When azimuth angle θ is ±75°, the deflection amplitude of the HF along the CCFs becomes smaller. The propagation pattern of HF mainly penetrates or repenetrates the CCFs after a short extension distance along the CCF accompanied by a small number of branching fractures. Figure 10 shows the injection pressure-step curves of numerical models. The initiation pressure of HF tends to increase gradually with the increase of the azimuth angle θ. The smaller the θ is, the easier the HF will deflect along the CCFs, and the lower the required propagation pressure is. On the contrary, the larger the θ is, the easier the HF will cross through the CCFs and propagate in the rock matrix.
In Figure 11, the total equivalent hydraulic fracture length (EL) and propagation area of complex fracture network (PA) increase with the increase of the azimuth angle θ from 15°to 60°. When θ exceeds 60°, the EL and PA decrease again. When θ ranges from 45°to 60°, the normal stress acting on the CCF is small and the shear stress along the CCF is large. The CCFs have complex fracture behaviors such as dilation, opening, and shear sliding, and the HF has a variety of propagation patterns such as penetration, deflection, and offsetting. Therefore, when θ ranges from 45°to 60°, both the EL and PA can meet the requirements of increasing oil and gas production. Table 3. Under the condition that other parameters remain unchanged, the tensile strength of CCFs is set to 0.1 MPa, 0.4 MPa, 0.7 MPa, 1.0 MPa, 1.3 MPa, and 1.6 MPa, respectively. Figure 12 shows a comparison of four representative fracture propagation patterns. The HF network's complexity after hydraulic fracturing decreases with the increase of tensile strength of CCFs from 0.1 MPa to 1.6 MPa. When the strength of CCFs is small, CCF is easier to be activated during the hydraulic fracturing process, which promotes the formation of deflecting or branching fractures. When the tensile strength of CCFs is large, the stress produced at the tip of HF is not enough to drive the activation and opening of CCFs. In this case, the possibility and distance of HF deflection are small, and the complexity of the fracture network after hydraulic fracturing is low.

Influence of CCFs Tensile Strength. The establishment of DFN based on the parameters in
In Figure 13, the tensile strength of CCFs has little effect on HF's initiation pressure when the azimuth angle of CCFs is unchanged. The initiation pressure ranges from 36.3~37.9 MPa. However, the CCFs' tensile strength has an impact on HF's propagation pressure. The propagation pressure increases with the increase of CCFs' tensile strength. In naturally fractured reservoirs, the HFs mainly deflect and propagate along the CCFs or branch at the intersection point. Thus, the greater the tensile strength of the CCFs, the greater the resistance to fracture propagation that needs to be overcome.

Geofluids
In Figure 14, the EL and PA have a negative correlation with the tensile strength of CCFs. Under the same stress condition, the resistance of CCFs open and shear slip decreases as the tensile strength of CCFs decreases. The smaller the tensile strength of CCFs is, the greater the activation degree of CCFs is, and HFs are easier to communicate with more branching propagation fractures. Table 3. The numerical model sets the minimum horizontal principal stress (σ h ) to 30 MPa and changes the maximum horizontal principal stress (σ H ) to adjust the anisotropy of in-situ stress. The horizontal in-situ stress difference of the model is as follows: 0 MPa, 1.5 MPa, 3.0 MPa, 4.5 MPa, 6.0 MPa, and 7.5 MPa. Other parameters of the numerical model remain unchanged. Figure 15 shows a comparison of four representative fracture propagation patterns. The horizontal in-situ stress difference (Δσ) has a significant effect on HFs' propagation patterns. When Δσ is 0 MPa, the in-situ stress has no control over the fracture propagation direction. The HFs mainly propagate along the CCFs after the intersection, and the overall propagation direction of HFs is consistent with the distribution direction of CCFs. At this time, although the fracture morphology is complicated, the propagation length of HF is limited. When Δσ is 1.5 MPa, the control of in-situ stress anisotropy on the propagation direction of HFs increases, but most HFs still propagate along the CCFs with a large deflection angle. When Δσ is between 3 MPa and 6 MPa, the HF and the CCF have obvious interaction behaviors, and the geometry of fracturing fractures is complex and the propagation length is long. When Δσ exceeds 6 MPa, insitu stress has a strong control of fracture propagation. The deflection behavior of HFs decreases while the penetration and offsetting behaviors increase. Although HF has a good propagation length, the complexity of fracture morphology decreases after hydraulic fracturing.

Influence of In-Situ Stress Anisotropy. The establishment of DFN based on the parameters in
In Figure 16, the initiation pressure of HF has a slight increase from 37.7 MPa to 38.6 MPa with the increase of Δσ.

Geofluids
The initiation of HF mainly overcomes the sum of horizontal minimum principal stress and rock tensile strength. In this group of simulations, σ h is a fixed value, so the variation of initiation pressure is very small. However, the propagation pres-sure of HF increases with the increase of Δσ. The increase of σ H will increase the normal stress acting on the surface of the CCF, and the resistance to be overcome when the HF propagates along CCFs will increase. Step  In Figure 17, both the EL and PA increase first and then decrease as Δσ increases from 0 MPa to 7.5 MPa. The EL reaches the maximum when Δσ is 4.5 MPa, while the PA reaches the maximum when Δσ is 3 MPa. This shows that when Δσ is less than 3 MPa, the propagation of HF is dominated by the direction of CCFs, which cannot guarantee the effective propagation length. However, when Δσ exceeds 6 MPa, the propagation of HF is dominated by in-situ stress, and the degree of HF deflecting and branching is low. In Figure 18, when the injection rate is 4~6 m 3 /min, the HF interacts with CCFs and forms branching and deflecting propagation patterns, but the degree of branching and deflecting of the fracture is small. When the injection rate increases to 8 m 3 /min, the interaction between HFs and CCFs is strengthened, and HFs' geometry morphology is more complicated. On the whole, the hydraulic fracture network after fracturing becomes more complicated with the increase of the injection rate of fracturing fluid.
In Figure 19, the initiation pressure of HF increases from 39.7 MPa to 44.3 MPa with the increase of the injection rate. Besides, the higher the injection rate, the faster the pressure rise rate in perforation and fracture initiation rate. The HF's propagation pressure is positively correlated with the injection rate, while the hydraulic fracturing duration is negatively correlated with the injection rate.   13 Geofluids In Figure 20, the EL and PA gradually increase with the increase in the injection rate. However, when the injection rate is increased from 12 m 3 /min to 14 m 3 /min, the increase rate of the EL and PA decreases. It shows that an appropriate increase in the injection rate of fracturing fluid is conducive to improving the hydraulic fracture network's complexity and sweep    Step   In Figure 21, when the fluid viscosity is 20 mPa·s and 40 mPa·s, the HFs' propagation pattern is intricate, and the crisscrossing network fractures are formed in the model.  But when the fluid viscosity exceeds 60 mPa·s, the complexity of HFs gradually decreases with the continuous increase of fluid viscosity. In Figure 22, the initiation pressure of HF gradually increases from 36.15 MPa to 40.05 MPa with the increase of fluid viscosity. The higher the fluid viscosity is, the faster the initiation speed of HF is. There is a good positive correlation between fracture propagation pressure and fracturing fluid viscosity. When the fracture is stably extending, the high viscosity fracturing fluid will encounter greater resistance when flowing in the fracture, and larger injection pressure is required to maintain the propagation of the HF.
In Figure 23, both the EL and PA increase when the fluid viscosity increases from 20 mPa·s to 40 mPa·s. But when the fluid viscosity exceeds 60 mPa·s, the EL and PA gradually decrease after hydraulic fracturing. This shows that in the naturally fractured tight oil reservoir, appropriately increasing viscosity can reduce the filtration loss of fracturing fluid and promote fracture propagation. However, when the viscosity is too high, the resistance of fracturing fluid in the fracture is too large, which limits the fracture propagation and eventually leads to the complex fracture network only form in a limited range around the perforation.

Conclusions
In this paper, a new laboratory fracturing experimental model is designed, which not only has mechanical properties closer to natural tight sandstone but also contains several groups of closed cementing preexisting fractures. Based on exploring the influence of fracture network direction on hydraulic fracture propagation through hydraulic fracturing experiments, a field-scaled hydraulic fracturing numerical model of naturally fractured tight sandstone formation is established by combining FSD coupling model and Monte-Carlo simulation method. The effects of natural fracture distribution direction, natural fracture tensile strength, insitu stress anisotropy, fracturing fluid injection rate, and fracturing fluid viscosity on the propagation of complex fracture networks are studied by this numerical model. It can be concluded as follows.
(1) The experiment results show that when θ is 45°± 15°; the HF has two kinds of propagation behaviors when interacting with CCFs: penetration and deflection, which leads to the most complicated HF geometry. The numerical simulation results show that when θ ranges from 45°to 60°, the CCFs on the propagation path of HF include a variety of complex fracture behaviors such as dilation, opening, and shear sliding. This makes the HF appear in many propagation patterns such as penetration, deflection, branching, and offsetting. The numerical simulation results are in good agreement with the experimental results (2) The tensile strength of CCFs has a significant influence on HF's propagation patterns. As the decrease   effective propagation length of hydraulic fracture HF is limited. When the horizontal stress difference exceeds 6 MPa, the HF propagation is dominated by in-situ stress anisotropy, and the complexity and spread area of HFs are limited. Therefore, when the horizontal stress difference of tight sandstone formation is within the range of 3~6 MPa, the propagation length and spread area of the HF network can be guaranteed simultaneously (4) The increase of the injection rate of fracturing fluid is helpful to increase the complexity and spread range of HF network, but the improvement degree of fracturing modification effect by increasing the injection rate is no longer obvious when it exceeds 12 m 3 /min (5) An appropriate increase in the viscosity of fracturing fluid can reduce the filtration loss and promote the expansion of the HF. However, when the fluid viscosity is too high, the resistance of fracturing fluid in the fracture is too high, which will lead to the complex fracture network only appearing in a limited range around the perforation Nomenclature HF: Hydraulic fracture NF: Natural fracture CCF: Closed cemented fracture FSD: Flow-stress-damage coupling model DFN: Discrete fracture network θ: The angle between horizontal maximum principal stress and preexisting fractures Δσ: Horizontal in-situ stress difference EL: Total equivalent hydraulic fracture length PA: Propagation area of the complex fracture network K: Permeability coefficient D: Damage variable ξ: Permeability damage factor ε t0 : Elastic ultimate strain ε t1 : Ultimate tensile strain σ t0 : Tensile failure strength σ c0 : Compressive failure strength σ rt : Residual tensile strength of the element.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declared that they have no conflicts of interest to this work.