Numerical simulation study on three‐dimensional fracture propagation of synchronous fracturing

A new three‐dimensional fluid‐solid coupling fracture propagation model of synchronous fracturing was established to study the dynamic propagation process of simultaneous hydraulic fracturing. In this paper, we used an improved three‐dimensional boundary element method and a finite difference method, respectively, for calculating the deformation of the reservoir rock and the flow of the fracturing fluid during the fracturing process. The interaction between hydraulic fractures and nature fractures was also considered. Sensitivity analysis showed that during synchronous fracturing, the fractures between the wells attracted each other and injected with more fluid due to interference between the fractures. Affected by natural fractures, the trajectory of hydraulic fractures was more tortuous. Besides, the width and height of hydraulic fractures captured by natural fractures were reduced. In order to reduce the difference in length between individual fractures, smaller well spacing and larger cluster spacing should be used.


| INTRODUCTION
In order to improve the development and utilization efficiency of fossil in unconventional reservoirs, hydraulic fracturing technology is widely used. 1,2 In unconventional reservoirs, conventional hydraulic fracturing technology cannot meet the development needs, so the new hydraulic fracturing technology and the new fracturing concept also have developed. For improving the fracturing efficiency, in recent years, the unconventional reservoir fracturing has been transformed from a single horizontal well fracturing to a multi-well fracturing in the well group, and the concept of synchronous fracturing has been proposed. 3 Synchronous fracturing is a new hydraulic fracturing technology which means that there are two parallel or more horizontal wells fracturing simultaneously.
Field practice shows that synchronous fracturing technology can increase production significantly and has a small impact on the environment of the fracturing zone. 4 Synchronous fracturing is essentially a problem of multiple fracture propagation at the same time. When multiple fractures are propagating simultaneously, the fractures will interfere with each other because of the existence of the induced stress field. [5][6][7] Hong studied the change of stress field under different fracturing processes, but the dynamic propagation of fractures was neglected. 8 To study the multi-fracture propagation problem, 9 established a multi-fracture propagation model based on the displacement discontinuity method. And the results showed that the fracture width will be significantly reduced due to the interaction between multiple fractures. 10  of stress interference between multiple fractures in shale formations and found that the stress interference between multiple fractures is obviously affected by fracture spacing and in situ stress, and the reduction of fracture spacing increases the stress interference between the two fractures. Based on the two-dimensional plane strain, 11 used the finite element method to study the fracture propagation path and plastic strain with different fracture cluster numbers. 12 analyzed the synchronous fracturing and zipper fracturing process according to commercial software Abaqus. However, those models ignore the effects of natural fractures. The existence of natural fractures has a significant impact on the propagation trajectory of hydraulic fractures. 13,14 It is found that when the hydraulic fractures propagate, the nearby natural fractures can be induced to open, and the closer between the hydraulic fractures and the natural fractures, the more obvious of the induction effect will be. 15 Weng 16,17 established an unconventional fracture propagation model (UFM) based on the classical fixed high fracture propagation model. The model considered the influence of natural fractures on the hydraulic propagation path through the fracture intersection criterion and used the displacement discontinuity method to characterize the fracture and stress interference during dynamic propagation. Since the hydraulic fracture and the natural fracture are randomly intersected, the UFM model can accurately simulate the propagation of the hydraulic fracture after encountering the natural fractures at any angle, which is close to the actual situation. 18,19 used the extended finite element method to simulate the hydraulic fracture propagation of shale reservoirs under the influence of natural fractures. Through simulation, it is found that when the angle between the natural fracture and the maximum horizontal principal stress is small, it is difficult to form a fracture network. And the greater the horizontal principal stress difference, the more difficult to form a complex fracture network. 20 established a fracture propagation model based on the boundary element method, and the interaction between random natural fractures and hydraulic fractures was also considered. Wang used RFPA software to analyze the fractures propagation law during synchronous fracturing under the influence of natural fractures. 21 It is worth noting that the above scholars did not study the propagation of hydraulic fractures in the direction of the fracture height. For hydraulic fractures, the height change during hydraulic fracturing was very obvious, 22,23 so it is necessary to study the three-dimensional fracture propagation model. Some scholars 24,25 have studied the three-dimensional fracture propagation, but they have not considered the influence of natural fractures. At present, there are models that can simulate the multi-fracture propagation, the interference of hydraulic fractures and natural fractures and the three-dimensional fracture propagation, but few models have considered the influence of all the above factors. Especially for hydraulic fracturing in unconventional reservoirs, the effects of natural fracture and the vertical propagation of fractures need to be considered together.
As the difficulty of fossil resource development increases, the accuracy of calculation of fracture height in hydraulic fracturing process is gradually increased. In order to study better the synchronous fracturing fracture propagation trajectory and fracture morphology in unconventional fossil reservoirs, a three-dimensional complex fracture propagation model which is suitable for unconventional reservoirs is established based on the theory of linear elastic mechanics and fracture mechanics. And the characteristics of unconventional reservoirs rich in natural fractures and interference between fractures have also been considered. After the numerical solution of the model, the synchronous fracturing fracture propagation was studied in the unconventional reservoir, and the influencing factors in the fracture propagation process were also analyzed.

| PHYSICAL MODEL
Simultaneous fracturing generally refers to the simultaneous propagation of multiple fractures in a well or multiple wells. But synchronous fracturing refers to the simultaneous fracturing of two wells or more than two wells, and multiple fractures propagation at the same time, as shown in Figure 1. The induced stress generated near the fracture will change the magnitude and direction of the in situ stress of the reservoir which can affect the propagation direction and fracture morphology of the adjacent fracture. Due to the change of ground stress, the morphology of the two wings of the fracture is no longer symmetrical. Therefore, we divide the fracture of each cluster into two asymmetric fractures. In order to facilitate the study, the model has been properly and reasonably simplified and the following assumptions are made: (a) Natural fractures are vertical fractures, 26

| Fracture propagation model
In this paper, an improved boundary element method is used to simulate fracture propagation, because the boundary element method only needs to discretize the boundary, which can reduce the dimension of the research object and the calculation cost, and improve the calculation accuracy. For the convenience of calculation, a fixed size discrete fracture grid is used, as shown in Figure 2.

| Rock deformation
Each vertical section of the fracture can be reduced to a line fracture in the plane strain problem, and these line fractures are independent with each other. The coordinate system is established as shown in Figure 3. Where h u and h l are the upper and lower seam heights of the fracture, h is the half height of the reservoir, H is the height of the fracture half, and z u and z l are the coordinate values of the intersection of the fracture cross section and the cover layer and the bottom layer, respectively. z d is the distance between the center of the fracture and the center of the reservoir. The net pressure acting on the fracture surface can be decomposed into the following forces: net pressure at the center of the fracture ( p net ), the frictional pressure drop of the fluid along the fracture height, (g v |z|), pressure difference generated by the gravity(g z), the geostress gradient (g s z), the stress difference between the cover layer and the reservoir (), and the stress difference between the bottom layer and the reservoir. (). The fracture width variation caused by the above various stress components can be expressed as w 1 , w 2 , w 3 , w 4 , w 5 , and w 6 , respectively. The actual width of the fracture is obtained by the superposition principle: The width of the fracture at any position on the vertical profile of the fracture can be calculated from the England and Green formulas:  E is the modulus of elasticity; and w is the fracture width. For the specific calculation equation of the component of each fracture width, please see Appendix A.

| Fluid flow
The hydraulic fracturing problem is a typical hydro-mechanical coupling problem, so the flow of fracturing fluid in the fracture also needs to be studied. The pressure in the fracture is unevenly distributed due to the existence of flow resistance in the fracture, and the pressure drop equation in the fracture can be obtained by the Navier-Stokes equation, 28 which can be expressed as 16 : Here, Q is the injection rate; w is the fracture width; n' is the fluid power-law index; and k' is the fluid viscosity index. Perforating friction is a crucial factor in fracturing construction and an important design parameter in the fracturing process, which seriously affects the distribution of fluids in the wells. According to the Bernoulli equation, the calculation formula of perforating friction 29,30,31 can be obtained as follows: Here, Q i is the injection rate of the i-th fracture; ρ s the density of the fracturing fluid; n p,i is the number of perforations of the i-th perforation section; dp i is the perforation diameter of the i-th perforating section; and C d,i is the hole correction factor of the i-th fracture. Due to the frictional resistance of the wellbore, there is also a pressure difference between the fracture clusters, which also causes the difference in flow injection between the fractures. The pressure drop due to the friction of the wellbore can be calculated by 32 : Here, C cf is the friction coefficient; x j is the distance from the fracture j to the end of the wellbore; Q k is the flow in the fracture k; and D is horizontal wellbore diameter. According to Kirchhoff's first and second laws, when multiple fractures are simultaneously propagation, it is necessary to satisfy the pressure balance and the total fluid injection mass conservation. Then, the flow conservation equation and the pressure balance equation of the bottom hole 33 are as follows: Here, Q T and Q i are the total fluid injection rate and the fluid injection rate of the fracture i, respectively; p o horizontal well bottom fluid pressure; p wi is the fluid pressure at the entrance of the i-th fracture; p pf,i is the friction of the perforation of the i-th fracture; and p cf,j is the j-stage wellbore friction. By the law of conservation of mass, the global mass conservation equation can be expressed as the sum of the volume of the fluid in the fracture and the volume of leak-off, as follows: Here, C L s the leak-off coefficient; t is the current injection time; and τ is the fracture element start leak-off time.

| Fracture propagation
In order to facilitate the calculation, the stress intensity factor theory 34,35 is used to calculate the fracture propagation in this paper, rather than the implicit level set algorithm. 36 Since the angular change of the fracture in the vertical propagation is neglected, the type II stress intensity factor is not zero at z = 0. Then, the length of the fracture propagation can be calculated by:

Where
Here, ▵ d max is the maximum propagation length of an element; K is a stress intensity factor; K I is the type I stress intensity factor; K II is the type II stress intensity factor; K max is the maximum stress intensity factor; K IC is fracture toughness of reservoir rock; E is Young's modulus; d is the half length of fracture element; D Tip n is the normal displacement of tip element; D Tip s is the shear displacement of tip element; and is the fracture steering angle. If Δd is greater than or equal to the length of the fracture element in the direction of fracture propagation, the adjacent elements in the direction of fracture propagation will be opened.

| Crossing criterion
The interaction between hydraulic fractures and natural fractures is the most common mode of mechanical action. In the process of hydraulic fracture propagation, the damage at the fracture tip is not a simple tensile failure, so the stress field at the fracture tip calculated by a single type I fracture factor is not accurate. The fracture tip stress field can be obtained by introducing a type II fracture factor as follows 37,38 : Here, K I and K II are type I and type II stress intensity factors, respectively; S Hmax and S hmin are the horizontal maximum and minimum principal stress, respectively; and θ is the angle between hydraulic fracture and natural fracture, as shown in Figure 4. When tensile failure occurs in natural fractures, the fluid pressure in the natural fracture is the same as the fluid pressure in the fracture tip, and the following formula is required: Here, p nf is the fluid pressure in natural fractures; σ nf is the normal stress on the surface of natural fractures; and S n is the tensile strength of natural fractures. When shear failure occurs in natural fractures, the natural fracture element has no normal displacement, and the tangential stress boundary condition needs to satisfy the friction law 39 : Here, τ nf is the shear stress of natural fractures; τ n is the shear stress of natural fracture surface caused by far field stress; and K f is the friction coefficient of fracture surface.

| The fluid-solid coupling approach
Hydraulic fracture propagation is a problem of fluid-solid coupling. In this paper, the displacement discontinuity method and the Palmer model are used for the fracture propagation and the fracture width calculation, and the finite difference method is used for the calculation of the fluid flow in the fracture. In order to ensure the accuracy of the model calculation, the full coupling method is adopted. And in order to improve the convergence of the model, the weak coupling form is used for the coupling calculation of the flow field and (13) the stress field. The solution process of the hydraulic fracturing fracture propagation model in this paper is shown in Figure 5.

| Model comparison and verification
The hydraulic fracturing fracture propagation model proposed in this paper is a three-dimensional model. In order to verify the correctness of the model, the new model is compared with the three-dimensional analytical model 40,41 and numerical model. 24 The length of the fracture and the width of fracture at bottom hole in the three-dimensional analytical model can be calculated by the following formula: Where, Q is the injection rate; and μ is the viscosity. Since the analytical model cannot account for vertical in situ stress changes, the reservoir thickness is set to be greater than the fracture length and the detailed input parameters can be found in the paper by. 24 The mathematical model and the solution process established in this paper are programmed and solved, and the fracture width distribution obtained is shown in  it can be found that the shape of the fracture obtained by the model in this paper is more accurate. This is because the Yang model ignores the propagation of the vertical fracture elements in the horizontal direction. Analytical solutions are generally considered to be exact solutions. And the comparison between the model in this paper and the analytical solution is shown in Figure 7. It can be found that the improved model in this paper has a higher degree of agreement with the analytical solution.

| Number of clusters
In the development of fossil energy resources, horizontal well staged fracturing technology is widely used. The horizontal section of the horizontal well is divided into stages by splitters and successively fracturing, and each stage can be designed to contain one or more clusters of fractures. In this section, the single-stage single-cluster and singlestage two-cluster synchronous fracturing of two wells are simulated and analyzed, and the input parameters are shown in Table 1. Figure 8 shows the fracture geometry under different number of clusters. It can be found that in reservoirs with natural fractures developed, fractures between wells during synchronous fracturing attract each other. And in the case of a single-stage two-cluster, the fractures of cluster 1 and cluster 2 are longer than the fractures of cluster 3 and cluster 4, which is caused by the resistance in the wellbore. On the other hand, if there are more than one fracture in a single section, the hydraulic fractures will not only be affected by the stress generated by the fractures of adjacent well, but also by the adjacent fractures in the same stage, which will inhibit the propagation of the fracture cluster at the far end. It also can be found that the shape of the hydraulic fracture affected by the natural fracture will be more tortuous and the width and height distribution of the fracture will also change significantly. It is worth noting that the width of the fracture in the single-stage multi-cluster fracturing is greater than that in the single-stage single-cluster. This is because the existence of the stress interference between the fractures can increase the fracture deflection and the fluid pressure. Figure 9 shows the percentage of fluid injected into each fracture under different number of clusters. It can be found that the inter-well fractures are not only attracted to each other due to the stress interference between the fractures, but also the fluid injection amount between the wells is significantly higher than that of the outer fractures (but also injected with more fracturing fluid than the outer fractures). This is because the fracture-induced stress affects the inter-well fractures and is easier to propagate. And this phenomenon is obvious in single-stage single-cluster and single-stage multi-cluster.

| Well spacing
Well spacing settings often have a significant impact on the economics of development. The effect of well spacing on synchronous fracturing fracture propagation will be studied with three different well spacing: 300 m, 400 m, and 500 m. For the case where the well spacing is 400 m, see Section 4.2.2. Figure 10 shows the fracture geometry under different well spacing. It can be found that the shape of hydraulic fractures affected by natural fractures will be more tortuous and difficult to control, and this phenomenon will be more obvious at larger well spacing. Although fractures are easier to control at smaller well spacing, the reduction in well spacing necessarily results in increased development costs. It also can be found that the larger the well spacing is, the larger the fracture width will be when the fracture is completed. This is because the larger the spacing of the wells means the larger fluid injection time, and the length and width of the fractures increase with the injection time increase. Figure 11 shows the percentage of fluid injected into each fracture under different well spacing. It can be found that the larger the well spacing, the greater the difference in the amount of fluid injected between the fractures. This is because the larger of the well spacing, the larger of the fracture length and width will be. And the interference between the fractures will be greater with large fracture length and width. Therefore, in order to prevent the difference between the fractures from being too large, a reasonable well spacing should be selected.

| Cluster spacing
The selection of cluster spacing has always been a hot issue in the discussion of horizontal well staged fracturing. The effect of cluster spacing on synchronous fracturing fracture propagation will be studied, and three different cluster spacing: 20 m, 30 m, and 40 m. (The effect of cluster spacing on synchronous fracturing fracture propagation will be studied with three different cluster spacing: 20 m, 30 m, and 40 m). Figure 12 shows the fracture geometry under different cluster spacing. It can be found that the smaller the cluster spacing is, the larger the fracture deflection is, and even the phenomenon that two fractures gather in the same stage occurs. This is because small cluster spacing can cause great stress interference between adjacent fractures in the same stage. Therefore, when selecting the cluster spacing, it should not be too small or too large, which will cause some reservoirs to be stimulated insufficiently and excessively respectively. It also can be found that as the cluster spacing decreases, the fracture width also decreases, which is because the reduction of the cluster spacing reduces the stress interference between the fractures. Therefore, appropriately increasing the cluster spacing can reduce the construction pressure, which can not only reduce the construction risk but also reduce the construction cost. Figure 13 shows the percentage of fluid injected into each fracture under different cluster spacing. It can be found that as the cluster spacing increases, the difference in the cumulative injection volume between the individual fractures is reduced, which is due to the reduction of the stress interference caused by the increased cluster spacing. Therefore, the cluster spacing is not as small as possible, and a reasonable cluster spacing allows each fracture to be fully developed.

| Density of natural fracture
The existence of natural fractures has a significant impact on the direction of fracture propagation and the shape of fractures. Then, the cumulative effect of multiple natural fractures on artificial fractures will have a huge impact on its trajectory. The effect of well spacing on synchronous fracturing fracture propagation will be studied with three different density of natural fracture: 8.3 × 10 -3 Number/m 2 , 16.7 × 10 -3 Number/m 2 , and 25.0 × 10 -3 Number/m 2 . Figure 14 shows the fracture geometry under different well spacing. It can be found that the density of natural fractures has a significant effect on the trajectory of the fracture. The greater of the natural fracture density is, the more tortuous of the fracture trajectory will be, which means that the fracture trajectory is more difficult to predict and control. It also can be found that large width of artificial fractures can be obtained with large natural fractures density, which is caused by the more tortuous fractures when the density of natural fractures is greater. Figure 15 shows the percentage of fluid injected into each fracture under different well spacing. It can be found that the reservoir with a higher natural fracture density has a larger difference in fluid injection volume between the artificial fractures.

| CONCLUSIONS
1. In unconventional reservoirs, which contain a large number of natural fractures, the geometry of artificial fractures can be very complicated, because the propagation