Perforation Optimization of Intensive-Stage Fracturing in a Horizontal Well Using a Coupled 3D-DDM Fracture Model

: Intensive-stage fracturing in horizontal wells is a potentially new technology for reservoir stimulations of deep shale oil and gas. Due to a strong stress interaction among the dense fractures, the fracture geometry and stress ﬁeld are very complicated, which are the bottlenecks of this technology. Aiming at simulating the intensive-stage fracturing, a coupled three-dimensional (3D) fracture model of multiple-fracture simultaneous propagation is proposed. The dynamic behavior of the fracture propagation and stress ﬁeld was analyzed using this model. The perforation parameters were optimized for improving the fracture geometry equilibrium. The results showed that the exterior fractures of the multiple fractures penetrated by the horizontal well become the main fractures, while the interior fractures are drastically restrained. The exterior fracture widths increased with increasing injection time, while the interior fracture widths decreased with increasing injection time. An extruded region was created among the multiple fractures, which restrained the propagation of the interior fractures. Only increasing the perforation cluster number did not improve the fracture geometry equilibrium in the intensive-stage fracturing. To improve the fracture geometry equilibrium, we suggest designing more perforation numbers in each perforation cluster and ensuring that both the perforation number and diameter in the interior perforation cluster are greater than those of the exterior ones.


Introduction
Unconventional oil and gas developments have been the most rapidly expanding trend due to a huge demand for fossil energy. The development of deep shale oil and gas has been an important area over the past few years. Hydraulic fracturing technology is widely used in enhancing shale oil and gas productions. Due to the large in-situ stress field, high clay content and low fracability in the deep area, hydraulic fracturing technology faces bottlenecks. Aiming at deep shale oil and gas reservoir simulation, field engineers try to reduce the length of each fracturing section and increase the perforation cluster number compared with those of conventional multi-stage fracturing technology. This new technology is called intensive-stage fracturing [1][2][3]. The horizontal wellbore length in each fracturing section decreases from 80-100 m to 20-25 m, while the perforation cluster spacing decreases from 30-50 m to 5-10 m [1,2]. However, the efficiency of intensive-stage fracturing technology urgently needs to be improved. Previous studies indicated that some perforation clusters could not create high-conductivity fractures or even initiate fractures [2,4,5]. This phenomenon increases the perforation cost and wastes the reservoir volume. Each fracturing section has many perforation clusters in the multi-stage fracturing technology. The stress interaction among the fractures created from the perforation cluster predominantly affects the fracture configuration [6][7][8]. Thus, building a fracture model for In intensive-stage fracturing of a horizontal well, each fracture has its own injection rate Qk at the entry point, and the total injection rate is Qc. The fluid partitioning is [5]: where Nc is the perforation cluster number. Qk is the injection rate into the kth fracture, m 3 /s. Qc is the injection rate into the wellbore, m 3 /s. Due to the flow resistance inside the fracture and the pressure drop in the perforation cluster, the fluid pressure in the wellbore is: , , , + 1, 2,3,..., w k r k in k c P P P k N = = (2) where Pr,k is the pressure drop of the flow through the kth perforation cluster, Pa. Pin,k is the flow resistance along the fracture, Pa. Pw,k is the pressure of the kth fracture at the wellbore inlet, Pa.
The perforation cluster is very closed to each other, and the frictional resistance of flowing through the casing wall is much smaller than the perforation hole friction such that the pressure is equal at every entry point [5,36]: where Pw is the fluid pressure in the wellbore inlet, Pa.
To build a 3D fracture model for intensive fracture propagation, the main assumptions are as follows: ① The fracture opening and closure satisfy the linear elasticity; ② The fracture tip propagation satisfies the linear elastic fracture mechanics such that the maximum tensile stress strength criterion is selected for judging the fracture initiation and termination; ③ The fracture geometry is a planar 3D plane. Hence, a rectangular mesh can be used to discretize the fracture plane; ④ The leak-off of fracturing fluid satisfies Carter's leak-off equation because of low permeability; ⑤ The flow resistance along the wellbore is ignored because it is much smaller than the perforation friction; ⑥ The horizontal well is assigned along the minimum horizontal principal stress, which is the general design.

Fracture Deformation
In general, the hydraulic fracture propagates along a plane perpendicular to the minimum principal stress. This plane is discretized into many rectangular meshes ( Figure 2). With the injection of fracturing fluid, the meshes are activated sequentially. All the activated meshes constitute the shape of hydraulic fracture (Figure 3). In intensive-stage fracturing of a horizontal well, each fracture has its own injection rate Q k at the entry point, and the total injection rate is Q c . The fluid partitioning is [5]: where N c is the perforation cluster number. Q k is the injection rate into the kth fracture, m 3 /s. Q c is the injection rate into the wellbore, m 3 /s. Due to the flow resistance inside the fracture and the pressure drop in the perforation cluster, the fluid pressure in the wellbore is: where P r , k is the pressure drop of the flow through the kth perforation cluster, Pa. P in , k is the flow resistance along the fracture, Pa. P w,k is the pressure of the kth fracture at the wellbore inlet, Pa. The perforation cluster is very closed to each other, and the frictional resistance of flowing through the casing wall is much smaller than the perforation hole friction such that the pressure is equal at every entry point [5,36]: where P w is the fluid pressure in the wellbore inlet, Pa.
To build a 3D fracture model for intensive fracture propagation, the main assumptions are as follows: 1 The fracture opening and closure satisfy the linear elasticity; 2 The fracture tip propagation satisfies the linear elastic fracture mechanics such that the maximum tensile stress strength criterion is selected for judging the fracture initiation and termination; 3 The fracture geometry is a planar 3D plane. Hence, a rectangular mesh can be used to discretize the fracture plane; 4 The leak-off of fracturing fluid satisfies Carter's leak-off equation because of low permeability; 5 The flow resistance along the wellbore is ignored because it is much smaller than the perforation friction; 6 The horizontal well is assigned along the minimum horizontal principal stress, which is the general design.

Fracture Deformation
In general, the hydraulic fracture propagates along a plane perpendicular to the minimum principal stress. This plane is discretized into many rectangular meshes ( Figure 2). With the injection of fracturing fluid, the meshes are activated sequentially. All the activated meshes constitute the shape of hydraulic fracture ( Figure 3).   A hydraulic fracture is treated as a discontinuous displacement plane. The displacement discontinuity method (DDM) provides the displacement and stress solutions in an infinite elastic medium [12,37]. The 3D DDM is based on the elastic solution to the problem of a discontinuous displacement over a finite area in an infinite region [38]. The fundamental solutions of 3D DDM are [24,39]: ikl U κ ζ are the influential coefficients of the kth displacement discontinuity of the lth fracture element on the ijth stress and the ith displacement at the location κ , respectively. kl D ζ ( ) is the discontinuous displacement at the location ζ .
Γ represents the fracture surface.
All of the hydraulic fractures have N rectangular fracture elements in total, and each element has three discontinuous displacement components Dlk (k = 1, 2, 3; l = 1, 2, 3, …, N) in the local coordinate system originated the fracture element center. The discrete form of Equation (4) is [39,40]: where l ζ is the center of the lth fracture element.  A hydraulic fracture is treated as a discontinuous displacement plane. The displacement discontinuity method (DDM) provides the displacement and stress solutions in an infinite elastic medium [12,37]. The 3D DDM is based on the elastic solution to the problem of a discontinuous displacement over a finite area in an infinite region [38]. The fundamental solutions of 3D DDM are [24,39]: where ( ) Γ represents the fracture surface.
All of the hydraulic fractures have N rectangular fracture elements in total, and each element has three discontinuous displacement components Dlk (k = 1, 2, 3; l = 1, 2, 3, …, N) in the local coordinate system originated the fracture element center. The discrete form of Equation (4) is [39,40]: where l ζ is the center of the lth fracture element. A hydraulic fracture is treated as a discontinuous displacement plane. The displacement discontinuity method (DDM) provides the displacement and stress solutions in an infinite elastic medium [12,37]. The 3D DDM is based on the elastic solution to the problem of a discontinuous displacement over a finite area in an infinite region [38]. The fundamental solutions of 3D DDM are [24,39]: where Π ijkl (κ, ζ) and U ikl (κ, ζ) are the influential coefficients of the kth displacement discontinuity of the lth fracture element on the ijth stress and the ith displacement at the location κ, respectively. D kl (ζ) is the discontinuous displacement at the location ζ. Γ represents the fracture surface. All of the hydraulic fractures have N rectangular fracture elements in total, and each element has three discontinuous displacement components D lk (k = 1, 2, 3; l = 1, 2, 3, . . . , N) in the local coordinate system originated the fracture element center. The discrete form of Equation (4) is [39,40]: where ζ l is the center of the lth fracture element.
A hydraulic fracture suffers both local and far-field stresses. The discontinuous shear displacement is assumed to be zero in this model. The physical meaning of the normal discontinuous displacement is the fracture aperture. According to Equation (5), the hydraulic fracture opening satisfies [41]: where p(κ) is the fluid pressure acting at the point κ of the fracture surface, Pa. S h is the minimum horizontal stress, Pa. w l is the fracture aperture of the lth fracture element, m. C(κ, ζ l ) is the influential coefficient of the lth fracture opening on the normal stress at the fracture surface, Pa/m.

Flow of Fracturing Fluid
The friction pressure loss of a single perforation hole is [8,31]: where ρ is the density of fracturing fluid, kg/m 3 . C k is the dimensionless perforation coefficient. n k is the perforation number of the kth perforation cluster. d k is the perforation diameter of the kth perforation hole, m. The fluid equation inside the fracture is [42]: where w is the fracture width, m. q l is the leak-off velocity of the fracturing fluid, m/s. u is the viscosity of the fracturing fluid, Pa·s. Fracturing fluid leak-off velocity is [43]: where C lv is the leak-off coefficient, m/s 0.5 . t 0 is the moment when the fracture opens, s. Using the finite difference method, Equation (8) is divided as: If the square mesh (∆x = ∆y) and central difference are used, Equation (10) is simplified as:

Fracture Propagation Condition
The maximum tensile stress [24] at the type-I fracture tip is: where K I is the stress intensity factor of the type-I fracture tip, Pa·m 0.5 ; r, θ represent the local coordinate system originating at the fracture tip. σ θ is the tensile stress at the fracture tip, Pa. According to the maximum tensile strength criterion [24], the fracture will propagate when K I > K IC : where K IC is the fracture toughness of type-I fracture, Pa·m 0.5 .

Fluid-Solid Coupling Method
The combination of Sections 2.1-2.4 constitutes a fluid-solid coupled 3D fracture propagation model for the simulation of intensive-stage fracturing in the horizontal well. Newton-Raphson iteration method is used to solve it. The residual function vector is: where The solution vector is: The iteration format is: The iteration convergence is checked by: The analytical solution of a penny-shaped fracture is used as x (0) . Once the iteration reaches convergence, the fracture element widths are used to judge the fracture tip propagation. New fracture elements are added to the fracture tip once the previous fracture tip propagates. The fracture footprint is built step-by-step (Figure 4).

Benchmarking
In the absence of fracture fluid leak-off, the analytical solutions of the fracture radius and the central fracture width of a penny-shaped fracture are [24,27]: where E is Young's modulus, Pa. t is the injection time, s. W 0 (t) is the centric width of hydraulic fracture, m. v is the Poisson's ratio, dimensionless. R(t) is the radius of hydraulic fracture, m. A comparison of the analytical and numerical solutions is given in Figure 5, using the input parameters in Table 1.

Benchmarking
In the absence of fracture fluid leak-off, the analytical solutions of the fracture radius and the central fracture width of a penny-shaped fracture are [24,27]: where E is Young's modulus, Pa. t is the injection time, s. W0 (t) is the centric width of hydraulic fracture, m. v is the Poisson's ratio, dimensionless. R(t) is the radius of hydraulic fracture, m. A comparison of the analytical and numerical solutions is given in Figure 5, using the input parameters in Table 1.

Fracture Geometry and Stress Field Analysis
The traditional perforation cluster spacing in the multi-stage fracturing technology is about 30-50 m, and the number of perforation clusters is about 2-3. The perforation cluster spacing declines by 5-10 m, and the number of perforation clusters increases by 5-6,

Fracture Geometry and Stress Field Analysis
The traditional perforation cluster spacing in the multi-stage fracturing technology is about 30-50 m, and the number of perforation clusters is about 2-3. The perforation cluster spacing declines by 5-10 m, and the number of perforation clusters increases by 5-6, which is called intensive-stage fracturing. The following cases will aim at simulating the fracture propagation and stress field in the intensive-stage fracturing technology.

Fracture Geometry
The related parameters in a deep shale oil reservoir of the Qianjiang formation, Jianghan Basin are given in Table 2. The hydraulic fracture geometry at an injection time of 5 min is given in Figure 6. The exterior fractures (fracs 1 and 5) become the main fractures with the largest widths. The interior fractures (fracs 2-4) are hardly restrained and have the smallest widths. The exterior fractures' radii are much larger than those of the interior fractures because the interior fractures stop after a short time. The hydraulic fracture geometry is symmetrical along the x and y-axes. Wu [30] and Gunaydin et al. [4] observed multiple transverse fracture propagations in the PMMA sample, respectively. The fracture geometries in their test are very similar to our simulation results.
The time variation of the hydraulic fracturing is given in Figure 7. Due to the symmetry, the exterior fractures' widths (fracs 1 and 5) have the same variation versus the injection time. With increasing injection time, the fracture width increases drastically in the beginning and then gradually increases. The widths of the interior fractures (fracs 2-4) decline with the injection time. In the beginning, the declining rate is much greater than that later. Due to the symmetry, fracs 2 and 4 have the same width variation. Their widths both decline to zero because of the strong stress interaction induced by the adjacent main fractures (fracs 1 and 5). Frac 3 is located in the middle area, which is far from the main fractures; hence, the stress interaction strength is weaker, and the width of frac 3 is wider than that of fracs 2 and 4.
The fluid pressure distribution inside the hydraulic fracture is given in Figure 8. The inlet pressure is the largest, while the pressure at the fracture tip is the lowest. The fluid pressure declines along the fracture radius, which indicates a penny-shaped fracture geometry. The injection pressure at the wellbore versus the injection time is given in Figure 9. The injection pressure declines rapidly in the beginning and then gradually declines later. The final pressure tends to the original minimum horizontal stress. 4) decline with the injection time. In the beginning, the declining rate is much greater than that later. Due to the symmetry, fracs 2 and 4 have the same width variation. Their widths both decline to zero because of the strong stress interaction induced by the adjacent main fractures (fracs 1 and 5). Frac 3 is located in the middle area, which is far from the main fractures; hence, the stress interaction strength is weaker, and the width of frac 3 is wider than that of fracs 2 and 4.  The fluid pressure distribution inside the hydraulic fracture is given in Figure 8 inlet pressure is the largest, while the pressure at the fracture tip is the lowest. The pressure declines along the fracture radius, which indicates a penny-shaped fractur ometry. The injection pressure at the wellbore versus the injection time is given in Fi 9. The injection pressure declines rapidly in the beginning and then gradually dec later. The final pressure tends to the original minimum horizontal stress.

Stress Field
The stress distribution of S zz is symmetrical along a plane z = 10 m (Figure 10). The S zz in the middle area (red area in Figure 10) is much greater than the original minimum horizontal stress (S h = 58.5 MPa). This indicates a strong compression along the z-axis in the middle area. In the area far from the horizontal well (blue area in Figure 10), S zz is less than the original minimum horizontal stress (S h = 58.5 MPa). This indicates a tension along the z-axis induced by the opening fracture. Due to the stress interaction, fracs 2-4 suffer strong compressive stress created by the exterior fractures, and their propagation is strongly restrained. The fluid pressure distribution inside the hydraulic fracture is given in Figure 8. The inlet pressure is the largest, while the pressure at the fracture tip is the lowest. The fluid pressure declines along the fracture radius, which indicates a penny-shaped fracture geometry. The injection pressure at the wellbore versus the injection time is given in Figure  9. The injection pressure declines rapidly in the beginning and then gradually declines later. The final pressure tends to the original minimum horizontal stress.

Stress Field
The stress distribution of Szz is symmetrical along a plane z = 10 m ( Figure 10). The Szz in the middle area (red area in Figure 10) is much greater than the original minimum horizontal stress (Sh = 58.5 MPa). This indicates a strong compression along the z-axis in the middle area. In the area far from the horizontal well (blue area in Figure 10), Szz is less than the original minimum horizontal stress (Sh = 58.5 MPa). This indicates a tension along the z-axis induced by the opening fracture. Due to the stress interaction, fracs 2-4 suffer strong compressive stress created by the exterior fractures, and their propagation is strongly restrained.

Stress Field
The stress distribution of Szz is symmetrical along a plane z = 10 m ( Figure 10). The Szz in the middle area (red area in Figure 10) is much greater than the original minimum horizontal stress (Sh = 58.5 MPa). This indicates a strong compression along the z-axis in the middle area. In the area far from the horizontal well (blue area in Figure 10), Szz is less than the original minimum horizontal stress (Sh = 58.5 MPa). This indicates a tension along the z-axis induced by the opening fracture. Due to the stress interaction, fracs 2-4 suffer strong compressive stress created by the exterior fractures, and their propagation is strongly restrained. The stress shadowing along the z-axis induced by the opening fracture is Due to the symmetry of the x-axis and y-axis, the stress shadowing is also symmet- The stress shadowing along the z-axis induced by the opening fracture is Due to the symmetry of the x-axis and y-axis, the stress shadowing is also symmetrical along the x-axis and y-axis. The ∆S zz on 1/4 of the x-y plane is given in Figure 11. S zz is also symmetrically distributed at a plane z = 10, so the ∆S zz on the planes z = −2.5 m, 2.5 m, 7.5 m are given in Figure 11a-c, respectively. ∆S zz on the planes z = −2.5 m, z = 2.5 m, z = 7.5 m are equal to that on the planes z = 22.5 m, z = 17.5 m, z = 12.5 m, respectively. ∆S zz is greater than zero in the area adjacent to the wellbore (locating on x, y = 0), while it is less than zero in the far-field area. The stress fluctuation in the area adjacent to the wellbore is caused by the stress singularity of the rectangular element in DDM. When the region is far from the hydraulic fracture surface, the stress fluctuation will disappear. The stress distribution of Sxx is symmetrical along the plane z = 10 m, as shown in Figure 12. Sxx in the middle area (red area in Figure 12) is greater than maximum horizontal stress (SH = 62.5 MPa). This indicates compressive behavior along the x-axis. In the area far from the wellbore (blue area), Sxx is slightly less than SH. This indicates tensile behavior along the x-axis in the fracture tip area. The stress distribution of S xx is symmetrical along the plane z = 10 m, as shown in Figure 12. S xx in the middle area (red area in Figure 12) is greater than maximum horizontal stress (S H = 62.5 MPa). This indicates compressive behavior along the x-axis. In the area far from the wellbore (blue area), S xx is slightly less than S H . This indicates tensile behavior along the x-axis in the fracture tip area. The stress distribution of Sxx is symmetrical along the plane z = 10 m, as shown in Figure 12. Sxx in the middle area (red area in Figure 12) is greater than maximum horizontal stress (SH = 62.5 MPa). This indicates compressive behavior along the x-axis. In the area far from the wellbore (blue area), Sxx is slightly less than SH. This indicates tensile behavior along the x-axis in the fracture tip area. The stress shadowing along the x-axis induced by the opening fracture is ∆S xx is symmetrically distributed at the x-y plane; the ∆S xx on 1/4 of the x-y plane is given in Figure 11. ∆S xx is also symmetrically distributed at a plane z = 10, so the ∆S xx on the z = −2.5 m, 2.5 m, 7.5 m planes are given in Figure 13a-c, respectively. ∆S xx is greater than zero in the area adjacent to the wellbore, while it is less than zero in the far-field area.
Therefore, an open hydraulic fracture creates stress shadowing on both the maximum and minimum horizontal stress directions. The horizontal stress contrast coefficient after hydraulic fracturing is defined as: The horizontal stress contrast coefficient is less than 1 in the middle area (blue area in Figure 14). It indicates that the horizontal stress contrast in the middle area declines after the hydraulic fracturing process. The horizontal stress contrast coefficient decreases when the area is close to the wellbore. The horizontal stress contrast coefficient is greater than 1 in the area far from the wellbore (red area in Figure 14). This indicates that the horizontal stress contrast in the far-field increases after the hydraulic fracturing process. Previous studies have shown that the probability of creating a complex fracture network increases with decreasing horizontal stress contrast. Hence, the area close to the wellbore is more likely to include a complex fracture network after the intensive-stage fracturing process.  The horizontal stress contrast coefficient after hydraulic fracturing is defined as: The horizontal stress contrast coefficient is less than 1 in the middle area (blue area in Figure 14). It indicates that the horizontal stress contrast in the middle area declines after the hydraulic fracturing process. The horizontal stress contrast coefficient decreases when the area is close to the wellbore. The horizontal stress contrast coefficient is greater than 1 in the area far from the wellbore (red area in Figure 14). This indicates that the horizontal stress contrast in the far-field increases after the hydraulic fracturing process. Previous studies have shown that the probability of creating a complex fracture network increases with decreasing horizontal stress contrast. Hence, the area close to the wellbore is more likely to include a complex fracture network after the intensive-stage fracturing process. The horizontal stress contrast coefficient after hydraulic fracturing is defined as: The horizontal stress contrast coefficient is less than 1 in the middle area (blue area in Figure 14). It indicates that the horizontal stress contrast in the middle area declines after the hydraulic fracturing process. The horizontal stress contrast coefficient decreases when the area is close to the wellbore. The horizontal stress contrast coefficient is greate than 1 in the area far from the wellbore (red area in Figure 14). This indicates that the horizontal stress contrast in the far-field increases after the hydraulic fracturing process Previous studies have shown that the probability of creating a complex fracture network increases with decreasing horizontal stress contrast. Hence, the area close to the wellbore is more likely to include a complex fracture network after the intensive-stage fracturing process.

Perforation Number Optimization
It is assumed that a horizontal well section is perforated in 3 clusters with perforation

Perforation Number Optimization
It is assumed that a horizontal well section is perforated in 3 clusters with perforation cluster spacing of 10 m. The injection volume is 50 m 3 . We only change the perforation number of each perforation cluster into Table 3. The rest parameters are listed in Table 2. Fracture width and surface area are given in Table 3. According to the perforation number distribution, the data in Table 3 are divided into 4 groups, and each group includes 3 simulation results: group 1 (no. 1, 2 and 3); group 2 (no. 1, 4 and 5); group 3 (no. 6, 2 and 7); group 4 (no. 8, 9 and 3).
Group 1: The perforation numbers in this group are equally distributed. The exterior fracture surface area and its centric width increase with the increase of perforation number. The interior fracture surface area and centric width have the tendency to decrease with the increase of perforation number.
Group 2: We retain the exterior perforation number and increase the interior perforation number. The centric width of interior fracture decreases, and the centric width of exterior fracture increases with the increase of interior perforation number. The interior fracture surface area increases, and the exterior one decreases with the increase of interior perforation number. Increasing the perforation number of the interior perforation cluster, the equilibrium of fracture geometry is better. The higher the ratio of interior perforation number to exterior perforation number, the better the fracture geometry equilibrium. A similar rule can be observed in groups 3 and 4.
To improve the fracture geometry equilibrium, it is suggested to design enough perforation numbers in each perforation cluster and make sure that the perforation number in the interior perforation cluster is greater than that of the exterior ones.

Perforation Cluster Number Optimization
It is assumed that the perforation cluster spacing is 10 m. The total injection volume of each simulation case is 50 m 3 . We only change the perforation cluster number, and the rest parameters are listed in Table 2. Fracture surface areas under different perforation cluster numbers and different injection volumes are listed in Table 4. The area of the exterior fracture is still greater than that of the interior one. This phenomenon was proven by the experiments [4,30].
As shown in Table 4, the total fracture surface area in each fracturing stage is given in Figure 15. The total fracture surface area reaches the peak value when the perforation cluster number is 5. The total fracture surface area decreases linearly when the perforation cluster number is greater than 5. The total fracture surface area increases with the increase of injection volume of fracturing fluid. Only increasing the perforation cluster number cannot improve the fracture geometry equilibrium. It is suggested to optimize the perforation cluster number and perforation number in each cluster, aiming at improving the injection rate of the interior fractures. Pushing the fracturing fluid into the interior perforation hole by dropping a degradable ball to seal the exterior perforation hole may be a good choice for improving the fracture geometry equilibrium. As shown in Table 4, the total fracture surface area in each fracturing stage is given in Figure 15. The total fracture surface area reaches the peak value when the perforation cluster number is 5. The total fracture surface area decreases linearly when the perforation cluster number is greater than 5. The total fracture surface area increases with the increase of injection volume of fracturing fluid. Only increasing the perforation cluster number cannot improve the fracture geometry equilibrium. It is suggested to optimize the perforation cluster number and perforation number in each cluster, aiming at improving the injection rate of the interior fractures. Pushing the fracturing fluid into the interior perforation hole by dropping a degradable ball to seal the exterior perforation hole may be a good choice for improving the fracture geometry equilibrium. As shown in Table 4, the centric width of exterior fracture in each stage is given in Figure 16. The centric width of exterior fracture reaches the minimum value when the perforation cluster number is 6. The more the injection volume, the wider the centric width of the exterior fracture. A wider width contributes to the proppant migration, and the perforation cost of the 5-cluster is cheaper than that of the 6-cluster. Therefore, the best choice in this study site is to use the 5-cluster perforation for improving the stimulation efficiency. As shown in Table 4, the centric width of exterior fracture in each stage is given in Figure 16. The centric width of exterior fracture reaches the minimum value when the perforation cluster number is 6. The more the injection volume, the wider the centric width of the exterior fracture. A wider width contributes to the proppant migration, and the perforation cost of the 5-cluster is cheaper than that of the 6-cluster. Therefore,

Perforation Hole Diameter Optimization
It is assumed that a horizontal well section is perforated in 3 clusters with a spacing of 10 m. The injection volume is 100 m 3 . The perforation number of each cluster is 10. The rest parameters are listed in Table 2. We only use the perforation diameter in Table 5 to research the relationship among the fracture width, fracture area and the perforation diameter distribution. According to the perforation diameter, the data in Table 5 are divided into 4 groups, and each group includes 3 simulation results: group 1 (no. 1, 2 and 3); group 2 (no. 1, 4 and 5); group 3 (no. 6, 2 and 7); group 4 (no. 8, 9 and 3). Group 1: The perforation diameters in this group are equally distributed. The exterior fracture area increases with the increase of perforation diameter. The centric width of the exterior fracture decreases, and the centric width of the interior fracture increases with perforation diameter. The surface area of exterior fracture increases with the increase of perforation diameter. Under equal diameter of perforation, the fracture geometry equilibrium is better when the perforation diameter is bigger.
Group 2: We retain the exterior perforation diameter and increase the interior perforation diameter. The centric width of interior fracture decreases, and the centric width of exterior fracture increases with the increase of interior perforation number. The surface area of exterior fracture increases with the increase of interior perforation diameter. Increasing the perforation diameter of the interior perforation cluster, the equilibrium of fracture geometry is better. The higher the ratio of interior perforation diameter to exterior

Perforation Hole Diameter Optimization
It is assumed that a horizontal well section is perforated in 3 clusters with a spacing of 10 m. The injection volume is 100 m 3 . The perforation number of each cluster is 10. The rest parameters are listed in Table 2. We only use the perforation diameter in Table 5 to research the relationship among the fracture width, fracture area and the perforation diameter distribution. According to the perforation diameter, the data in Table 5 are divided into 4 groups, and each group includes 3 simulation results: group 1 (no. 1, 2 and 3); group 2 (no. 1, 4 and 5); group 3 (no. 6, 2 and 7); group 4 (no. 8, 9 and 3). Group 1: The perforation diameters in this group are equally distributed. The exterior fracture area increases with the increase of perforation diameter. The centric width of the exterior fracture decreases, and the centric width of the interior fracture increases with perforation diameter. The surface area of exterior fracture increases with the increase of perforation diameter. Under equal diameter of perforation, the fracture geometry equilibrium is better when the perforation diameter is bigger.
Group 2: We retain the exterior perforation diameter and increase the interior perforation diameter. The centric width of interior fracture decreases, and the centric width of exterior fracture increases with the increase of interior perforation number. The surface area of exterior fracture increases with the increase of interior perforation diameter. Increasing the perforation diameter of the interior perforation cluster, the equilibrium of fracture geometry is better. The higher the ratio of interior perforation diameter to exterior perforation diameter, the better the fracture geometry equilibrium. A similar rule can be observed in groups 3 and 4.

Conclusions
This paper proposes a coupled fracture model for the intensive-stage fracturing technology based on the 3D-DDM. It can simulate the fracture geometry and stress field in the intensive-stage fracturing of horizontal wells. Some interesting findings are: (1) Under the uniform perforation parameters, the exterior fractures become the main fractures with the largest widths, while the interior fractures are strongly restrained and have the smallest widths. The fracture geometry equilibrium is very low under the uniform perforation parameters. Due to the symmetry, the exterior fractures' widths have the same variation versus the injection time. With increasing injection time, the fracture width increases drastically in the beginning and then gradually increases. (2) The inlet's fluid pressure is the highest, while the fluid pressure at the fracture tip is the lowest. The injection pressure declines rapidly in the beginning and then gradually declines later. The final injection pressure tends to the original minimum horizontal stress. (3) The S zz in the middle area is much greater than the original minimum horizontal stress, indicating strong compression along the z-axis. In the area far from the horizontal well, S zz is less than the original minimum horizontal stress, indicating tension along the z-axis induced by the opening fracture. Due to the stress interaction, the interior fractures suffer strong compressive stress, and their propagations are strongly restrained. (4) The horizontal stress contrast in the middle area declines after the hydraulic fracturing process. The horizontal stress contrast in the far-field increases after the hydraulic fracturing process. The area close to the wellbore is more likely to include a complex fracture network after the intensive-stage fracturing process. (5) Only increasing the perforation cluster number in each stage cannot improve the fracture geometry equilibrium in the intensive-stage fracturing. In this study, site 5-cluster perforation is the best choice for improving the stimulation efficiency. To improve the fracture geometry equilibrium, it is suggested to design more perforation numbers in each perforation cluster and ensure that both the perforation number and diameter in the interior perforation cluster are greater than those of the exterior ones.