Effects of Fracturing Parameters on Fracture Unevenness During Large-Stage Multi-Cluster Fracturing in Horizontal Wells

Horizontal wellswith multi-cluster fracturing technology is an effective approach to exploit unconventional hydrocarbon reservoirs. The on-site diagnosis results indicate that multi-cluster fractures always tend to propagate unevenly due to stressinterference, therefore it is very essential to study the effect of fracturing parameters on fracture propagation unevenness. In this paper, the unconventional fracturing model (UFM, Unconventional Fracturing Model) is used to study the effect of multi-cluster fracturing parameters on fracture unevenness in a large stage. This model has been validated with the actual fracturing case on-site in the Longmaxi shale. The investigated parameters include completion parameters (cluster spacing, number of perforations per cluster), pumping parameters (fluid injection intensity and proppant injection intensity). Our simulation results show that firstly reducing fracture spacing will increase stress interference, andhydraulic fractures exhibit a “radial” pattern. Secondly, reducing the perforation number of a single cluster can promote the more uniform propagation of multi-cluster fractures. Thirdly, increasing the fluid injection intensity will increase the fracture length, but will also increase the fracture unevenness. Besides, the injection strength of the proppant has a little effect on the average fracture length and the unevenness of the fracture length. Finally, setting a reasonable cluster spacing and injection fluid strength can obtain a more uniform fracture propagation. Meanwhile reducing the number of perforations per cluster can also reach the goal of propagating evenly. This paper provides a certain reference for the optimization of multi-cluster fracturing parameters in large-stage and multi-cluster wells.


INTRODUCTION
The unconventional reservoirs such as shale have nanoscale pore-spaces and very low permeability,so the traditional reservoir stimulation methods cannot satisfy their commercial exploitation (Jia et al., 2012). In recent decades, horizontal well drilling with multicluster fracturing technology has become the mainstream technology for the exploitation of unconventional reservoirs.
The implementation of this technology is to shoot multicluster perforations in the casing of a horizontal well and then to inject high-rate and high-pressure fracturing fluid into the wellbore. Multiple hydraulic fractures can initiate and propagate simultaneously, which increases the stimulation volume and greatly improves reservoir seepage conditions (Chen et al., 2010;Wu et al., 2011a,b).
With the increase in the burial depth of shale gas reservoirs, problems such as complex stratigraphic structure and large insitu stress differences will emerge in the reservoirs, which makes it difficult to form multiple uniform and effective hydraulic fractures in the target formation (Xie, 2018;Xie et al., 2019;Ma et al., 2020). With the background of the low oil prices, the large-stage multi-cluster fracturing technology, increasing the stage length to reduce the amount of bridge plugs and reducing the cluster spacing to improve the reservoir contact, is the key technology to reduce the cost and improve production (Fan et al., 2019). However, recent studies have shown that reducing fracture spacing will perform stronger stress interference and more uneven fracture length are obtained (Pan et al., 2014;Zhao et al., 2015;Li et al., 2017;Zhou et al., 2019). Pan et al. (2014) established a three-dimensional finite element model to simulate the simultaneous propagation of multiple fractures in horizontal wells. Their results show that the number of perforation clusters and fracture spacing is the main factors affecting fracture propagation. The closer perforation clusters are distributed, the more serious the stress interference in one stage can be obtained. Zhao et al. (2015) established the simultaneous propagation model of multiple fractures using the displacement discontinuity method, and they studied the influence of the stress interference between fractures on the local stress field and initial fracture pressure. Their results show that induced stress may lead to the asymmetric propagation of the two wings of the fractures. Li et al. (2017) studied the influence of perforation erosion on the fracture propagation of multiple clusters by a comprehensive model that couples perforation erosion effect and hydraulic fracture propagation. Their results show that the proppant in the sand slurry will destroy the perforations, resulting in a huge limit-entry capacity reduction, which in turn leads to a wider variation among the length of fractures in one stage. Liu et al. (2020) explored multi-fracture propagation in shale strata based on the Lattice model, and their results show that middle fractures are suppressed by other adjacent fractures. Fracture spacing and in-situ stress difference are the key parameters to keep uniform fracture length from their results. Meanwhile their results (Liu et al., 2019) show that the application of zipper fracturing technology increases SRV when two horizontal wells exist. These above studies indicate that uneven fracture propagation is very common during the propagation of multiple fractures.
In multi-cluster fracturing field operations in the large stages, fracturing fluids with temporary plugging agents are usually used to plug the over-treated fractures and promote the uniform propagation of multi-cluster fractures (Wang et al., 2018;Murphree et al., 2020;Zhou et al., 2020). Zhou et al. (2020) studied the effect of temporary plugging ball on fracture propagation using a fully coupled model of "wellbore-perforation-fracture propagation" based on the boundary element method. Their results show that the addition of temporary plugging agents helps to improve the non-uniform propagation of multi-cluster fractures. Li et al. (2020) studied the temporarily plugging staged fracturing (TPSF) using the cohesive zone method, and their results show that lateral fractures can break through the suppression from previous fractures and TPSF enhance the complexity of fracture network. Chen et al. (2020a) numerically investigated the optimization of the number of diverters and diverting time during nearwellbore diversion in a heterogeneous stress reservoir, and their results indicate that more ball sealers and the earlier diverting time are required for creating a new fracture in the formations with the high-stress zone. Some advanced monitoringmethods in the oil-field also show the effectiveness of temporary plugging technology. Rahim et al. (2017) adopted post-job diagnostics and analyzed well performances in the cases in carbonate and sandstone reservoirs, and they proved that the diverters can successfully plug perforations. Panjaitan et al. (2018) carried out a comprehensive method including water hammer analysis, step-down test, and micro-seismic data in Haynesville Shale and their results show the diverters can divert the fluid to intended perforations and TPSF can promote uniform propagation of multiple fractures. Some other field researches also get similar results (Huang et al., 2018;Senters et al., 2018;Weddle et al., 2018). In TPSF operations, a key question is how uneven is the fractures at this point when the temporary plugging agents are injected. Fracturing parameters have an important impact on the unevenness of multi-cluster fracture propagation. Therefore, studying fracturing parameters on multi-cluster fracture length is of great significance to the design of the timing and dosage of the temporary plugging agents.
At present, there are few studies on the influence of multicluster fracturing parameters on the unevenness of multicluster fracture lengthin a large stage. In this paper, numerical simulation is conducted using the unconventional fracturing model (UFM) to analyze the sensitivity of fracturing design parameters.The first part of this study will briefly describe the mathematical model of the UFM model; In the second part, a large-scale multi-cluster fracturing simulation will be performed on the shale formation of the Longmaxi Formation in Changning, and compared with the real treatment pressure curve to verify the accuracy of the model; In the third part, numerical simulations will be carried out on different fracturing parameters to study the influence of fracturing parameters on the unevenness of fracture length. The main fracturing parameters studied include completion parameters (cluster spacing and the number of perforations per cluster); pumping injection parameters (fluid injection intensity and proppant injection intensity).

MATHEMATICAL MODEL
The UFM model is used to simulate the propagation of multiple clusters of fractures. The UFM model is a complex fracture model developed by Weng et al. based on a pseudothree-dimensional plane hydraulic fracture model. The model can simulate the propagation of multiple fracture tips. The mathematical model is as follows (Nolte, 1991;Weng et al., 2011Weng et al., , 2012:

Fluid Flow in Fractures
The flow equation of power-law fluid in a fracture follows Poiseuille's law (Nolte, 1991): where, p is the fluid pressure; q is the local flow velocity in the fracture; h fl is the fluid height in the fracture; w is the average fracture width; w(z) is the depth-related fracture width; s is the distance along the fracture; n and K is the fluid power-law exponent and consistency coefficient, respectively.
The fluid in the fracture obeys the conservation of mass, and its continuity equation is: where c 1 is the total filtration coefficient; h 1 is the filter loss zone height; τ(s) is the time when the fracture unit first contacts the fracturing fluid. Global volume balance equation: where Q(t) is the injection rate; L(t) is the total length of the entire fracture at time t; h is the fracture height.
In the wellbore, the total injection rate should be equal to the sum of the rates flowing into the fractures, and its equation satisfies (Mack and Elbel, 1992): where, q i (t) is the injection rate of the ith fracture. The 2D KPN model is used to describe the relationship between fracture width and pressure, and the equation is satisfied   (Perkins and Kern, 1961;Nordgren, 1970): where E is the plane strain rock modulus. At the crack tip, the following boundary conditions are met:

Fracture Height Control Model
For the control of the fracture height, the UFM model uses a model like the pseudo-three-dimensional model (Weng et al., 2011). The pressure p in the fracture is: where p cp is the fracture pressure at the perforation depth (measured at the bottom of the fracture); ρ f is the fluid density. The fracture height is obtained by matching the stress intensity factor and fracture toughness. The intensity factors at the top and bottom of the fracture tip can be obtained from the pressure in the fracture, fracture height, and local stress. The calculation equation for the stress intensity factor and fracture width at the top and bottom of the fracture tip is: The fracture spacing is 5 meters The fracture spacing is 10 meters The fracture spacing is 15 meters The fracture spacing is 20 meters where, K Iu and K Iu is the stress intensity factors at the top and bottom of the slit tip, respectively; δ n and δ i is the local stress at the top of the fracture tip and the bottom of the fracture tip, respectively; h is the fracture height; h i is the height from the bottom of the seam tip to the top of the seam tip.

Stress Interference Model
The 3D modified boundary element method is used to calculate the induced stress field (Crouch and Starfield, 1983;Olson, 2004): where A ij is the 3D correction factor; h is the fracture height; d ij is the distance between i th element and j th element; α and β is the fitting parameters (α = 1, β = 3.2); C ij is the influence coefficient of the discontinuous displacement of element j on the induced stress on element i; C ij ns indicates the influence coefficient of the shear discontinuous displacement of element j on the normal stress of element i; C ij nn indicates the influence coefficient of the discontinuous displacement of the element j on the normal stress of the element i.

Proppant Transport Model
The model determines the type of fluid and proppant by volume concentration, and the calculation equation for the average volume concentration in the fracture is as follows: horizontal advection migration. The fluid concentration change equation is described as (Adachi et al., 2007): where q fl is the flow rate in the fracture; f leakoff is the velocity of the fluid passing through the fracture wall; The sedimentation velocity of the proppant in the vertical direction is described by the Stokes equation (Weng et al., 2011): Where v set,k is the sedimentation velocity of the k proppant; ρ prop,k is the density of the k proppant; D k is the diameter of the k proppant; n and K is the rheological index and consistency coefficient weighted by the fluid concentration, respectively.

Interaction Criterion Between Hydraulic Fractures and Natural Fractures
Natural fractures can be seen as the friction interface, Renshaw and Pollard (Weng et al., 2012) established the orthogonal interaction criterion of hydraulic fracture and natural fracture, and Gu et al. (Nolte, 1991) extended to the nonorthogonal criterion on Renshaw and Pollard's study. The detailed work can be seen in Gu's study. With a certain cohesion interface standard, the extended orthogonal formula is: where S 0 is the cohesive force; µ is the interfacial friction coefficient; T 0 is the tensile strength of the rock; σ H , σ h 1 2 3 4 5 6 7 8 9 1 0 1 2 3 4 5 6 7 8 9 1 0 1 2 3 4 5 6 7 8 9 1 0 1 2 3 4 5 6 7 8 9 1 0 are the horizontal maximum and minimum principal stresses, respectively.

Model Establishment and Model Verification
The UFM has been integrated into the Mangrove platform of the geological integrated software Petrel. The verification of the model and the analytical solution has been conducted by Weng et al. (Murphree et al., 2020), so this study will not repeat.
In this section, based on the same/real logging data, injection procedures, completion measures, and fracturing design as well A, three-dimensional fracturing numerical simulation results are generated and compared with the on-site fracturing treatment pressure curve to verified the accuracy of the geological model and fracturing model.
As shown in Figure 1, the target zone of Well A is located in the Long 1 layer of the Longmaxi Formation in Changning, with a burial depth of 3060-3175 m. Logging interpretation results show that organic carbon content (TOC) is between 4 and 8%. is between 71.34 and 83.69 MPa. And the vertical stress (Sv) is between 71.83 and 73.17 MPa. The stress state belongs to a strike-slip fault. The measured depth of Well A is 4891.47 m, the vertical depth is 3100∼3200 m, the horizontal completion section length is 1500 m, the average stage length is designed to be 120 m, the total number of perforations in one stage is 48 holes, and the number of clusters is 11 clusters and the cluster spacing is 10.9 m. The average minimum horizontal, maximum horizontal and vertical stress is 64.19 MPa, 77.52 MPa and 72.5 MPa, respectively. X-axis, Y-axis, and Z-axis are the directions of horizontal minimum stress σh, horizontal maximum stress σH and vertical stress σV, respectively. In addition, the horizontal wellbore direction in this model is parallel to the direction of the minimum horizontal principal stress, so that multiple fractures can be generated in the horizontal wellbore perpendicular to the wellbore. In fact, this is a common practice in multi-cluster fracturing of horizontal wells. Sensitivity analysis is proposed to determine how uniform will the fractures distribute based on changes in certain fracturing parameter while the other parameters remain the basic model parameters. The basic model parameters are shown in Table 1. Figure 2 shows the on-site treatment pressure and the simulated pressure are in good agreement. It should be noted that there are still some pressure disturbances in the pumping pressure curve, which are mainly ascribed to stratigraphic causes. However, it can be proved that simulating the Longmaxi section of the Changning block by this model can be fairly representative to some extent.
This section will study the fracturing parameters from the following two aspects: completion parameters (cluster spacing and the number of perforations per cluster); pumping parameters (fluid injection intensity and proppant injection intensity).
The unevenness of fracture length is defined as the length difference in the geometrical morphology of multiple clusters of fractures in the same stage. In this study, the fracture propagation geometry mainly refers to the fracture length, so we define the value of the fracture unevenness as the standard deviation of the fracture length of multiple clusters. The larger the value is, the more uneven the distribution of multiple clusters of fractures in the stage; otherwise, the more uniform fracture propagation can be obtained when the value is smaller. Its mathematical expression is: where σ (L) is the length unevenness of multiple clusters of fractures; L i is the fracture length of the ith fracture; L is the average fracture length of multiple fractures.

Fracture Spacing
The fracture spacing will be changed to 5 m, 10 m, 15 m, and 20 m, respectively, to explore the influence of the cluster spacing on the length unevenness of multi-cluster fracture propagation. As shown in Figure 3, when the fracture spacing is 5 m, 10 clusters of fracture propagate in a "radial shape", and the strong stress interference under the narrow cluster spacing promotes the deflection of the clusters. With the continuous increase of the fracture spacing, when the fracture spacing is 20 m, the morphology of each fracture is almost parallel, which indicates that the increase of the cluster spacing can reduce the stress interference and the fracture morphology transitions from a strong fracture deflection to a mutual parallel fracture group. Fracture lengths of the cases in Figure 3 are analyzed, and the length unevenness of the cracks and the average value of the fracture lengths in each case according to formula (20) are calculated. Figure 4A shows the fracture lengths in different fracture spacing cases. It can be seen from the figure that as the fracture spacing increases, the histogram of fracture length distribution gradually turns from uneven to even. Figure 4B shows the variation of the average fracture length and unevenness with the different fracture spacing. The results showed that the fracture spacing increased from 5 to 20 m, the unevenness of the multi-cluster cracks decreased, and the average fracture length did not change much.

Number of Perforations per Cluster
With the basic parameters unchanged, the number of singlecluster perforations is modified to 2 holes, 4 holes, 6 holes, and 8 holes, respectively. Influence of the number of perforations per cluster on the fracture propagation morphology and length unevenness are explored. The fracture propagation pattern is shown in Figure 5. It can be seen from this figure that most of the fractures have a fracture width of less than 2 mm, and the overall crack width is small. The side fractures are affected by stress interference and the fracture length is reduced. The middle fractures are deflected due to stress interference, and multiple fractures merge into the main fracture within 100 m from the wellbore and continue to expand, which further increases the induced stress and promotes uneven fracture propagation.
The fracture lengths of the cracks in each case in Figure 5 are counted, and the fracture length results are shown in Figure 6A. The results in Figure 6A do not reflect the changes in the uniformity of the cracks. Therefore, the unevenness calculation is performed, and the calculation result is shown in Figure 6B. It can be seen from Figure 6B that the average length of the fracture has a slight downward trend with the increase in the number of perforations per cluster. The unevenness of the fractures increases with the increase in the number of perforations per cluster. When the number of perforations per cluster is six, it is an abnormal point. This may be due to the influence of the calculation of the fracture length when multiple clusters of fractures merge, which interferes with the calculation result of unevenness. These results show that reducing the number of single-cluster perforations can help reduce the unevenness of cracks and improve the uniformity of cracks; at the same time, increase the average length of fractures slightly.

Fluid Injection Intensity
With the basic parameters unchanged, the pump injection schedule is changed to alter the fluid injection intensity as: 10 m 3 /m, 20 m 3 /m, 30 m 3 /m, respectively. Influence of the injection fluid intensity on the fracture propagation morphology and length unevenness are explored. Figure 7 is the fracture morphology diagram. From Figures 7A-C, as the fluid injection intensity increases, the difference in fracture length of multiple clusters of fractures gradually increases, and the fracture propagation gradually becomes uneven, such as the longest fracture length in Figure 7A is 319.05 m, and the shortest fracture length is 141.894 m. In Figure 7B, they are 497.189 m and 102.252 m, respectively. This shows that the increase of the injection strength increases the unevenness of the fracture propagation. Figure 8A shows the fracture length distribution under different fluid injection intensity. Figure 8B shows the variation of unevenness and fracture length with fluid injection intensity. When the fluid injection intensity is 10 m 3 /m, the average fracture length is 221.7 m, and the unevenness is 51.53 m; while the injection intensity is 20 m 3 /m, 30 m 3 /m, the average fracture length is 278.7 m, 327.2 m, the unevenness of the cracks increased to 104.58 m and 88.02 m; This shows that the increase of the injection intensity increases the length of the fracture stimulation as a whole, but the unevenness of the fracture gradually increases, and there is the extreme value of the maximum unevenness of the fracture. In our simulations, the injection intensity in the case with the maximum unevenness is 20 m 3 /m.

Proppant Injection Intensity
In this section, we will keep the basic parameters unchanged, and change the proppant injection intensity to 1.8 t/m, 2.4 t/m, and 3.0 t/m, respectively, and explore the influence of the proppant injection intensity on the fracture propagation pattern and the unevenness of fracture length. The fracture morphology diagram is shown in Figure 9. It can be seen from Figures 9A-C that when the proppant injection intensity increases from 1.8 to 3.0 t/m, the average fracture width has been increased, and the average fracture width is 2.1 mm, 2.4 mm, and 3.0 mm, respectively. There is not much change in other geometric forms of fractures. Figure 10A shows the distribution of fracture lengths under different proppant injection intensity. Figure 10B shows the variation of average fracture length and the unevenness with the proppant injection intensity. The results show that when the injection proppant intensity increases, the average length of the fractures in each cluster hardly changes, which are 278.7 m, 285.8 m, and 285.4 m, respectively; the unevenness of the cracks gradually decreases to 104.6 m, 90.9 m, and 85.9 m, respectively. This shows that the change of the proppant injection intensity will not affect the average length of fractures or the scale of stimulation and meanwhile have a little effect on improving the uniformity of the fracture length.

Effect of Approach Angle When Natural Fractures Exist
Complicated natural fractures (NF) and bedding usually exist in shale formations. Therefore, complex fracture networks are created when hydraulic fractures (HF) encounter natural fractures. The approach angle refers to the angle between hydraulic fractures and natural fractures, which is an important parameter that affects the interaction between hydraulic fractures and natural fractures. This section mainly describes the complex fracture network when the natural fractures are distributed under different approach angles. The approach angles for setting natural fractures are 30, 60, and 90 degrees. The basic setting parameters of natural fractures are shown in Table 2. Figure 11 shows the fracture networks with different natural fractures approach angles. Compared with the no-HF case, the fracture morphology with NFs shows a more complicated fracture network. Due to the high fluid loss of natural fractures, the local fracture width shows wider hydraulic fractures. In addition, the fracture morphology of different approach angles is also quite different. When the natural fracture angle is 30 degrees, hydraulic fractures divert easily to the natural fractures, which show a highly complex fracture network in the nearwellbore region and far-field region. When the approach angle is 60 degrees, hydraulic fractures are captured by the natural fracture only at the near-wellbore region and a small amount of the fracture network was disturbed in the far-field region. This phenomenon is more obvious when the approach angle is 90 • . There are only a few hydraulic fractures in the farfield region, and the figure shows that hydraulic fractures do not communicate with the natural fractures. The complex fracture network in the near-wellbore region is mainly due to the strong stress interference created by the simultaneous propagation of multiple fractures, which makes it easier for hydraulic fractures to divert to natural fractures. However, the stress interference is relatively weak in the far-field region. At low approach angles, hydraulic fractures in the far-field region can still divert to natural fractures to produce a more complex fracture network, as shown in Figure 11A; At high approach angles, the hydraulic fractures in the far-field region will directly cross the natural fractures, as shown in Figure 11C. Figure 12A shows the fracture length under the conditions of different approach angles. It can be seen from this figure that natural fractures with low approach angles promote the propagation of hydraulic fractures. Although natural fractures cause the high fluid loss of fracturing fluid, natural fractures 30-degree approach angle 60-degree approach angle 90-degree approach angle  provide a fast propagation path, which caused hydraulic fractures to spread to the far-field region. However, when the approach angle increases, natural fractures cut off the normal propagation path of hydraulic fractures to the far-field. Hence the shape of hydraulic fractures near the well section is complex, the length of hydraulic fractures gradually decreases. Figure 12B proves the above conclusion. When the approach angle is low, the average fracture length is the highest, the reservoir stimulation zone is the best, and the unevenness is the highest; when the approach angle is high, the average fracture length is the shortest and the unevenness is the lowest.

CONCLUSION
This study uses an UFM to carry out geological modeling and fracturing simulation of real Sichuan Changning shale. Comparing with field fracturing treatmentpressurecurvewith the simulationtreatmentpressurecurve, the accuracy of this simulation can be verified. Through a series of numerical simulations of different fracturing parameters, the main conclusions as follows: 1. Small fracture spacing has a strong stress interference during fracturing in the horizontal wells, which will cause obvious uneven propagation of fractures.Its influencing factors include completion parameters (cluster spacing and the number of perforations per cluster) and pumping parameters (fluid injection intensity and proppant injection intensity). 2. Reducing the fracture spacing will increase the stress interference, causing the fracture to appear "radial style, " Meanwhile the unevenness of the fracture's length increases with the decrease of the fracture spacing. 3. Reducing the number of perforations per cluster can promote more uniform propagation of multiple clusters of fractures and reduce the unevenness of fracture length. 4. Increasing the fluid injection intensity will increase the length of the fractures, but at the same time, it will increase the unevenness of the fracture length and cause the fractures to propagate more unevenly. 5. The increase of proppant injection intensity has little effect on the average length of fractures, while it has a little effect on improving the uniformity of the fracture length. 6. When natural fractures exist, the uniformity of fracture length increases with the decrease of approach angle in horizontal wells. 7. In horizontal large-stage multi-cluster wells, reasonable cluster spacing, and fluid injection intensity should be set to achieve a more effective stimulation effect, meanwhile reducing the number of perforations per cluster is also helpful to improve the uneven propagation.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.