Staged multicluster hydraulic fracture propagation mechanism and its influencing factors in horizontal wells for CBM development

Staged multicluster fracturing (SMCF) is crucial for the low‐permeability coalbed methane (CBM) reservoirs development. To reveal the fracture propagation mechanism and characteristics in SMCF and select appropriate hydraulic fracturing parameters, this study considers CBM exploitation in the Shizhuang block, Qinshui Basin in China as the engineering background. Then, XSite hydraulic fracturing software based on the discrete lattice theory and synthetic rock mass method is used to construct a three‐dimensional model. The influence mechanisms of the cluster spacing (CS) and lateral pressure coefficient (LPC) on the fractures propagation characteristics were investigated. With an increase in the CS, the stress interference between fractures gradually weakens, and the fracture deflection also decreases. Both the stimulated reservoir area (SRA) and fracture aperture of a single cluster increase with an increase in CS, but the SRA and fracturing efficiency exhibit limited growth when the CS is greater than 15 m. The SRA decreases exponentially with an increase in LPC. When the LPC increases from 0.4 to 1.4, the SRA decreases by 74.2%. The LPC plays a decisive role in the fracture propagation direction. When the LPC reaches 1.4, the fractures propagate along the horizontal direction. With an increase in LPC, the fracturing fluid pressure increases linearly, and the fracturing efficiency decreases exponentially. The results of this study reveal the fracture propagation law of SMCF and can provide a reference for engineering applications. In the SMCF at the 3# coal seam of the Shizhuang block, the CS can be set to approximately 15 m, and the CS can be appropriately reduced when the LPC is large.


| INTRODUCTION
Hydraulic fracturing technology is crucial for improving coalbed methane productivity, and it is widely used in China. [1][2][3][4] Staged multicluster fracturing (SMCF) methods have gradually become crucial for low-permeability coalbed methane development. [5][6][7] In SMCF, the hydraulic fracture morphology is complex, 8 and has a stress interference. Thus, understanding the hydraulic fracture propagation behavior and its influencing factors during SMCF is important to optimize the engineering construction design.
With the wide application of hydraulic fracturing technology, scholars have conducted numerous studies on fracture propagation in horizontal well fracturing using numerical simulation methods. For example, the boundary element method (BEM), finite difference method (FDM), finite element method (FEM), displacement discontinuity method (DDM), discrete element method (DEM), and lattice modeling method (LMM) have been applied to hydraulic fracturing. [9][10][11][12][13][14] Olson and Taleghani 15 first confirmed that fractures generate a surrounding induced stress field and proposed the concept of a stress shadow. Cheng et al. 16 analyzed the characteristics and basic principles of horizontal well fracturing from the perspective of rock mechanics, which provided theoretical guidance for horizontal well fracturing. Wong et al. 17 confirmed that there is an obvious stress interference between adjacent fractures. The minimum cluster spacing (CS) without stress interference is further determined by Morrill and Miskimins. 18 East et al. 19 proposed that induced stress generated by fractures could be reasonably used to promote fracture network formation to improve the fracturing range. Pan et al. 20 simulated the fracture morphology in hydraulic fracturing and analyzed its influence factors. Saberhosseini et al. 21 investigated fracture propagation in SMCF by a fully fluid-solid coupled model. Qu et al. 22 established a two-dimensional (2D) model of SMCF to analyze the effects of horizontal stress and CS on fracture propagation. Li et al. 23 established a mathematical model for SMCF simulation in terms of the initiation pressure, fracture width, and CS. Zhao et al. 24 and Liu et al. 25 used XSite to simulate the 3D fracture propagation behavior under different CS, stresses, and other parameters.
In summary, many scholars have studied the hydraulic fracturing of oil-gas wells and rock masses, but relatively few studies have considered the effect of SMCF in horizontal wells on CBM hydraulic fracturing. Due to the influence of burial depth and geological conditions, the distribution of reservoir stress varies greatly. 26,27 However, few studies have been conducted on the effect of lateral pressure coefficient (LPC) and CS on coalbed hydraulic fracturing. This study considers the 3# coal seam of the Shizhuang block in the Qinshui Basin. XSite hydraulic fracturing software developed based on the discrete lattice model (DLM) and synthetic rock mass (SRM) method is used to simulate the 3D hydraulic fractures propagation to analyze the fracture propagation mechanism under different LPC and CS conditions. This study can help to optimize the design of CBM fracturing parameters in this area.
The study area was selected from the Shizhuang Block, Qinshui Basin. The location and structural outline of the region are shown in Figure 1. The main developing coal seam, 3# coal seam, is stable with an average thickness of 6.1 m, and its. The average vertical stress, σ v , is 20 MPa in the region, and the variations are small and stable. The average maximum horizontal principal stress is 17.21 MPa and the variation range is 6.42-41.86 MPa. The minimum horizontal principal stress is 3.30-26.40 MPa, with an average of 11.39 MPa. The overall LPC is between 0.4 and 1.4, with an average of 0.8. In-site stress distribution of studied coal seam in different depths is shown in Figure 2. In the production process, it was found that the hydraulic fracturing effect is very sensitive to horizontal principal stress. Therefore, the influence of LPC and CS on hydraulic fracture propagation was studied in this paper.

| LATTICE MODELING METHOD
It should be noted that coal seam has obvious anisotropic characteristics. [28][29][30] However, in terms of engineering scale simulation, it is difficult to give the distribution of real fractures in coal seams, and the discrete fracture network fracture model is mostly used. [31][32][33] Lattice modeling method used in this paper was developed based on the SRM in DLM. 34 SRMs ( Figure 3A) consist of two components: the smooth joint model for the discontinuities mechanical behavior, and the bonded particle model describing the intact rock fracturing and deformation. The DLM ( Figure 4B) is a simplified version of the particle flow code, it consists of nodes with mass and springs, which has a higher computational efficiency than the particle flow code. 36

| Mechanical model
Two springs between lattice nodes in Figure 2 are representing the shear contact stiffness and the normal contact stiffness, respectively. And it can break and form microcracks when its strength is exceeded. These microcracks develop continuously and finally form macroscopic fractures. An explicit numerical method was used to solve the motion for all nodes, the Equation (1), (2), and (3) can be used to describe model displacement and spring force change. 35,37  ω ω where u i t ( ) and u̇i t ( ) are the position and velocity, I is the moment of inertia, ω i is the angular velocity, M i t ( ) is the sum of all moment components, k S and k N are the shear and normal stiffnesses of spring, F is the spring force, m is the node mass, Δt is the time step. "S" denotes "shear," and "N" denotes "normal," respectively in the equations.

| Flow model
Fluid elements linked by pipes are used for flow simulation in hydraulic fracturing, which is positioned at the centers of the broken springs or springs overlapped by the joint (Figure 4). The flow pipe network is dynamic and automatically updated by connecting newly formed microcracks to the existing flow network. lubrication equation is used to describe the fluid flows in pipes 38 (Ivars et al., 2011).
where β is the calibration parameter depending on the resolution and pipe connectivity; a is the fracture aperture; ρ w is the fluid density; g is the gravity acceleration; µ f is the fluid viscosity; s is the saturation, Z A and Z B are the heights of node A and B, respectively; Δp are the fluid pressures difference between nodes A and B, it can be calculated as: where Q is the sum of the flow rates; V is the volume of the fluid element, and K f is the bulk modulus of a fluid element.

| Model validation
Lattice modeling method was used to simulate the multifracture hydraulic fracturing test described by Xia et al. 39 In the process of hydraulic fracturing laboratory, the compressive strength of fractured sandstone is 59.3 MPa, the tensile strength is 9.5 MPa, the Elastic modulus is 21.5 GPa, and the Poisson's ratio is 0.17. The dense sandstone collected on-site is processed into a 300 mm × 300 mm × 300 mm test sample, with drill holes at the center of the test sample with a drill bit with an outer diameter of 25 mm. The CS in the test was set to 15, 30, and 45 mm, respectively. The three-dimensional stress difference is set to 2 and 6 MPa, respectively. In the model validation simulation, the model size and parameters are consistent with the laboratory model. The comparison between laboratory tests and numerical simulation results of the fracture propagation with different CSs were shown in Figure 5. It shows that the numerical simulation results are consistent with the laboratory results. When the CS is small, owing to the large stress interference between fractures, the fracture direction deviates in the maximum principal stress direction. With increasing CS, the stress interference and the fracture deflection decreases gradually, which indicates that the deflection angle of the fractures is clearly affected by the CS. The fracture propagation morphology of the numerical simulations and laboratory tests under different stress differences were shown in Figure 6. At specific CS, the fracture deflection degree decreases as the horizontal stress decreases from 8 to 4 MPa. The fracture propagation radius of the three clusters is basically consistent, and the fracturing efficiency is clearly improved. This indicates that the decrease in horizontal stress causes a decrease in stress interference. Therefore, the lattice modeling method proposed in this study can accurately reflect the actual hydraulic fracture propagation behavior in the laboratory and can provide a basic method for multicluster fracturing research.

| Model setup
The hydraulic fracturing well is modeled as a group of Ushaped horizontal wells, which are connected by a directional production well (LDP-02D) and a horizontal well (LDP-02H). The horizontal borehole strikes 70°NW, and the wellbore structure and main parameters are shown in Figure 7. The basic model established according to the fracturing well is shown in Figure 8. As shown in Figure 8, the length of each segment of the horizontal well is generally within 100 m, and about 3-5 fracturing clusters are arranged, and the CS is less than 30 m. In consideration of the simulation efficiency and the purpose of this paper is to compare the effects of different CS and pressure coefficients on the fracturing effect, the model size was 30 m × 6 m × 6 m to simulate one segment of the horizontal well. The model contains 96,960 nodes and 660,367 springs. The horizontal well diameter in the model was 0.14 m. Multiple clusters were established in the horizontal well, and its radius was 0.2 m. The initial fracture radius of the fracturing clusters was 0.5 m. Water was used as the fracturing fluid for synchronous fracturing. The fluid was injected into the clusters at a constant rate of 6.0 m 3 /min according to the actual parameters on site. 22,40 In the simulation process, the mechanical calculation for 0.1 s has been carried out first to reach the mechanical equilibrium.

| Effect of CS
The development of ultra-low permeability shale gas has benefited from the development of horizontal-well SMCF technology. 41,42 This study aims to introduce this technology in low-permeability CBM production. CS is the key parameter for dense SMCF in horizontal wells.
Research on the CS of horizontal wells has mainly focused on the stress redistribution state of different CSs caused by hydraulic fracturing under different reservoir conditions. [43][44][45] The distance between shale gas horizontal wells in the United States is generally 60-100 m, with five fracturing clusters in each stage. The distance between the clusters is generally 10-15 m, and in the densest clusters, the distance is 3-4 m. The CS in domestic shale gas horizontal wells is generally 20-30 m, with a minimum of 17.8 m, a maximum of 43 m, and 2-4 clusters in each stage. [46][47][48][49] The daily CBM production of ultra-low permeability shale gas horizontal wells generally increases with a decrease in CS and an increase in sand addition. Therefore, to study the fracture propagation mechanism for different CSs in the SMCF process in horizontal wells, CSs of 5, 10, 15, 20, and 25 m were simulated in this study. In the simulation, when the CS is 5 m, the model is arranged with five clusters, and three clusters for the 10 m CS. Because the average LPC of the 3# coal seam is 0.8, the stress states selected for this section were σ v = 20 MPa and σ x = σ y = 16 MPa. The other parameters were the same as those listed in Table 1.
During SMCF in horizontal wells, there is a certain relationship between the induced stress field and the CS. Different CSs will affect the induced stress field distribution, which then affects the fracture morphology. When the CS is small, more fracturing clusters are required in the horizontal well with the same length, and thus the fracture will be influenced by the superimposed stress during fracture propagation. Branch fractures will deviate from the original direction, increasing the fracture propagation network complexity. With increasing CS, the number of clusters will be reduced, and the stress interference will gradually be weakened. The hydraulic fracture propagation is thus affected by the original stress, and its direction tends toward the original maximum principal stress direction.
To analyze the interaction between hydraulic fractures, the fracture aperture ( Figure 9) is projected onto the XZ and the YZ plane ( Figure 8). Thus, the YZ plane cloud diagrams in the paper are a superposition of several cluster apertures. Figure 10 shows the XZ projection of the fracture aperture with different CSs. When the CS is small, the stress superposition effect between multiple fractures is strong, the fractures exhibit obvious deflection, and the fracturing efficiency is relatively low. With increasing CS, the stress superposition effect is weakened, and the stress interference has  Figure 11 shows the fracture aperture projection onto the YZ plane, the black solid line is the boundary of the fracture aperture of 0.1 mm. As shown in Figure 11, when the CS is less than 15 m, both the fracture aperture and fracturing range increase gradually. However, when the CS is greater than 15 m, the increase in the fracture aperture and fracture propagation range is limited. During hydraulic fracturing, changes in fluid pressure and ground stress can cause fracture development, resulting in damage to the area around the fracture (including tensile failure and shear failure), and the area of the failure area is the "stimulated reservoir area" (SRA). SRA for the fracturing area is adopted as the evaluation index for the hydraulic fracturing effect in this study. Figure 12 shows the relationship between the SRA of a single cluster and the CS. The SRA increases with an increase in the CS. When the CS is less than 15 m, the SAR of a single cluster is significantly affected by the CS. As the CS increases from 5 to 15 m, the SAR increases by 126.67%, from 8.85 to 20.06 m 2 . When the CS increases from 15 to 25m, the increase in the amplitude of the SAR is only 9.17%. Thus, when the CS is greater than 15 m, the SAR of a single cluster is not significantly affected by the CS, and with increasing CS, the permeability enhancement of the coal seam between clusters decreases gradually, which affects the gas extraction of the entire coal seam. Figure 13 shows that the initiation fracturing time of the coal matrix under different CSs is essentially equal. Fractures begin to propagate in a very short time at the beginning of hydraulic fracturing. Figure 14A shows that a large fracturing pressure is required at the beginning of fracture propagation. With the continuous constant fluid injection, the fracture continues to propagate. However, the propagation of the fracture also reduces the increasing amplitude of the fracturing water pressure, which gradually stabilizes. When the hydraulic fracture begins to fracture, the fracturing area is linearly correlated with the time (Figure 13). The fracture tip fluid pressure initially increases and then decreases with an increase in the CS ( Figure 14B). This is because a smaller CS leads to a larger interference force between the fractures and deflection of the fractures, and thus a smaller fracturing fluid pressure can promote fracture development. As the CS increases, the stress interference between the fractures weakens, and the pressure required for fracturing increases. As the CS increases to a certain value, the fracturing pressure should gradually stabilize. However, in Figure 14B, when the CS is greater than 15 m, the fracturing pressure decreases to a certain extent, which may occur because the cluster is relatively close to the model boundary and there is a certain boundary effect. Therefore, in this study, a CS of 15 m was selected for the subsequent analyses, considering the fracturing efficiency and cost ( Figure 14).  Table 1.

| Effect of the LPC
As shown in Figure 15, different LPCs cause significant differences in the fracture aperture. When the LPC is 0.4 or 0.6, the fractures are vertical, and the fractures range in the XZ plane is exhibiting a roughly rugby morphology. With an increase in the LPC, the range of the vertical fractures becomes narrow and begins to exhibit an obvious deflection. In particular, when the LPC reaches 1.4, the horizontal stress is significantly greater than the vertical stress, restricting vertical fracture propagation. Thus, the fractures exhibit visible deflection and development in the horizontal direction. Figure 16 shows the projection of the YZ plane fracture aperture under different LPCs. The results show that the larger the LPC, the smaller the fracturing range will be. In addition, the basic characteristics are similar to those of the XZ section. For LPCs of 0.4 and 0.6, the fracturing range is roughly the same, but with an increase in the LPC, the fracturing range decreases markedly. However, with LPC increases from 0.4 to 1.4, the fracture aperture in the middle region of the model increases gradually. It mainly occurs because the increase in the LPC results in an increase in the fracturing pressure. This causes the near fracture to be subjected to higher water pressure, and thus the fracture aperture increases significantly. Figure 17 shows the variation in the SAR with the different LPCs. The variation in the SRA decreases exponentially with an increase in the LPC, as shown in Equation (6).
where S is the SRA, and γ is the LPC. Compared with the CS, the LPC has a greater impact on the SRA. When the LPC increases from 0.4 to 1.4, the SRA decreases by 57.93 m 2 , representing a decrease in amplitude of 74.2%. Figure 18 shows the SAR-time curves corresponding to different LPCs. The initiation fracturing time is slightly different with different LPCs. With larger LPCs, the corresponding initiation fracturing time lags slightly. The fracture initiation times of the different LPCs are 1.0646,  Figure 18). This is mainly because this paper adopted the constant flow rate injection method. The higher the initiation fracturing pressure, the longer the corresponding initiation time will be. With an increase in the LPC, the curve slope of the SAR-time also decreases exponentially from 31.24 to 20.60, 16.04, 13.36, 11.68, and 10.44, respectively. This indicates that the LPC, especially the horizontal stress, has a considerable influence on hydraulic fracturing efficiency. However, unlike the effect of CS, the fracturing tip fluid pressure is proportional to the LPC ( Figure 19). Compared with Figure 18, it can be seen that the increase in the LPC increases the fracturing pressure linearly, while the fracturing efficiency decreases exponentially. Therefore, in the case of a high LPC, to improve the fracturing effect, it is necessary to increase the fracturing pressure or decrease the CS.

| CONCLUSION
(1) Based on the engineering background of CBM exploitation in the Shizhuang block, Qinshui Basin in China, a 3D hydraulic fractures model in a horizontalwell SMCF process is constructed based on the SRM in DLM. This model is verified by the hydraulic fracturing test results of a rock mass with different CSs.
(2) With the increase in the CS, the stress interference between the clusters gradually weakens, and the fracture deflection also decreases. The SRA and fracture aperture increase with the increase in CS. When the CS increases from 5 to 15 m, the SRA increases by 126.67%; however, it increases by only 9.17% when the CS increases from 15 m to 25 m. The fracturing efficiency (SRA/s) also shows that when the CS increases from 5 to 15 m, the fracturing efficiency increases by 126.5%, while the growth is limited for CSs of greater than 15 m.
(3) The LPC significantly affects the SRA of the SMCF. The SRA decreases exponentially with an increase in the LPC. When the LPC increases from 0.4 to 1.4, the SRA decreases by 74.2%. The fracture direction is also affected by the LPC. When the LPC is small, the fractures propagate along the vertical direction with no obvious deflection and have a good fracturing effect. When the LPC reaches 1.4, the excessive horizontal stress leads to the fractures propagating in the horizontal direction. In addition, an increase in the LPC causes the fracturing fluid pressure of the coal seam to increase linearly, and the fracturing efficiency decreases exponentially.
(4) Considering the fracturing efficiency and fracturing cost comprehensively, to ensure that the coal seam within the CS has a sufficient pressure relief range, the CS can be set to approximately 15 m in the 3# coal seam of the Shizhuang Block with horizontal well SMCF. The CS is appropriately adjusted according to the LPC of the location, and the CS can be appropriately reduced when the LPC is large.

AUTHOR CONTRIBUTIONS
Cun Zhang performed the lattice modeling method and drafted the manuscript. Yixin Zhao conducted the numerical simulation and conceived, designed, and coordinated the study. Ziyu Song mainly carries out numerical simulations. All authors gave their final approval for publication.