Abstract
An appropriate well spacing plan is critical for the economic development of shale gas reservoirs. The biggest challenge for well spacing optimization is interpreting the subsurface uncertainties associated with hydraulic and natural fractures. Another challenge is the existence of complex natural fractures. This work applied an integrated well spacing optimization workflow in shale gas reservoirs of the Sichuan Basin in China with both hydraulic and natural fractures. The workflow consists of five components: data preparation, reservoir simulation, estimated ultimate recovery (EUR) analysis, economic calculation, and well spacing optimization. Firstly, the multiple realizations of thirteen uncertain parameters of matrix and fractures, including matrix permeability and porosity, three relative permeability parameters, hydraulic fracture height, half-length, width, conductivity, water saturation, and natural fracture number, length, and conductivity, were captured by the assisted history matching (AHM). The fractures were modeled by the embedded discrete fracture model (EDFM) accurately and efficiently. Then, 84 AHM solutions combining with five well spacing scenarios from 517 ft to 1550 ft would generate 420 simulation cases. After reservoir simulation of these 420 cases, we forecasted the long-term gas production for each well spacing scenario. Gas EUR degradation and well interference would imply the critical well spacing. The net present value (NPV) for all scenarios would be calculated and trained by -nearest neighbors (KNN) proxy to better understand the relationship between well spacing and NPV. In this study, the optimum well spacing was determined as 793 ft, corresponding with a maximum NPV of 18 million USD, with the contribution of hydraulic fractures and natural fractures.
1. Introduction
There is no doubt that the development of unconventional reservoirs has changed the oil and gas industry. However, many challenges, such as heterogeneity, nanopore, proppant distribution, multiphase flow, and complex fractures, have existed in unconventional reservoirs [1–3]. Reservoir simulation is a rigorous method applied in unconventional reservoirs. Among the worldwide unconventional resources, shale gas and oil are the main components. Optimum well spacing is one of the key parameters for shale reservoir development. It is essential to find a well spacing that can balance the recovery and economics. Many studies have focused on this area both numerically and analytically [4, 5]. Some of them have investigated the controlling factors for the well spacing determination, such as fracture half-length, reservoir permeability, rock properties, and natural fractures [6–8]. However, due to shale reservoirs’ complexity, it is still challenging to quantify the subsurface uncertainty-associated hydraulic and natural fractures [9, 10]. Several methods, including the microseismic method, well test method, and rate-transient analysis (RTA) method, have applied to capture the uncertain parameters [11–14]. Nevertheless, the high data requirement of these methods makes it not easy to be performed in new wells.
Another method widely used in calibrating uncertainties is history matching, which is an efficient and inexpensive approach. Many authors have performed single history matching to obtain one solution of shale reservoir and fractures [15–19]. However, the nonuniqueness of history matching should be considered. Therefore, multiple history matching was applied. Cao et al. [20] determined the optimal well spacing for Delaware Basin by multiple history matching. What is more, they did not take the uncertainty of natural fractures into account, which is another challenge for the well spacing optimization in shale reservoirs.
Several researchers have investigated that natural fractures could impact the fracture’s propagation during hydraulic fracturing by microseismic event patterns and complex fracture propagation models [21–24]. [2] modeled the complex natural fractures by the embedded discrete fracture model (EDFM), a modeling method with accuracy and efficiency. They indicated that the two-set natural fractures could increase the gas recovery by 23.2% after 30 years.
In this study, we applied an integrated AHM and EDFM workflow for well spacing optimization in shale gas reservoirs of Sichuan Basin in China with complex natural fractures. The hydraulic fractures and natural fractures were modeled by the EDFM method [2]. According to 84 AHM solutions for a shale gas well in this reservoir, the multiple realizations of thirteen uncertain matrix and fracture parameters can be calibrated. It is worth pointing out that the uncertain parameters of natural fractures include the number, the length, and the conductivity of natural fractures. We also considered three uncertainties about relative permeability. Then, we compared the gas EUR in the long-term from the reservoir simulation results for five well spacing scenarios associated with these 84 solutions. The well spacing scenarios are distributed from 517 ft to 1550 ft, corresponding to 6 wells to 2 wells. The critical well spacing can be determined by analyzing the gas EUR degradation to minimize well interference. Then net present values (NPVs) of all cases can be evaluated and predicted by -nearest neighbors (KNN) proxy to identify the optimum well spacing for this shale reservoir with natural fractures.
2. Well Spacing Optimization Workflow
The integrated well spacing optimization workflow consists of five components: data preparation, reservoir simulation, EUR analysis, economic calculation, and well spacing optimization. The framework is shown in Figure 1. First, we need to prepare the input data for the reservoir simulation. There are two main things we need to decide. One is which uncertain parameters of matrix and fractures are essential to the shale gas reservoir. According to the short-term production data, we apply AHM to calibrate the distribution of these uncertain parameters and screen the results of AHM solutions with the lowest global error. Another is to design the minimal and maximum numbers of well placed into the reservoir to figure out the optimum well spacing. If the well spacing is too small, the well interference will reduce the gas production per well; if the well spacing is too large, the recovery may not satisfy. Therefore, it is essential to determine the range of well spacing and design several well spacing scenarios. The input of simulation cases can then be generated by integrating the AHM solutions and well spacing scenarios and modeled by EDFM considering hydraulic and natural fractures.

Subsequently, the reservoir simulation is performed for all cases to forecast the long-term gas production. Cumulative gas production and gas EUR of each case can be calculated and analyzed in a lognormal probability plot. By comparing gas EUR degradation, the influence of well interference is observed, and the corresponding critical well spacing minimizing well interference can be obtained. Next, we evaluate the NPV for all the cases and plotted them into a boxplot. The equation of NPV is discussed in our previous work [25]. The P50 NPV for each well spacing scenario can be obtained directly. Besides, to better understand the relationship between NPV and well spacing, we predict the NPV for more well spacing using the KNN proxy method. The calculated NPV is considered as a predictor, and the corresponding well spacing is added into the prediction features. Finally, we can identify the optimum well spacing, which leads to the maximum NPV.
3. Field Application
3.1. Reservoir Model
Our integrated AHM and EDFM workflow is applied to a shale gas reservoir in the Sichuan Basin in China with complex natural fractures to determine the optimum well spacing for hydraulic-fractured wells. First, it is essential to build a field-scale model to represent the shale gas reservoir. The schemes of the model are shown in Figure 2. The model is 5840 ft long in the -direction and 3100 ft long in the -direction. The thickness is 65 ft in the -direction. The red lines are the horizontal wells with a constant length of 4921 ft. We set 2 to 6 wells in the model to illustrate the different well spacings. Figure 2(a) represents the well spacing of 1550 ft at two wells per section, while Figure 2(b) shows the well spacing of 517 ft at six wells per section. The blue lines distributed in the wells represent the 54 hydraulic fractures. They were separated into 18 stages, which are 145 ft away from each other. And each stage contains 3 clusters with cluster spacing of 67 ft. The green lines are the natural fractures distributed at 45° or 135°. Although there are other degrees of fracture growth azimuth existing, we would only consider the ideal condition to simplify the method. One of the biggest reasons is lacking fracture diagnostic data. In this study, we assume the natural fractures have a constant height of 65 ft and a width of 0.1 ft. And the dip angle of natural fractures is 90°. All fractures were modeled by EDFM. It is worth noting that the model is ideal, and the heterogeneity of the reservoir is not considered. Other properties of our reservoir model are listed in Table 1.

(a)

(b)
After demonstrating all the certain parameters, we need to calibrate the uncertain parameters of matrix and fractures for the shale gas reservoirs by history matching. The first step is to determine the critical parameters and their range based on prior expert experience. This study chose 13 parameters as uncertainties: matrix permeability, porosity, exponent of relative permeability for gas, exponent, and endpoint of relative permeability for water, hydraulic fracture height, half-length, conductivity, water saturation, and width, natural fracture number, length, and conductivity. These three parameters of natural fractures reflect the ability of fluid transport within the complex natural fractures. We limited the number of natural fractures from 200 to 1200 in this reservoir. Moreover, the natural fracture length is between 100 ft to 500 ft, while the hydraulic fracture half-length is between 200 ft and 780 ft. Other parameters’ ranges are listed in Table 2. The assisted history matching was then performed to capture the multiple realizations of these parameters. Details of AHM workflow are discussed Tripopoom et al. [26]. For this reservoir, there are a total 84 assisted history solutions. Figures 3(a)–3(d) show the AHM results compared with short-term production data [26]. It is observed that the bottomehole pressure (BHP) and gas flow rate can match with the field production data properly. The water flow rate and water-gas ratio (WGR) match the production data overall, except for 50 days to 250 days, but we can find the water flow rate from production data during this period changes without the similar trend with BHP, which can be ignored for the overall results. Therefore, the history matching solutions are accurate enough for the following well spacing optimization.

(a)

(b)

(c)

(d)
The multiple realizations of uncertain parameters from AHM solutions are shown in Figure 4. Figures 4(a)–4(m) represent the posterior distribution of matrix permeability, fracture height, fracture half-length, fracture water saturation, fracture width, fracture conductivity, matrix porosity, three components of relative permeability, number of natural fractures, nature fracture length, and nature fracture conductivity, respectively. Each plot’s -axis reflects the range of this uncertain parameter, and the range is divided into ten bins. The -axis is the probability of each bin. We combined the AHM solutions with well spacing scenarios, which would generate 420 cases. One dot represents one possible case. It can be easily observed the highest probability of uncertain parameters with most points in a specific bin. The distribution of natural fracture number is shown in Figure 4(k), which implies that more than half of the points distribute between 200 and 500. Moreover, the bin of 400 to 500 has most points compared with other bins. Therefore, the possible value of the number of natural fractures would be 200 to 500, especially in the range of 400 to 500. Similarly, we can find that the possible natural fracture length is about 100 ft to 140 ft, and the possible natural fracture conductivity is 6.8 md-ft to 7.2 md-ft. Compared with hydraulic fracture half-length, concentrating on 316 ft to 374 ft, the natural fractures are much shorter. And the hydraulic fracture conductivity is about 67 md-ft to 105 md-ft, which is larger than that of natural fractures. It indicates that hydraulic fractures would contribute more to gas production than natural fractures. The probability distribution of other uncertainties can be analyzed in the same way.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)
To have an intuitive embodiment of natural fractures’ properties, we built the fracture model using EDFM for different natural fracture numbers and lengths. Figure 5 illustrates the minimal, average, and maximum number of natural fractures: 211, 453, and 1196. All other properties are the same. Moreover, Figure 6 shows minimal, average, and maximum length of natural fractures: 102 ft, 1855 ft, and 456 ft. The more and the longer natural fractures, the better communication through the fracture system.

(a)

(b)

(c)

(a)

(b)

(c)
In addition, we need to consider the relative permeability of different realizations. The relative permeability of water and gas can be obtained using the following equations: where is the water relative permeability, and is the gas relative permeability. is the endpoint of water relative permeability, is the water exponent, and is the gas exponent. where is the normalized water saturation, is water saturation, is irreducible or residual saturation of water, and is residual saturation of gas for a water/gas displacement.
The distribution of endpoint of water relative permeability and the exponent of relative permeability to gas and water is reflected in Figures 5(h)–5(j). The relative permeability curves related to water saturation of 84 solutions are plotted in Figure 7. Red solid lines represent gas relative permeability, and blue lines represent water relative permeability. It is reflected that water relative permeability distributes wider than that of gas, which is corresponding with the greater number of uncertain parameters used in the calculation equations.

After history matching, we captured the pressure distribution in 2 years based on the 690-day production data. The initial pressure is 8000 psi used in the reservoir simulation. It drops to 1000 psi dramatically in 2 years. Then, the pressure remains constant at 1000 psi in the following 18 years. It implies that most gas will be produced in the first five years, especially in the first two years. It is also worth pointing out that the pressure drawdown rate decreases slightly after 100 days compared to that within 100 days.
3.2. Production Analysis
After preparing the reservoir model, we generated the 420 input cases by integrating the 84 history matching solutions with five well spacing scenarios. Then, we performed reservoir simulation for the long-term production simulation of all cases.
Firstly, the gas flow rate and water flow rate change of all cases in 20 years can be obtained, as shown in Figure 8. Different cases of the same well spacing scenario were plotted with the same color. During the production, the gas flow rate decreases first due to the pressure drop. Then, it increases a bit. It is because the pressure decrease becomes slower than the beginning, shown in Figure 9. It would lead to a larger pressure difference between the reservoir and the wellbore. Therefore, more gas and water would be produced. If the pressure drawdown rate did not change, the flow rate would decrease continuously. After two years, the gas flow rate reaches a peak and turns to decrease, which means the bottomhole pressure drops to the target pressure of 1000 psi. The water flow rate has the same trend. This is the result of the pressure change. Moreover, it is observed that the six-well scenario with well spacing of 517 ft has the largest gas flow rate and water flow rate at the beginning. The gas flow rate drops to 0.2 MMscf per day averagely for the six-well scenario, which is the lowest among all scenarios. It implies the existence of well interference for the small well spacing.

(a)

(b)

The cumulative gas production and cumulative water production in 20 years can be obtained and plotted in Figure 10. The cumulative gas production increases dramatically in the first two years. Finally, the cumulative gas production is in the range of 10 billion cubic feet (Bcf) to 20 Bcf. The cumulative water production increases from 200 MSTB to 400 MSTB averagely with the well number increases. In addition, cumulative gas production is more concentrating than cumulative water production. It reflects the uncertainty of water and gas relative permeability.

(a)

(b)
To figure out how well spacing influences gas production, gas EUR per well for each well spacing scenario is calculated and plotted in the lognormal cumulative probability plot, as shown in Figure 11. The -axis is the lognormal cumulative probability, and the -axis shows the gas EUR per well. Five well spacing scenarios are represented in different colors. More dots falling on the corresponding dashed fitting line indicate a higher possibility of the gas EUR distributes in the lognormal format. What is more, the larger the well spacing, the larger the gas EUR per well. We can find the P50 gas EUR for each scenario, listed in Table 3.

Then, the P50 gas EUR degradation related to well spacing is calculated and plotted in Figure 12. It implies that when the well spacing is smaller than 775 ft, the gas EUR degrades significantly up to 37%. This means there is a strong well interference under this well spacing. Then, the well interference becomes smaller and smaller when then well spacing increases from 775 ft to 1033 ft. When the well spacing is larger than 1033 ft, the degradation curve changes is little, which means the well interference could be ignored. Therefore, we need to make sure the well spacing is larger than 775 ft to reduce the influence of well interference.

3.3. Well Spacing Optimization
Next, we calculated the net present value (NPV) to determine the optimal well spacing. The values of NPV calculation inputs are listed in Table 4. The operation cost is 4.5 million USD per well, and the gas price is 1.8 USD/MScf. Other values can be seen in the table.
The result of NPV is plotted in the boxplot, as shown in Figure 13. The -axis represents the NPV while the -axis represents the well spacing. The NPV result of each scenario is drawn in one box. The three lines of the box from higher to lower are the P25, P50, and P75 NPV of each scenario. The highest line outside the box means the maximum NPV, and the lowest line outside the box represents the minimum NPV. There exists a maximum NPV along with the well spacing. Therefore, the optimum well spacing is 775 ft at where the NPV reaches the highest. The highest P50 NPV is about 18 million USD.

Then, we applied the KNN proxy model to train the data and predict the NPV relationship with well spacing. Each point in Figure 14 represents one prediction result and is regressed on one polynomial curve, shown as a blue line. It is observed that the optimal well spacing is 793 ft with 18 million USD. Also, this well spacing can satisfy the requirement of critical well spacing. It is worth pointing out that this result is similar to the boxplot result, but this curve is smoother with more points. Therefore, the KNN proxy method could provide a more reliable result for the optimum well spacing decision.

3.4. Pressure Distribution Visualization
Finally, we combined the optimal well spacing of 793 ft with the lowest global error AHM solution to perform the reservoir simulation and predict the pressure distribution in the long time. The matrix permeability is 202 Nd, and porosity is 0.04. The water exponent, the gas exponent, and the endpoint of relative water permeability are 2.46, 3.91, and 0.55, respectively. Fracture height, fracture half-length, fracture conductivity, fracture water saturation, and fracture width are 44 ft, 289 ft, 194 md-ft, 0.82, and 0.5 ft, respectively. The number of natural fractures is 418, and its length is 270 ft with a conductivity of 8.09 md-ft. Next, we modeled the fractures by EDFM. After the reservoir simulation, the pressure distribution in the matrix is shown in Figure 15. The pressure drops significantly in the first five years, which implies that most gas would be produced in the first 5 years. However, we would like to perform long-term EUR prediction. Therefore, we would predict the 20-year production in this study. Figure 16. shows the pressure distributions in the fracture. We can find that the pressure drop is faster of the hydraulic fractures than natural fractures far away from the well. It drops from 8000 psi to 3000 psi. This implies that hydraulic fractures would play a more critical role in the gas production process in the early time. In addition, the hydraulic fracture pressure drops to low pressure in the first five years. And the pressure of the natural fractures far from the wellbores drops more than that of hydraulic fracture after five years, which implies natural fractures contribute more after five years. The drainage volume is illustrated intuitively in Figure 17, clearly showing the strong well interference after five years of production.

(a)

(b)

(c)

(a)

(b)

(c)

(a)

(b)

(c)
4. Conclusions
This study applied the well spacing optimization workflow for shale gas reservoirs with hydraulic and complex natural fractures in the Sichuan Basin by integrating AHM and EDFM. A total of 84 AHM solutions and five well spacing scenarios were used to predict gas EUR and NPV in 20 years. Then, the optimum well spacing was identified. We summarize the following conclusions from the study: (1)The maximum NPV is around 18 million USD, whether directly calculated or predicted by the KNN proxy model. The corresponding optimum well spacing is 775 ft and 793 ft, respectively, for the two methods. Moreover, the optimum well spacing from KNN proxy is more accurate as considering more well spacing scenarios(2)The gas EUR degradation reaches to 35% when the well spacing is 517 ft, which shows substantial well interference. And when the well spacing is more extensive than 775 ft, the degradation rate starts becoming slow. It indicates that the influence of well interference turns to small(3)The critical well spacing to avoid the influence of well interference is 775 ft. Therefore, optimum well spacing obtained from the two methods is satisfied with this critical spacing(4)The pressure drop of hydraulic fractures is faster than natural fractures, which implies that hydraulic fractures are more important for early time shale gas production. Then, for a longer time, the natural fractures contribute more to the gas production
Acronyms
AHM: | Assisted history matching |
BHP: | Bottomhole pressure |
EUR: | Estimated ultimate recovery |
EDFM: | Embedded discrete fracture model |
HM: | History matching |
KNN: | K-nearest neighbor |
LGR: | Local grid refinement |
NPV: | Net present value. |
: | Fixed well maintenance cost |
: | Total cost of a specific month |
: | Water disposal cost |
: | Total individual well cost |
: | Distance associated with this connection. |
: | Gross income of a specific month |
: | Permeability |
: | Number of wells |
: | Net present value of a specific scenario |
: | Gas price (dollars per million standard cubic feet) |
: | volume flow rate between two cells in a NNC pair |
R: | Annual discount rate |
: | Miscellaneous tax rate |
: | Gas tax rate |
: | Transmissibility factors in each of NNC pair |
: | Gas production in a specific month (million standard cubic feet) |
: | Water production in a specific month (barrels) |
: | Relative mobility. |
ft×3.048e-01 = m: | |
ft3×2.832e-02 = m3: | |
psi×6.895e+00 = kPa: | |
md×1e-15e+00 = m2: |
Bcf: | Billion standard cubic feet |
MMscf: | Million standard cubic feet |
MSTB: | Thousand stock tank barrel |
md: | Millidarcy |
nd: | Nanodarcy. |
Data Availability
Data are available upon request.
Additional Points
Highlights. (1) Optimum well spacing for shale gas reservoirs was obtained. (2) Both the influences of natural fractures and hydraulic fractures were considered. (3) AHM calibrated thirteen uncertainty parameters of fractures and matrix. (4) EDFM was performed to establish fracture models. (5) Maximum NPV predicted from the KNN model determined the optimum well spacing.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors would like to acknowledge Sim Tech LLC for providing EDFM software for this study.