A Modeling Study of the Productivity of Horizontal Wells in Hydrocarbon-Bearing Reservoirs: Effects of Fracturing Interference

State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum, Beijing 102249, China Beijing Gas Group Company Limited, Beijing 100035, China Geological Research Institute of CNPC Logging Company Limited, Xi’an, Shaanxi 710000, China Supervision Center of CNPC Tarim Oilfield Company, Korla, Xinjiang 841000, China Department of Petroleum Engineering, China University of Petroleum-Beijing at Karamay, Karamay, Xinjiang 834000, China


Introduction
Wells in low permeability reservoirs bearing hydrocarbons typically have low productivity as it is hard for hydrocarbons to efficiently flow. Therefore, horizontal wells and hydraulic fractures are often used to enhance the contact between the wellbore and the matrix [1][2][3][4][5][6]. Since multistage and multicluster hydraulic fracturing facilitates the establishment of complex fracture networks, this type of fracturing technique is widely used in the development of unconventional reservoirs such as shale oil reservoirs [7][8][9]. However, due to the stress changes induced by hydraulic fracture initiation and propagation, multistage and multicluster hydraulic fracturing is affected by stress interference, and the geometry of the fracture network can be negatively impacted. It can lead to nonuniform hydraulic fracture growth and unevenly distributed fracture networks [9]. As fracture quality is directly related to horizontal well productivity, it is meaningful to quantitatively understand the relationship between hydraulic fracturing stress interference and horizontal well productivity. In order to quantify the relationship between stress interference in hydraulic fracturing and the production performance of the fractured horizontal well, it is important to describe the fracture mechanics during hydraulic fracturing and the fluid flow in porous media during horizontal well depletion.
To build an efficient and reliable modeling method for the temporal and spatial evolution of stress, hydraulic fracture width and path, pressure, and hydrocarbon production, adequate modeling techniques should be employed. In geologic media, multiphase fluids flow in the porous media. Hydrocarbon-bearing reservoirs are usually characterized by rock heterogeneity, where fractured and unfractured media play an important role in the evaluation of pressure and stress evolutions [10,11]. In such problems, the coefficient matrices in multiple time steps in the numerical system are usually huge, and the use of compositional and multiphase flow models and multiporosity assumptions further increases the complexity and leads to higher computational costs [12][13][14]. Therefore, the accurate modeling of the stress interference and the multistage fractured horizontal well productivity should be specifically investigated.
Hydraulic fracturing involves the use of high-pressure fluid injection into the reservoir to break the rock, forming fractures with high conductivities. This makes it easier for reservoir fluids to flow into the wellbore, thereby increasing oil and gas production. This process can activate natural fractures and increase the complexity of fracture networks. To understand and optimize hydraulic fracturing parameters, researchers have proposed many mathematical models for the simulation of hydraulic fractures. The early researches generally assumed simplified fracture geometries using the 2D plane strain assumption. The widely used Khristinaovic-Geertsma-de Klerk (KGD) model and Perkins-Kern-Nordgren (PKN) model calculate fracture geometries in two-dimensional planes, and fracture width distributions are calculated as a ratio of length to height [15][16][17][18]. 2D and 3D models were then derived, and the effect of in-situ stress and fluid flow was incorporated [19][20][21].
Industrialized exploitation of shale gas and shale oil reservoirs requires more sophisticated models describing the process of multistage and multicluster hydraulic fracturing. In this process, it is necessary to take into account the effect of the evolution of stress fields as there are interstage and intercluster interferences. As a result, hydraulic fractures with nonuniform half-lengths can be generated. Nonplanar and asymmetric fracture geometries can be generated under the influence of stress interference as well [22,23]. The interference between simultaneously growing fractures can be intensified by decrease fracture spacing, and a minimum spacing should be determined for good fracturing quality [24]. Optimized hydraulic fracture-related parameters including fracturing timing and cluster location were determined in a reservoir-geomechanics-fracturing workflow where the effect of production-induced stress state changes is also considered [25,26]. In addition to numerical modeling, triaxial tests are usually used in the lab to physically understand the fracturing initiation and propagation process [27,28]. Monitoring techniques are also developed to better understand the shape and geometry of fractures in shale reservoirs [29,30].
After the establishment of hydraulic fractures in the horizontal well, numerical simulation techniques are used to compute the fluid flow from the low permeability reservoir to the wellbore [31]. Stress sensitivity was sometimes incorporated in reservoir simulators, and the production prediction was affected by poroelasticity. Typically, the consideration of geomechanical effects tends to decrease the predicted production [32,33]. In another reservoir simulation model, Moradi et al. [34] pointed out that changes in fracture aperture significantly alter the simulated production rates.
Previous studies usually focused on the stress interference phenomenon during the propagation of multiple hydraulic fractures, and the quantitative impact of this phenomenon on the horizontal well productivity has not been thoroughly investigated. This study employs hydraulic fracturing modeling and reservoir simulation techniques and proposes a modeling workflow quantifying the relationship between stress interference and horizontal well productivity. Parametric studies are also conducted to investigate the stage and cluster parameters on horizontal well productivity. This workflow provides a reference for the optimization of hydraulic fracturing parameters based on fracturing quality and horizontal well productivity modeling.

Methodology
In the methodology, a combined modeling workflow consisting of hydraulic fracturing modeling and porous media flow is presented. Finite element methods are used to establish the numerical models. Fracture propagation, geomechanical response, and fluid flows are all considered.
The momentum balance in the stress tensor is used to describe the rock deformation.
The boundary condition for the rock deformation problem has three types. They are the stress boundary, the pressure boundary, and the displacement boundary: where n represents the unit normal vector, Γ t represents the traction boundary, Γ + p ∩ Γ − p denotes the pressure boundary, and Γ u represents the displacement boundary.
The Poiseuille's law is employed to compute the incompressible fracturing fluid flow in fractures as follows: In Equation (5), w is the width of fracture at location l; p is the pressure; μ is the viscosity; q t is the flow rate [35]. 2 Geofluids Then, a lubrication equation can be used to characterize the mass balance in the fracture flow as follows: where w is the fracture width, q f i is the flow rate, and q ts is the fluid loss in the fracture flow.
The fluid leak-off into the host rock can be calculated by the following: where c is the leak-off coefficient, p i is the pressure in the fracture, and p t is the pressure in the formation.
In the fracture mechanism, the normal and shear stresses are depicted by the traction-separation method as follows: where t n , t s , and T t represent the normal and shear components for stresses; D represents the damage. Once the fracture is initiated, the damage factor gradually increases and the fracture propagation can be described [36]. A two-phase black oil model is used to calculate the production in the hydraulically fractured horizontal well. The mass balance equations are as follows: where s o and s w are saturation values for oil and water, ρ o and ρ w are densities for oil and water, ϕ is porosity, q is the sink/source term, and t is time. The terms of ∂ðρ o ϕs o Þ/ ∂t and ∂ðρ w ϕs w Þ/∂t describe the accumulation of oil and water flows in porous media. ∇⋅ ðρ o v o Þ and ∇⋅ ðρ w v w Þ terms represent the fluxes. Darcy's law is widely used for fluid flows in porous media with low flow rates. It is used in this model as follows: where k r is the relative permeability, μ is the viscosity, v is the velocity, k is the permeability, g is the gravitational acceleration, and H is the depth. In this model, it is assumed that the hydraulic fracture network is fully propped and a fracture permeability is prescribed to denote the fracture conductivity.
Putting Darcy's law and the mass balance together, where more detailed forms of fluid flow diffusivity are obtained.
In the two-phase black oil model, the relationship between water and oil saturations and the initial sink/source rates can be written as follows: Based on the assumption of slightly compressible fluids in the reservoir, the compressibility water and oil can be defined as C w = 1/ρ w ∂ρ w /∂p w and As a result, the flow diffusivities in Equations (12) and (13) can be extended as follows: Neglecting capillary pressure terms, Equations (16) and (17) become the following: For numerical treatment, a matrix form can be obtained as follows: An initial condition is as follows: The boundary condition is as follows: where v is a velocity tensor. Equation (21) represents a noflow boundary condition for the flow problem.

Modeling Study
A synthetic scenario in the development of a shale oil reservoir is established for the numerical study. Modeling parameters in the synthetic case are based on a realistic field in Junggar Basin, northwestern China [37]. Table 1 records the modeling parameters related to hydraulic fracturing while Table 2 shows the parameters for the modeling of fractured horizontal well production in reservoir simulation.
In the modeling workflow, the fracturing and production from a segment in the horizontal wellbore are considered. In the base case, totally four stages are hydraulically fractured in the horizontal wellbore. In each stage, there are three fracture clusters. A fracturing fluid injection rate of 12 m 3/ min per stage is used. The base case fracture stage spacing is 60 m. In the parametric study for hydraulic fracturing, the effects of stage number, cluster number, stage spacing, and injection volume are quantified. Then, the horizontal well production and pore pressure depletion in the reservoir corresponding to each hydraulic fracturing scenario are simulated. Thus, the relationship between the stress interference in hydraulic fracturing and the productivity can be quantified using this workflow.
In addition to the hydraulic fracturing modeling and productivity modeling in the aforementioned base case, several parametric studies are also carried out. In the parametric study for the effect of fracture stage and cluster, two scenarios are considered. One scenario has two stages and six clusters in each stage, and the other scenario has three stages with four clusters in each stage. In the parametric study for the effect of fracture stage spacing, another two scenarios of 40 m and 60 m stage spacings are considered. In the parametric study for the effect of fracturing fluid injection rate per stage, two scenarios of 10 m 3 /min and 14 m 3 /min are modeled.
3.1. Hydraulic Fracturing Modeling. In this section, the effects of fracture stage/cluster number, stage spacing, and fracturing fluid injection rate on the propagation of the hydraulic fracture networks are modeled. The resultant fracture length, fracture width, fracture area, and fracture   Figure 1 shows the nonuniform hydraulic fractures in the horizontal well with the three different stage and cluster numbers. In general, the four-cluster design leads to the most nonuniform hydraulic fracture growth in the wellbore, with the two outer fracture clusters longer than the inner fracture clusters. Note that the stage spacing is kept as 60 m in all the scenarios. To better quantify the effect of stage and cluster design on the nonuniform fracture lengths, Figure 2 is plotted where the length of each individual fracture cluster is presented as bars. In Figure 2, the labels A to D represent the four fracture stages, and 1 to 3 represent the clusters in a stage. Note that A represents the first fractured stage while D is the last stage that is fractured in the multistage fracturing job. Similarly, labels E to J represent the six fracture stages in another design, where 1 to 2 are the two clusters in each stage. K to M are the three stages in the 4-cluster scenario. Based on the detailed fracture length results, it is noted that when the cluster number exceeds two in each stage, the inner fractures are always shorter than the outer fractures. This indicates that the stress interference in simultaneously growing clusters inhibits the growth of the inner fractures and makes the outer fracture more competitive in terms of fracture propagation [38]. The 2-cluster scenario results in a more evenly distributed fracture length. This indicates that reducing the clusters in each stage can lower the stress interference effect on the heterogeneous growth of fracture clusters. In this scenario, only two fractures grow simultaneously each time. Therefore, the interference between clusters is reduced, leading to a more uniform growth of fractures. In the 4-cluster scenario, since four fracture clusters are competing in the simultaneous growth, the inhibition on the inner fractures becomes more noticeable. Also, the average fracture lengths are the shortest in the 4cluster scenario. Generally, the longest fracture length of an individual fracture is obtained in the first stage in the 3-cluster scenario. This indicates that the first stage is less affected by interstage interference than the stages fractured later on. In addition, in the 2-cluster, 3-cluster, and 4-cluster scenarios, it is noted that the first stage has very symmetric fractures, indicating that the interference does not take effect unless there are sequentially fractured stages. The interstage stress interference is the most significant in the 4-cluster scenario, as the second and third stages have lower fracture lengths. These results show that the interstage interference on hydraulic fracturing interference increases with the cluster number in each stage. However, it is noted that although reducing the number of clusters in each stage helps to form uniform fracture networks, it may increase the cost of the fracturing operations and fracturing time. Figure 3 presents the temporal evolution of fractured volume and fracture area during the multistage fracturing operation with different fracture stage and cluster designs. This helps to improve the understanding of the correlation between stage and cluster designs and fracture quality.
Step-by-step trends are obvious in these results, as stages are sequentially fractured in the operation. The final fractured volumes for the 3-cluster and 2-cluster scenarios are very similar, while the selection of 4-cluster design leads to the lowest fractured volume. This is direct evidence that increased fracturing interference reduces the final fracture  5 Geofluids network quality. Based on the fractured area results, the 2cluster scenario has the best fracture quality and the 3cluster scenario has the intermediate quality. Again, the 4-cluster scenario leads to the lowest fracture quality. Note that the change from 3-cluster to 2-cluster design increases the contact area between fractures and the low permeability matrix, implying an improved fracturing performance.

Fracture Stage Spacing.
In this section, the effect of the fracture stage on the hydraulic fracture quality is investigated. Three stage spacings of 40 m, 60 m, and 80 m are studied while other fracturing-related parameters are the same as the base case as in Table 1. Figure 4 describes the fracture geometries from three different stage spacings. Note that the cluster spacing within each stage is kept the same during the parametric study. In general, when the stage spacing is reduced to 40 m, a more nonuniform fracture length pattern is observed. In contrast, when the stage spacing is increased to 80 m, a rather uniform distribution of fracture length is obtained. However, it is noted that the inner fracture growth is always inhibited by fracturing interference even when the spacing is large. It means that the inhibition on the inner fracture is caused by intercluster interference instead of interstage interference. Figure 5 records the comparison of fracture length of each cluster in scenarios with different stage spacings. Since all three scenarios have the same stage number of four, labels A to D are used to represent the four stages and each stage has three clusters. Based on the fracture length results, the effect of stage spacing is not quite monotonic. The correlation between stage spacing and fracture cluster length is not clear. Figure 6 shows the temporal changes in the fractured volume and the fractured area during the 4-stage fracturing. The fracturing of each individual stage leads to a sharp increase in the fractured volume and area. The fractured volume results indicate that the 80 m spacing leads to the greatest volume while the 60 m spacing has the lowest fractured volume. The fractured area results show that the 40 m spacing leads to the highest fractured area.

Geofluids
The geometries are shown in Figure 7. Compared to the fracturing fluid injection rates 10 m 3 /min and 12 m 3 /min, the fracturing fluid injection rate of 14 m 3 /min leads to a more nonuniform fracture geometry, and the nonuniform lengths are more significant between the first and the second stages. This shows that an increase in injection rate leads to a more unevenly distributed fracture network, and an increased injection rate corresponds to an elevated fracturing interference between hydraulic fracturing stages. Figure 8 shows the cluster-by-cluster comparison of fracture length with various fluid injection rates. The correlation between injection rate and fracture length is well in the first stage (stage A). In the first stage, the longest fracture lengths are obtained when the injection rate is 14 m 3 /min, and the fracture lengths decrease with the decrease in injection rate.

Geofluids
The general trend is the same for stages B to D, while oscillations in fracture length are observed. This is caused by the nonuniform results obtained in the simulation for the fracture mechanism.
In Figure 9, the stepwise increases in the fractured volume and area are plotted against fracturing time. In these results, the correlation between the injection rate and the fractured volume/area is clear. With the increase in fracturing fluid injection, the step-by-step fractured volume and area and the final fractured volume and area both increase. Since the total injection time is constant, a higher fracturing fluid injection rate corresponds to a greater volume of fluid injected into the fractures for fracture initiation and propagation.
Based on these parametric studies for the effect of fracture stage/cluster design, stage spacing, and fracturing fluid injection rate, stage spacing is less influential than stage/cluster design and injection rate. The reduction of cluster number per stage and the increase in stage spacing help to establish more evenly propagated fractures. The increase in injection rate can elevate the interstage fracturing interference by making the fracture lengths more uneven, while a greater injection rate improves the overall fractured volume and area.
The study in this section quantifies the hydraulic fracture geometries in the horizontal well, and the effect of fracturing interference is investigated in terms of fracture geometry. However, the investigation of horizontal well productivity requires results more than fracture geometry, and the productivity-related parameters such as pore pressure and hydrocarbon production should be quantified.

Productivity Modeling.
In the previous section, fracture geometry, fractured volume and area, and fracture length are used as the variates to denote the effect of fracturing interference. To further investigate the horizontal well productivity, reservoir simulation techniques are used to calculate the production of the horizontal well with different hydraulic fracture geometries obtained in the previous fracturing modeling. Thus, the relationship between fracturing interference and horizontal well productivity can be established. In the modeling process, the fracture conductivity is assumed to be constant in each fracture and it does not change with time. The fracture conductivity is calculated as the product of the fracture width and a constant permeability value for the fracture. Modeling parameters are shown in Table 2. Productivity over 2 years of production is reported. Figures 10 and 11 compare the cumulative production of oil over 2 years from horizontal wells with three different cluster and stage designs, stage spacings, and fluid injection rates.
In Figure 10, the use of six stages with two clusters in each stage leads to the highest cumulative production. This is because this strategy has the lowest fracturing interference, and it can lead to even depletion within the low permeability reservoir. The use of four stages with three clusters in each stage (the base case) leads to intermediate cumulative production performance. In contrast, the use of four clusters in each stage results in a much lower cumulative oil production curve. Combined with results in Figures 1 and 2, it is noted that the use of four clusters in a stage largely inhibits the growth of fractures, especially for the two inner fractures in each stage. In this scenario, the average fracture length is the lowest, leading to lower production performance.
In Figure 12, how stage spacing affects cumulative production from a horizontal well is presented. In general, the  8 Geofluids effect of stage spacing is not as great as stage and cluster design, and the differences in cumulative production results are smaller. The base case of a 60 m spacing yields the highest cumulative production, while the stage spacing of 40 m has the lowest cumulative production. Based on the observations in Figures 4 and 5, reducing the stage spacing to 40 m strengthens the fracturing interference and leads to uneven fracture lengths, where the inhibition on the growth of certain fractures is also increased. Cumulative oil productions from horizontal wells with three fracturing fluid injection rates during the hydraulic fracturing process are compared in Figure 11. Intuitively, the highest injection rate corresponds to the greatest cumulative production, and the lowest injection rate leads to the smallest cumulative production, which is also correlated with the fracture geometries and lengths in Figures 7 and 8.
Based on productivity modeling, the fracture geometry and the horizontal well productivity can be generally correlated: horizontal wells with longer fracture lengths and weaker fracturing interference usually yield higher cumulative production. However, the differences in horizontal well productivity from multiple hydraulic fracturing design   Pressure (Pa) Figure 13: Pressure distribution after three months of production for multistage scenarios of three stages, four stages, and six stages. 10 Geofluids scenarios cannot be directly quantified by hydraulic fracturing modeling, and reservoir simulation has to be used to obtain the detailed differences in cumulative oil production.

Discussion
Note that in this modeling study, planar hydraulic fractures are considered while it is also possible to have nonplanar hydraulic fractures generated in multistage hydraulically fractured horizontal wells. Nonplanar fracture modeling can better quantify the effect of stress interference on resultant hydraulic fracture geometries. For example, curved fractures with uneven lengths can be obtained using this modeling technique. However, the focus of this study is on the correlation between stress interference and horizontal well productivity. As shown in Figure 13, hydraulic fracture geometries directly govern the drainage area. After three months of depletion, the pressure drop fronts are around 10 m away from the fractures, indicating that the hydrocarbons within this area are being produced. Therefore, the fracture length and the drainage area jointly govern the resultant horizontal well productivity. In this workflow, the effects of curvatures of nonplanar hydraulic fractures on horizontal well productivity are weakened. In consequence, although planar fracture modeling used in this study cannot characterize how fractures are curved by stress interference, the obtained fractures are unevenly distributed and the effects on fracture lengths and the resultant drainage area are still honored.

Conclusion
In this modeling study, based on a numerical modeling workflow, the relationship between the fracturing interference during multistage hydraulic fracturing and horizontal well productivity is established, which is the primary contribution in the study. Using hydraulic fracture geometries to quantify the effect of fracturing interference on productivity is not comprehensive. Based on the modeled fracture geometries including fracture length, fractured volume, and fractured area, an estimate can be obtained for productivity. However, quantitative understanding should be obtained using a more comprehensive workflow including fracturing modeling and productivity modeling. In conclusion, (1) Reducing the cluster number within a stage, increasing the spacing between two stages, and reducing the fracturing fluid injection rate help to decrease the negative impact of fracturing interference on multistage and multicluster fracturing (2) Fractured volume and fractured area are parameters indicating the magnitude of the stimulated reservoir volume/area. They are generally correlated with fracture lengths. The temporal evolution patterns of fractured volume and area are highly stepwise, which corresponds to the sequentially fractured stages in the field operation (3) Productivity modeling results are generally correlated with fracturing modeling results. However, productivity modeling is capable of providing the quantitative differences in cumulative production, which directly shows the effect of fracturing interference on horizontal well production performance

Data Availability
The data are available by contacting the corresponding author.

Conflicts of Interest
The authors declare that they have no conflicts of interest