Study on vertical multifracture propagations in deep shale reservoir with natural fracture network

Building a complex fracture network in a deep shale gas reservoir during multicluster hydraulic fracturing is challenging. On the one hand, hydraulic fractures tend to propagate along shale bedding planes, thus leading to a limited hydraulic fracture height; on the other hand, under the effect of high in situ stress and weak plane of fractures, multihydraulic fractures can hardly initiate and propagate uniformly. All of these will jeopardize the performance of hydraulic fracturing. This paper considers the influence of natural fracture network, shale beddings, fracturing‐fluid leakage, and high in situ stresses, and builds a numerical model to simulate multifracture propagations along vertical direction during hydraulic fracturing in deep shale reservoirs. The model is developed based on the theory of solid–fluid coupling and finite element method with pore pressure cohesive element technique are used to simulate multicluster fracture initiation and propagation. Moreover, the model is verified by published experiment data, and a good agreement between model and experiment results is observed. The influences of pump rate and fracturing‐fluid viscosity on fracture geometry are discussed, and the flow dynamics around each perforation cluster under different perforation parameters are also studied. Finally, suggestions to improve overall outcome of multicluster hydraulic fracturing are also provided.


| INTRODUCTION
Building a complex fracture network in a shale gas reservoir is a common practice to improve shale gas production. 1 The section length and cluster spacing in multicluster hydraulic fracturing keep decreasing since 2012. In Haynesville gas field, 2 the section length was between 30 and 60 m and the cluster spacing was 6-15 m in 2016. Currently, some shale gas fields deeper than 3500 m apply a section length of 50-55 m and cluster spacing of 6-15 m in southern Sichuan. 3 However, since hydraulic fractures tend to grow along shale bedding planes, the vertical propagations of hydraulic fractures would be affected, thus leading to short fracture heights and eventually a limited effect in production enhancement after hydraulic fracturing. 4,41 Therefore, improving multifracture propagations along the vertical direction is significant. In this paper, we mainly consider two effects on vertical propagations of fracture: 3,5 (i) flow dynamics in the wellbore-perforation system and (ii) natural fracture network. Many studies have been done on flow dynamics. Lecampion and Desroches 6 studied the influence of frictional resistance of the wellbore on flow distribution in perforation clusters, then further studied the effect of such flow distributions on multiple hydraulic fracture propagations. Wu and Olson 7 conducted a numerical study on multicluster fracture propagation in horizontal Wells. The results show two ways that can contribute to a uniform growth of multiple fractures: one is to adjust perforation friction resistance by changing the number or diameters of perforations. The other is to balance stress field interferences by uneven perforation spacing. Li et al. 8 considered stress field interference, wellbore friction, perforation friction, and fluid flow resistance, and built a numerical model to determine the flow distribution in different fractures during multicluster fracturing. However, the natural fractures and extensively developed bedding planes of deep shale formation tend to lead to different growths of hydraulic fractures and correspondingly more complex flow distributions. Therefore, further research is needed.
Interactions between hydraulic and natural fractures are essential for building a complex fracture network. 9 Natural fractures are well developed in shale reservoirs, and under high in situ stresses, the geometry of interaction between hydraulic fractures and bedding planes is particularly important to improve fracture propagation in the vertical direction. Wu and Olson 10 studied the interaction between hydraulic and natural fractures, and believed there are three possibilities that might occur when hydraulic fractures intersect with natural fractures. First, the hydraulic fracture could cross the natural fracture and continues to propagate in the original direction. Second, the hydraulic fracture could be arrested by the natural fracture, whereby propagation diverts along the natural fracture and then comes out from the tip of the natural fracture. Third, the hydraulic fracture deflects into the natural fracture and reinitiates at some weak point along the natural fracture. Xie et al. 11 proposed a three-dimensional (3D) hydraulic fracture propagation model, which was applied to shale gas reservoirs with bedding planes. The model results show that pump rate and treatment size are external factors that affect the interaction between hydraulic fractures and natural fractures. Yushi et al. 12 conducted a series of true triaxial fracturing experiments, natural beddingdeveloped shale has been hydrofraced. The results show that partially opened bedding planes in shale are helpful in creating complex fractures. Through laboratory tests and numerical simulation, Zhang et al. 13 and Heng et al. 14 showed that bedding planes tend to cause insufficient vertical propagation height of hydraulic fractures, which also affects the geometry of fracture propagation. Considering the leakage of fracturing fluid and the opening of natural fracture, the interaction between hydraulic fracture and natural fracture will be more complicated. Therefore, these factors need to be taken into the model.
Many numerical methods have been used to study multiple fractures propagation. The finite element method (FEM) is an effective and widely used numerical method. Haddad and Sepehrnoori 15 and Dahi-Taleghani et al. 16 established an interaction model between hydraulic fractures and natural fractures based on the cohesive element of the ABAQUS platform. To simulate hydraulic fracture propagation with the cohesive crack model, an extended finite element method (XFEM) was developed by Mohammadnejad and Khoei. 17,18 Han et al. 19 use XFEM and simulating the fracture propagation and interactions between adjacent hydraulic fractures under different injection rate increments. Zhang et al. 20,21 used the Displacement Discontinuity Method (DDM) to simulate hydraulic fractures. In the DDM method, rock is generally supposed to be an impermeable medium; therefore, the solid-fluid coupling between matrix pores and fracture network cannot be considered. Nasehi and Mortazavi 22 used the discrete element method (DEM) to simulate hydraulic fracturing, but DEM needs to obtain approximate mechanical parameters through trial calculation, so the accuracy of analysis results is affected. The finite-discrete element method was proposed by Munjiza et al. 23,24 It has also been applied to simulate complex fracture propagation in 2D and 3D models. 25,26 But this method assumes that only the fluid flows in the fracture elements and no fluid leakage is considered in the simulation. Huang et al. 27 used finite element and cohesive zone methods, studying the effects of stress field, natural fracture properties, and approach angle on the geometry of interaction between natural and hydraulic fractures. In conclusion, FEM with cohesive element technique has better advantages to study the propagation of multicluster fractures in some aspects.
This paper built a model to simulate multifracture propagations in the vertical direction during multicluster hydraulic fracturing. A solid-fluid coupling in shale matrix is considered by the theory of poroelasticity, and FEM with cohesive element technique is employed to simulate fracture initiations and growths. In addition, the model also considers the initial opening of natural fractures, the fluid leakage into the shale matrix and fluid flow dynamics. Finally, pump rate, fracturing-fluid viscosity, perforation diameter, and perforation number are optimized based on simulation results. The study can help improve the overall performance of hydraulic fracturing.  2.1 | Flow in the wellbore-perforation system Figure 1 shows the connections between wellbore elements, perforation elements, and cohesive elements. The injected fracturing fluid flows along the wellbore element and enters the formation through perforation elements and cohesive elements. As shown in Figure 2, the perforation entry element consists of two nodes. Pore pressure degree of freedom is assigned to each node.
The fluid flow in wellbore and perforation elements is considered as a 1D flow. The relation between total injection rate Q 0 and flow rate into each perforation is where Q i is the flow rate into the perforation i, i = 1-3.
The relation between flow rate and pressure in a wellbore element is determined by Equation (2): 42 where ∆p is the pressure difference between the two nodes, h Δ is the elevation difference between two nodes, Q 0 is the fluid volume flow in the pipe, ρ is the fluid density, g is the gravity acceleration, Z is the friction loss coefficient, L is the pipe length, f is the friction factor, A is the cross-sectional area of the pipe, and S is the wetted perimeter of the pipe.
For perforation elements, we assume all perforation elements share the same pressure at connecting to wellbore elements, that is, And the relation between flow rate and pressure in a perforation element is determined by Equation (5): 28,29 where i is the number of perforation clusters, p i w is the fluid pressure of the fracture at the node where the fluid flows into the perforation element, p i m is the fluid pressure of the fracture at the node where the fluid flows out the perforation element, p Δ i fric is the pressure drop of perforation cluster i, and ρ is the fluid density. n p is the perforation number of cluster i, which is in the range of 6-30; D p is the perforation diameter typically ranging between 6 and 25 mm. C is a dimensionless discharge coefficient reflecting the shape change of the perforation hole due to erosion, which is typically between 0.5 (before erosion) and 0.9 (after erosion).
By the above equations, the flow rate and pressure drop at cluster i can be determined. In a word, fracturing fluid needs to overcome the resistance in wellbore elements and perforation elements to flow into the formation. For different perforation parameters, the resistance at different perforation clusters is different, eventually leading to different flow rates and fracture growths at different clusters.

| Modeling of fracture initiation and propagation by cohesive element
The bilinear cohesive model is applied to simulate the initiations and propagations of hydraulic fractures. As shown in Figure 3, the model assumes a process zone at the crack tip, and the fracture in the process zone is described by the traction-separation (T-S) constitutive model, which can avoid the stress singularity at the crack tip.
The T-S constitutive model exhibits linear elastic behavior when applied loads are in the elastic range. Once the applied loads are beyond the elastic range, the damage evolves. The initial linear elastic relation of T-S constitutive model can be expressed as 15,30 where k is the stiffness matrix of a cohesive element in elastic response. The initiation of fractures is modeled by a sequence of damages of cohesive elements, and the damage of a cohesive element is determined by the quadratic normal stress criterion, 31 that is, where t n , t s1 , and t s2 are the nominal tractions along the normal and tangent directions, respectively; t n 0 , t s1 0 , and t s2 0 are the corresponding prescribed damage initiation tractions;  is the Macaulay bracket; In 2D space, t s2 and t s2 0 disappear.
Once the corresponding initiation criterion is reached, the damage evolution model describes the rate at which the material stiffness declines. 32 The damage evolution formulae are as follows: where tn , t¯s , and t¯t are the stresses in the three loading directions with the assumption that the evolution does not happen and the model is still in linear elastic deformation process; t n , t s , and t t are the real stresses in the three loading directions; D is dimensionless damage coefficient ranges from 0 to 1. If D = 0 the material has not damaged yet and if D = 1 the material damaged completely. Dimensionless damage coefficient D as follows: 15 where δ m 0 and δ f m are the initial damage displacement and the complete damage displacement of the cohesive element, respectively; δ m is the equivalent displacement of the damage evolution stage controlled by the normal displacement and the tangential displacement, and the expression as follows: 32 where δ n , δ s , and δ t are the normal, first, and second tangential damage displacements of the cohesive element during the damage, respectively. Under the condition of plane strain, the critical fracture energy G n of rock normal (type I fracture) can be determined by the critical maximum stress T max that the rock can bear 14,34 where E is the elastic modulus, ν is Poisson's ratio, and K n is the type I fracture toughness. In the same way, | 2895 the critical fracture energy corresponding to type II and III fractures can be calculated. After the damage occurs, the damaged surface will undergo tensile or shear failure, as shown in Figure 3A. To consider the combination of normal and tangential extension. Benzeggagh-Kenane (BK) mixed propagation mode is adopted as fracture propagation criterion: 35 where G n c is the critical fracture energy in type I; G s c is the critical fracture energy in type II; G n is the accrued fracture energy in mode I; G s is the accrued fracture energy in type II; G t is the accrued fracture energy in type III; η is the exponent, dimensionless; G c is the total critical energy release rate of the cohesive element in the mixed mode. The actual total energy release rate

| Flow in fractures
Momentum equation for incompressible Newtonian fluid in fractures is 36 And the corresponding continuity equation is given by 35 where q is the fluid flux of flow in fracture, p is flow pressure in the fracture, w represents the fracture width, and μ is the fracturing-fluid viscosity. The total flux of flow leaking into shale matrix q tot is determined by 36 where q t and q b are the flow flux normal to the top and bottom fracture surfaces (i.e., the interfaces between a cohesive element and its neighboring shale matrix element), and are determined by where p i is the pressure at the center of a cohesive element; p t and p b are the pressure on the top and bottom surfaces of a cohesive element, respectively.

| Flow in shale matrix
The constitutive model of shale matrix is 37 where σ ij is the stress; λ and G are the Lamé parameters of the porous material. While C and M are poroelastic moduli required to describe a two-phase medium. δ ν is the volumetric strain; δ ij is the Kronecker operator; ξ is a strain parameter describing the volumetric deformation of the fluid relative to that of the solid; p is the pore pressure.
The stress of rock matrix should satisfy the stress equilibrium equation Pore fluid flow should satisfy mass conservation equation where q i is the fluid flux.

| HYDRAULIC FRACTURES PROPAGATIONS BY COHESIVE ELEMENT TECHNIQUE
To simulate rock dynamic fracture without presetting fracture propagation path, we have established a method of zero-thickness pore pressure cohesive element (PCE) inserting into the global finite element (PCE-FEM). PCE-FEM divides the target area into a finite number of matrix elements. The matrix element is used to solve the stress and strain of the rock, and the pore PCE is used to solve the solid-fluid coupling of hydraulic fractures. The method is mainly divided into four stages: (1) Using an unstructured advanced algorithm divides the target area with a quadrilateral grid, as shown in Figure 4A. The unstructured mesh is easy to control the size of the mesh and the density of nodes, at the same time, it can divide irregular elements, which can realize the effect of hydraulic fractures diversion and propagation between matrix elements.
(2) Circulating all elements in the target area, according to the number of parent node occurrences N in these elements. Then copying the coordinate information of parent nodes generate N − 1 child nodes with the same coordinate information. (3) After parent nodes split into child nodes, a simple topology data structure is used to save the node number and its coordinate information. Thus, this forms child node information consisting of the number of repetitions and coordinates of the parent node. (4) Constructing rock matrix elements and cohesive elements based on the new node information. The numbers of parent node and child node are reconstituted to form rock matrix elements, and their types are unchanged as shown in Figure 4C. To ensure that fluid flow is allowed between elements, a node with pore pressure degree of freedom is inserted in the middle of each adjacent element boundary. Then the pore pressure nodes with the same coordinates are merged into a common pore pressure node to make fluid continuity in the element. Taking the two solid nodes on the boundary as the starting point, combined with the pore pressure node, the cohesive element is constructed in a counterclockwise order ( Figure 4B), as shown in Figure 4D. In this way, cohesive elements are inserted into the global finite elements.

| MODEL VERIFICATION
We verified the feasibility of PCE-FEM to study the propagation of the interaction between hydraulic and natural fractures, based on Blanton's laboratory tests data, 38 as shown in Table 1. It is assumed that the value of each property of natural fractures is 10% of that of hydraulic fractures. The growth of hydraulic fracture and its interaction with a natural fracture are simulated. Figures 5 and 6 show the comparison between simulation results and Blanton's tests CT-7 and CT-8, respectively. In Figure 5, test CT-7 has a horizontal stress difference of 9 MPa. In this case, the pressure at the intersection of hydraulic fracture and natural fracture is greater than the normal stress on its fracture surface, which overcomes the normal stress on natural fracture and leads to the opening of natural fracture. In test CT-8, the horizontal stress difference is increased to 15 MPa. The larger difference in horizontal stress increases the difficulty of stress field interference, and the result is that hydraulic fracture crosses the natural fracture and continues to grow along the original direction. In addition, the vertical stress of the test CT-14 is 20 MPa, its maximum horizontal stress is 14 MPa and the minimum horizontal stress is 5 MPa, and the injection rate is 0.82 cm 3 /s. According to Figure 7, it can be seen that the fracture pressure and extension pressure of the test CT-14 are in good agreement with the simulation. In sum, the simulation results and experiment are in good agreement. Thus, it is feasible to use the PCE-FEM method to simulate the propagation of the interaction between hydraulic fractures and natural fractures.  Table 2.
In Figure 8, a 2D shale reservoir natural fracture model is established through statistical data. The length and the width of the model both are 80 m.

| Multifracture propagations in vertical direction
The horizontal well W204H34-5 in southern Sichuan is taken as an example to analyze. The completion horizon is Longmaxi Formation, and the horizontal section of this well is 3370-5501.1 m. There are 15 fracture development sections at the identification of fracturing section. The horizontal differential stress is 14.5 MPa, and the  maximum principal stress direction is at an angle of 70°-90°with the wellbore. On the one hand, the fracture width is narrowed, which hinders the migration of proppant in the fractures. On the other hand, the fluid leakage is very large, leading to hydraulic fractures propagating in the vertical direction insufficiently. In this paper, by shortening the distance of perforation clusters, three clusters are adopted to improve the performance of hydraulic fractures propagating in the vertical direction. The model size is the same as the established natural fracture model, and the cluster spacing is 15 m.
(1) Mesh model In this paper, the numerical simulation is conducted on ABAQUS platform. The solid-fluid coupling pore pressure plane strain element (CPE4P) with an average grid size of 1 m is used to divide the model, and cohesive elements were inserted into the finite elements, and its type is six-node displacement and pore PCE (COH2D4P). An independent set of cohesive elements of natural fractures (Figure 9) is established to give attribute parameters of natural fractures.
(2) Load conditions The vertical stress is applied in the Y-direction, and the minimum horizontal principal stress is applied in the X-direction. The boundary with wellbore is set as a symmetrical boundary, and the other boundaries are set as a fixed displacement (3) Model parameters Parameters of the rock matrix in the model are as follows in Table 3. Table 4 shows the parameters of ordinary cohesive element and natural fracture cohesive element. Table 5 is the parameters for fracturing and perforation parameters.

| RESULTS AND DISCUSSION
First, the influence of pump rate and viscosity on the geometric shape of multifracture propagation is analyzed. Then we discuss the influence of perforation parameters on the flow dynamic distribution of each fracture during multifractures propagation. The simulation experiment scheme is shown in Table 6.

| Influence of pump rate
This section studies the influence of pump rate on the vertical propagation of fractures. The viscosity, perforation diameter, and perforation number are the same as those in Table 5, while the pump rate is changed from 8 to 20 m 3 /min. As shown in Figure 10, the simulation results are amplified by 100 times. According to the four solid nodes of the cohesive element, the coordinates of the two solid nodes are used to obtain a midpoint coordinate, and then using the two midpoint coordinates obtain the length of cohesive element. Finally, determining whether the element has failed. If it fails, the length can be accumulated. In this way, the length of hydraulic fractures is obtained, and the calculation result is shown in Figure 11.
According to Figures 10 and 11, when the pump rate is small, the fracturing fluid mainly flows along with natural fractures and bedding planes, which is beneficial to open the natural fracture system. Thereby small pump rate increases the complexity of hydraulic fractures, but it is difficult to propagate in the vertical direction ( Figure 10A). As the pump rate is bigger, enhancing the ability of hydraulic fractures to cross bedding ( Figure 10D). However, it is not conducive to form a complex fracture network.  We can see from Figure 11 that the bigger pump rate increases the net pressure in the fracture, which makes bigger average width of the hydraulic fractures. Moreover, a bigger fracture width is helpful to solve the difficulty of adding proppant into the fractures. In Figure 12, the rate of hydraulic fracture initiation is positively related to the pump rate. Larger pump rate increases fracture initiation rates. When the pump rate is 8 m 3 /min, affected by the stress field interference, the middle hydraulic fracture cannot propagate, so the total length of the hydraulic fracture is the smallest. In other cases, the smaller pump rate makes the total length of hydraulic fractures larger. The reason is that it is easier for a small pump rate to help hydraulic fractures to propagate in the natural fractures system. Thereby resulting in a larger total length of hydraulic fractures. When the pump rate is 20 m 3 /min, hydraulic fractures directly across the bedding planes, which cannot build a complex fracture network and thus result in a poor fracturing effect. Considering the total length of the fractures and the complexity of the hydraulic fracture network, when the pump rate is 16 m 3 /min, the reservoir has a better improvement effect.

| Influence of viscosity
From the above results, the pump rate is set to 16 m 3 / min, and the fracturing-fluid viscosity is changed from 3 to 30 mPa⋅s, while the perforation diameter and perforation number are the same as those in Table 5. Simulation results are shown in Figures 13  and 14. Low viscosity fracturing fluid can better transmit pressure and communicate with more natural fractures. Therefore, it can build a complex fracture network, but has a poor vertical propagation of hydraulic fractures. In the case of high viscosity fracturing fluid have less leakage, which will significantly increase the net pressure in the fracture. It is beneficial for the hydraulic fracture to cross through bedding planes. According to Figures 13 and 14, the higher viscosity increases the width and length of the fractures. Therefore, in deep shale fracturing, in the early stage of pumping, high viscosity fracturing fluid is used to create main fractures, while in the later stage, low viscosity fracturing fluid is used to create branch fractures to communicate more natural fractures at the far end of the wellbore, and improving the stimulated reservoir volume. It is recommended to use fracturing fluid with a viscosity of 30 mPa⋅s in the early stage of pumping, and a fracturing fluid with a viscosity of 3 mPa⋅s in the later stage of pumping.

| Influence of perforation diameter
The pump rate is 16 m 3 /min, and the viscosity is 30 mPa⋅s. According to Equation (6), the diameter of perforations is changed from 6 to 18 mm. The perforation number is the same as those in Table 5. When three clusters of fractures propagate simultaneously, the height and the width of the middle fracture are smaller than the outer fractures, because it is compressed by the outer fractures. By reducing the perforation diameter, the compression effect of the middle fracture is weakened, and the height of each cluster of fractures in the vertical direction is closer ( Figure 15A). It can be seen from Figure 16 that a smaller perforation diameter can promote a more uniform flow distribution among the fractures. When the perforation diameter is 6 mm, the flow distribution between the fractures is more uniform. Equation (6) shows that reducing the perforation diameter will increase the perforation friction, which weakens the stress field interference between the clusters of fractures and promotes the uniform growth of each fracture. However, under the same injection conditions, the smaller perforation diameter leads a higher perforation friction. Eventually, increasing the formation fracture pressure and the difficulty of hydraulic fracturing. Therefore, when optimizing the perforation diameter to promote the uniform flow of fluid into each fracture, it is necessary to consider the situation of a negative effect of decreasing the perforation diameter will increase the formation fracture pressure. On the basis of this reason, suggesting perforation diameter is 10 mm.

| Influence of perforation number
The pump rate is 16 m 3 /min, viscosity is 30 mPa⋅s, and perforation diameter is 10 mm. According to Equation (6), the number of perforations is changed from 10 to 22 holes. Analyzing the influence of different perforation numbers on the flow distribution of each cluster of fractures. According to Figure 17A, it can be seen that the smaller perforations number makes the height of each hydraulic fracture closer in the vertical direction. In Figure 18, reducing the perforation number can promote the more uniform fluid inflow of the fractures. When the perforation number is 10 holes, the flow distribution among the fractures is the most uniform ( Figure 18A). Equation (6) shows that reducing the perforation number can increase the perforation friction, when the fluid passes through the perforation, the resistance will increase. It will reduce the difference in fluid pressure in each cluster of fractures. Thereby, promoting the uniform propagation of each cluster of fractures. However, the low number of perforations will also bring about high perforation friction and fracture pressure. Comprehensively, it is recommended that the number of perforations per cluster is 14 holes. In this paper, we build a solid-fluid coupling model to simulate the complex hydraulic fractures propagations along a vertical direction in deep shale reservoirs with natural fracture networks using PCE-FEM technique. In addition, fluid flow in fractures, fluid leakage, and opening of natural fractures can be considered in this model. By comparing Blanton's experimental results with the corresponding numerical simulation results, it is verified that this method can accurately simulate the propagation behavior of the interaction between hydraulic fractures and natural fractures. The influences of pump rate and fracturing-fluid viscosity on fracture geometry are discussed, and the flow dynamics around each perforation cluster under different perforation parameters are also studied. It provides a reference for deep shale gas fracturing exploitation in southern Sichuan. The main conclusions are as follows: (1) For the results of numerical simulation, when the pump rate is 8 or 12 m 3 /min, a low pump rate helps open the natural fracture system, increasing the complexity of the hydraulic fractures. When the pump rate is 20 m 3 /min, hydraulic fractures propagate directly across the bedding planes, which cannot build a complex fracture network resulting in poor a fracturing effect. Comprehensively, when the pump rate is 16 m 3 /min, and the reservoir has a better improvement effect. (2) Low viscosity forms a complex fracture network, while high viscosity is beneficial for the hydraulic fracture to cross through bedding planes. In the early stage of pumping, high viscosity fracturing fluid is used to create main fractures, while in the later stage, low viscosity fracturing fluid is used to create branch fractures. So it is recommended to use fracturing fluid with a viscosity of 30 mPa⋅s in the early stage of pumping, and a fracturing fluid with a viscosity of 3 mPa⋅s in the later stage of pumping. (3) By reducing the perforation number and the perforation diameter, the perforation friction will increase, which will weaken the interference of the stress field, and thus contribute to a uniform fluid inflow and uniform propagation of fractures at each cluster. However, excessively reducing the number of perforations and the diameter of perforations will increase the perforation friction and the fracture pressure. These factors should be considered when optimizing the number and diameter of perforation. Therefore, suggesting perforation diameter is 10 mm, and perforation number per cluster is 14 holes.