A Novel Optimization Workﬂow Coupling Statistics-Based Methods to Determine Optimal Well Spacing and Economics in Shale Gas Reservoir with Complex Natural Fractures

: In order to account for big uncertainties such as well interferences, hydraulic and natural fractures’ properties and matrix properties in shale gas reservoirs, it is paramount to develop a robust and e ﬃ cient approach for well spacing optimization. In this study, a novel well spacing optimization workﬂow is proposed and applied to a real shale gas reservoir with two-phase ﬂow, incorporating the systematic analysis of uncertainty reservoir and fracture parameters. One hundred combinations of these uncertainties, considering their interactions, were gathered from assisted history matching solutions, which were calibrated by the actual ﬁeld production history from the well in the Sichuan Basin. These combinations were used as direct input to the well spacing optimization workﬂow, and ﬁve “wells per section” spacing scenarios were considered, with spacing ranging from 157 m (517 ft) to 472 m (1550 ft). An embedded discrete fracture model was used to e ﬃ ciently model both hydraulic fractures and complex natural fractures non-intrusively, along with a commercial compositional reservoir simulator. Economic analysis after production simulation was then carried out, by collecting cumulative gas and water production after 20 years. The net present value (NPV) distributions of the di ﬀ erent well spacing scenarios were calculated and presented as box-plots with a NPV ranging from 15 to 35 million dollars. It was found that the well spacing that maximizes the project NPV for this study is 236 m (775 ft), with the project NPV ranging from 15 to 35 million dollars and a 50th percentile (P50) value of 25.9 million dollars. In addition, spacings of 189 m (620 ft) and 315 m (1033 ft) can also produce substantial project proﬁts, but are relatively less satisfactory than the 236 m (775 ft) case when comparing the P25, P50 and P75 values. The results obtained from this study provide key insights into the ﬁeld pilot design of well spacing in shale gas reservoirs with complex natural fractures.


Introduction
In 2013, the U.S. Energy Information Administration (EIA) gathered data of global shale gas and oil reserves data. The technically recoverable shale oil and gas are 47 billion tons (345 billion barrels) and 207 trillion cubic meters (7299 trillion cubic feet) for all the countries in the world, and the corresponding increase in total reserves with the inclusion of shale resources are 11% and 47%, respectively [1]. These numbers are expected to increase even more in the coming years, and the productions from shale resources are also skyrocketing. With the growing contributions of the hydrocarbon productions from these unconventional plays, maintaining a prolific and sustainable shale resources development is currently a crucial goal in the oil and gas industry. The key factors contributing to the success of shale resources development are horizontal drilling and hydraulic fracturing. A well engineered design of the well spacing that considers the complex and uncertain reservoir and fracture properties is a key point to achieve such a target.
There are several aspects of challenges that need to be tackled for the well spacing optimization problem. For a shale reservoir, well spacing can be complicated by the subsurface uncertainties-associated hydraulic and natural fractures, multiphase flow, different types of interferences between wells such as frac-hits or long-term production/pressure interferences, and economic cost structures [2][3][4][5]. Methods such as rate transient analysis and pressure transient analysis, field pilot tests, and decline curve analysis are all feasible ways to partially solve some of the challenges mentioned above [6][7][8][9][10]. However, these methods are considered either too time expensive or too cost ineffective under the current fast-paced shale development environment. Alternative approaches include the consideration of geomechanics [11] and microseismic analysis [12], but these face the issues of inadequate data, due to high data acquirement costs and extra uncertainty, which is especially observed in the noisy microseismic data. Extensive analytical and numerical simulation approaches have been conducted for the well spacing problem by various authors [13][14][15], since it is relatively cheap. However, these simulations rarely consider the subsurface uncertainties associated with the fractures. Although [16,17] utilize the time of flight and depth of investigation concepts to achieve an analytical solution, it could encounter astronomical difficulty when complex systems of natural fractures are present. Therefore, it is important to understand the implications of different challenges on the problem of well spacing optimization, and crucial to develop an approach that minimizes the cost and maximizes efficiency to solve this integrated problem.
In order to maximize the understanding of the inherited uncertainty associated with the problem of well spacing, a statistical approach is necessary [3]. Some of the methods propose a simple Monte Carlo simulation to tackle the uncertainties [18,19]. However, these methods fail to extract the posterior distributions of the uncertainties so that new, representative samples can be generated to calibrate the reservoir model comprehensively. Other approaches focus on the extensive work of history matching, and the concept of assisted history matching (AHM) has recently become popular. It is confirmed by several works that AHM is the most promising method to understand the effects of uncertainties, and shorten the learning curve for the numerical characterization of shale reservoirs [20][21][22]. In addition, the application of the Markov chain Monte Carlo algorithm coupled with AHM has proven a great potential in the reservoir realization with uncertainty [23,24]. Therefore, in this study, the primary objective is to characterize the optimal well spacing in the shale gas reservoir, using the statistical method/non-deterministic method. Combinations of reservoir/fracture properties, with the lowest global error, from the AHM are gathered [24] as input to the well spacing optimization workflow. These realizations are not randomly sampled, but are rather robustly validated by the actual production history. The applied non-intrusive embedded discrete fracture model (EDFM) can accurately model both hydraulic and natural fractures, which are impossible to be precisely described by the traditional local grid refinement (LGR) method. With the help of the EDFM and a third-party reservoir simulator, 20 years of gas and water production can be predicted and the distributions of the estimated ultimate recovery (EUR) per well for both gas and water phase can be gathered. These results are collected to calculate the economic performance distributions, such as the distributions of the total net present value (NPV), using specific equations that consider the cost and revenues. The box-plot concept in statistics is incorporated to easily determine the optimal well spacing. Finally, the pressure distributions of both the matrix stimulated rock volume (SRV) Energies 2020, 13, 3965 3 of 24 and fractures at selected timelines are plotted to visualize the drawdown characteristics. This entire workflow is then applied to a shale gas reservoir, and optimal well spacing among the designs of different well placement scenarios with the inclusion of natural fractures is obtained. Furthermore, the economic performance of this reservoir was studied and acquired the NPV distributions for all the well spacing scenarios. With this comprehensive workflow, some significant insights regarding the decision-making processes are gained for unconventional developments.

Methodology
The non-intrusive EDFM method is an efficient way to accurately model the flow behavior of complex fracture networks, which was originally developed by [25][26][27]. The schematics are shown below, conceptually demonstrating the mechanism of the non-intrusive EDFM method ( Figure 1). As shown, two fractures are explicitly embedded into three matrix cells in the computational domain and four additional fracture cells were generated in the computational domain to represent these two fractures. Three types of non-neighboring connections (NNCs) were proposed to include the flow connections between the matrix and fractures, between the fracture segments in a single fracture, and between the intersecting fracture segments [25]. These cells contain modified information, such as permeability, porosity, transmissibility, and relative permeability specifically for the direct characterization of these media. Structured matrix cells can be used together with additional fracture cells to easily handle 3D complex fractures. More details of NNC calculations and validations can be found in [25]. One advantage of the non-intrusive EDFM method is that it introduces a smaller number of total cells and NNCs to reduce the computational time when compared to the traditional local grid refinement method. Another advantage is that the non-intrusive EDFM method can be easily applied in any third-party reservoir simulator to model complex fractures.
The following flow chart ( Figure 2) demonstrates the steps of the well spacing optimization workflow. First of all, it was necessary to obtain the accurately characterized combinations of uncertain parameters in the shale gas reservoir that were validated by the known production history. This could be achieved by AHM. These combinations were not generated randomly, but are rather the samples from the true posterior distributions of uncertain parameters after AHM. For more details regarding the AHM, please refer to the work by [24], in which the AHM concept is elaborated in detail. A total of 100 combinations of uncertain parameters, such as the matrix petrophysical properties, hydraulic fractures' geometry and properties, as well as the natural fractures' geometry and properties were selected as input for the well spacing optimization workflow.
Energies 2020, 13, x FOR PEER REVIEW 3 of 23 determine the optimal well spacing. Finally, the pressure distributions of both the matrix stimulated rock volume (SRV) and fractures at selected timelines are plotted to visualize the drawdown characteristics. This entire workflow is then applied to a shale gas reservoir, and optimal well spacing among the designs of different well placement scenarios with the inclusion of natural fractures is obtained. Furthermore, the economic performance of this reservoir was studied and acquired the NPV distributions for all the well spacing scenarios. With this comprehensive workflow, some significant insights regarding the decision-making processes are gained for unconventional developments.

Methodology
The non-intrusive EDFM method is an efficient way to accurately model the flow behavior of complex fracture networks, which was originally developed by [25][26][27]. The schematics are shown below, conceptually demonstrating the mechanism of the non-intrusive EDFM method ( Figure 1). As shown, two fractures are explicitly embedded into three matrix cells in the computational domain and four additional fracture cells were generated in the computational domain to represent these two fractures. Three types of non-neighboring connections (NNCs) were proposed to include the flow connections between the matrix and fractures, between the fracture segments in a single fracture, and between the intersecting fracture segments [25]. These cells contain modified information, such as permeability, porosity, transmissibility, and relative permeability specifically for the direct characterization of these media. Structured matrix cells can be used together with additional fracture cells to easily handle 3D complex fractures. More details of NNC calculations and validations can be found in [25]. One advantage of the non-intrusive EDFM method is that it introduces a smaller number of total cells and NNCs to reduce the computational time when compared to the traditional local grid refinement method. Another advantage is that the non-intrusive EDFM method can be easily applied in any third-party reservoir simulator to model complex fractures. The following flow chart ( Figure 2) demonstrates the steps of the well spacing optimization workflow. First of all, it was necessary to obtain the accurately characterized combinations of uncertain parameters in the shale gas reservoir that were validated by the known production history. This could be achieved by AHM. These combinations were not generated randomly, but are rather the samples from the true posterior distributions of uncertain parameters after AHM. For more details regarding the AHM, please refer to the work by [24], in which the AHM concept is elaborated in detail. A total of 100 combinations of uncertain parameters, such as the matrix petrophysical properties, hydraulic fractures' geometry and properties, as well as the natural fractures' geometry and properties were selected as input for the well spacing optimization workflow.  The following flow chart ( Figure 2) demonstrates the steps of the well spacing optimization workflow. First of all, it was necessary to obtain the accurately characterized combinations of uncertain parameters in the shale gas reservoir that were validated by the known production history. This could be achieved by AHM. These combinations were not generated randomly, but are rather the samples from the true posterior distributions of uncertain parameters after AHM. For more details regarding the AHM, please refer to the work by [24], in which the AHM concept is elaborated in detail. A total of 100 combinations of uncertain parameters, such as the matrix petrophysical properties, hydraulic fractures' geometry and properties, as well as the natural fractures' geometry and properties were selected as input for the well spacing optimization workflow. Once the calibrated combinations of the uncertain parameters from the AHM was obtained, the model size should then be determined for the well spacing optimization. Depending on the size of the data input and also the available computational resources, for the x direction size of the model, one could either divide the model into sectors, or simulate the entire length of the well. In this study, the entire original well length is simulated to accurately depict the whole productivity in the reservoir. For the y direction model size, it is important to consider all the different well placement scenarios that will be studied. For example, if a reservoir model has 1524 m (5000 ft) of y dimension, and a well Energies 2020, 13, 3965 5 of 24 placement from 2 wells to 10 wells, the spacing scenarios will range 762 m (2500 ft) down to 152 m (500 ft). Finally, the z direction of the model was generally determined from the related thickness information regarding the specific reservoir at hand. In this study, well placement scenarios from 2 wells to 6 wells were modelled, and the model size was tuned to give a spacing distance ranging from 152 m (500 ft) to 457 m (1550 ft). After these steps, a clearer picture of the general model information can be gained, and this information was then easily transferred to the proposed workflow. The input files for both the in-house EDFM preprocessor and the commercial simulator will be automatically generated by the workflow, based on a formatted template file. An automatic fracture generator is also developed that can create natural fractures randomly in the given model. It is always good practice to closely examine some input files for each different spacing scenario, so that one can make sure that the input files are generated correctly before the simulation runs. If the input files are error-free, all the simulations are then run in parallel. Different scenarios are divided into balanced workloads for different computation powers, which can greatly improve the efficiency of the workflow. The total number of the parallel runs can be determined based on the available computational resources, up to a total of N runs.
Once the simulation runs were completed, the cumulative productions of gas and water were extracted for all the scenarios for the simulated time period (20 years in this study). Then, the economic measure of NPV can be calculated by the following equations: where NPV n is the total NPV of a specific well spacing scenario, C i is the cost of each month, I i is the gross income of each month and R Discount represents the annual discount rate. Then, each term in (1) is explained in detail: where C well is the cost of the individual well while N well is the total number of wells: where P gas is the gas price and V gas,j is the production of gas of a specific month j: where C water , T gas , T ex are the water disposal cost, gas tax rate and miscellaneous tax rate, respectively. V water,j is the production of water of a specific month j. To make the calculation as accurate and realistic as possible, it was suggested to glean the information regarding these economic parameters for the specific project of interest. The relevant information regarding the cost and gas price will be provided later in the application part. These calculations are repeated for every combination, and for every different well spacing scenario (100 combinations by five spacing scenarios would be a total of 500 cases). Then, all the NPVs calculated within the same spacing scenarios were collected. These five groups of results can be effectively represented as box-plots, which convey information regarding the maximum values, minimum values, P10, P50, and P90 of the NPV. From this generated plot, there usually exists a specific spacing scenario whose values are distributed greater than any other spacing scenario. This scenario then would be the optimal well spacing scenario for the specific reservoir. Finally, the pressure distributions in both matrix and fractures can be extracted and visualized, at any specific time step in the entire 20 years of production, to gain some insights regarding the drainage area and give a rough estimate about drainage efficiency.

Reservoir Model
In this study, the proposed well spacing workflow was applied to a shale gas reservoir of the Sichuan Basin in China. Field-scale compositional reservoir models including multiple horizontal wells and hydraulic and natural fractures were built, as shown in Figure 3. A single component of methane (CH4) was considered. The model dimensions were 1670 m (5480 ft) long and 945 m (3100 ft) wide with a thickness of 20 m (65 ft). There were 137, 31 and 1 grid blocks in the x, y, and z directions, respectively. Gas desorption was modeled using a Langmuir isotherm with a Langmuir inverse-pressure parameter of 0.0775 MPa −1 (0.000534 psi −1 ), maximum moles of adsorbed component per unit rock mass of 0.132 gmol/kg (0.06 gmol/lb), and a bulk density of 2.59 g/cm 3 (161.69 lb/ft 3 ) [28]. The assumption was made that the initial water saturation of the matrix was 39%, and the initial gas saturation was thus 61%. In addition, the residual gas and water saturations were 10% and 20% respectively. The reservoir and fracture information were mainly based on an automatic history matching study with the same actual shale-gas well from the Sichuan Basin in China by [24]. As shown in Figure 3, 2~6 horizontal wells with a well spacing of 158, 189, 236, 315, and 472 m (517, 620, 775, 1033, and 1550 ft), respectively, were investigated for well spacing optimization. The well length for each well was 1500 m (4921 ft). It was assumed that the hydraulic fractures' geometry followed a planar shape and the dip angle was completely 90 degrees. There were 18 stages of hydraulic fractures and three clusters per stage for each well. The distance of the last cluster and first cluster between the two neighboring stages was 44 m (145 ft), and the distance between the two clusters in the same stage was 20.4 m (67 ft). All the designed well spacing scenarios in Figure 3 clearly show that the six-well scenario had overlapping hydraulic fractures in a zipper-frac pattern, implying that there might be stronger well interference than other scenarios. It was also suspected that the two-well scenario had a very large separation between the wells, potentially meaning that the drainage efficiency was not maximized. Therefore, it is important to understand that the well spacing problem has an optimal solution, and neither too much spacing nor too tight spacing is the ultimate answer to this question. What is worth pointing out is that all the hydraulic fractures are assumed to have the same half lengths and heights. The main idea is to simulate equivalent half-length and height values that will be accurately calibrated by the production data. It is true that a hydraulic fracture usually exhibits uneven stimulation growth in the half-length direction (some fractures are longer than average, others lower). This could possibly increase the potential spacing that maximizes the project NPV, because some fractures will grow more and could impose a negative impact on well interference. Therefore, it is important to keep in mind that this study is an idealized scenario.
There are in total two different sets of natural fractures that are generated randomly in the model. The height of the natural fractures is fixed and assumed to be the same as the reservoir thickness of 20 m (65 ft). The strike angles for the first set and second set of natural fractures were 45 degrees to the north and 135 degrees to the north, respectively. The dip angles of both sets were assumed to be 90 degrees, or vertical. The natural fractures' apertures were fixed to 0.03 m (0.1 ft). The reservoir was assumed to be homogenous with a fixed thickness. The detailed reservoir and fracture properties are given in Table 1. It is suspected that some degrees of spatial variations of both hydraulic fractures and natural fractures are usually present in reality. Nevertheless, to accurately describe such a phenomenon, microseismic data, wellbore image logs and some other fracture diagnostic data are required. Without such information regarding the spatial distribution of any kinds of fractures, it would be safest to make the simplifications for both natural and hydraulic fractures.       In order to capture the effect of uncertainty on well spacing optimization, a total of 10 uncertain reservoir and fracture parameters were selected from the work of [24] during the history matching process. They included the matrix permeability, fracture height, fracture half length, fracture conductivity, fracture width, matrix porosity, fracture water saturation, number of natural fractures, natural fractures' length and the natural fractures' conductivity. Their effective ranges are listed in Table 2. One thing to point out here is that the matrix permeability and matrix porosity were assumed to be homogeneous throughout the model (again, the concept of equivalent permeability/porosity was adopted here). Although this is a simplification of the problem, it would be inappropriate to couple a random Gaussian field of such properties to our model. This method lacks realistic physical constraints and would also add another layer of uncertainty to the reservoir model (instead it is best practice to minimize the impact of uncertainties). Therefore, the best method to solve the permeability/porosity uncertainty problem was to use a homogeneous, equivalent model, unless the actual seismic-calibrated geological model was available for this study. A total of 100 samples were selected with the lowest global errors, which were generated by the AHM process (refer to [24] for greater details). To briefly summarize how these 100 samples were generated for the well spacing optimization problem, the following steps are provided. First, 50 initial samples were created based on the latin hypercube sampling method given the empirical range of all the above listed uncertainties. Then, these samples were inputted into the reservoir simulator to obtain the simulated responses, such as the bottomhole pressure, water flow rate, etc. Usually, two of these response parameters are selected for the objective function that measures the error between the simulated response and the real response. Subsequently, a proxy model using any algorithm of neural network was built upon the initial samples to reflect the relationship between the uncertainty parameters and the two response parameters. Furthermore, the Markov chain Monte Carlo (MCMC) algorithm was used to select 25 new and representative samples. These new samples are again inputted into the simulator to get the results, and the proxy models are updated with the added samples. This procedure continued until the maximum number of iterations was achieved. Then, from an additional MCMC process, thousands of representative samples were generated, and the 100 samples with the lowest global errors were filtered out and used for the well spacing optimization problem.

Reservoir Description Value Unit
The relative global errors for these 100 samples range from 19% to 23.5%, which is relatively low considering this as the match with the field data. The collected 100 combinations of uncertainty variables from the history matching solutions were plotted as the form of probability distribution functions, as shown in Figure 4. These plots are best-fitted by an approximating normal distribution, based on the data distribution of the selected 100 samples. The x axis represents the value of the uncertainty variable, while the y axis corresponds to the normalized probability of the uncertainty variable at a specific value. Each of the uncertainty parameters was gathered, sorted and plotted Energies 2020, 13, 3965 10 of 24 as below. These validation distributions are following the distributions obtained from the history matching results (see [24]), and thus they are representative of this shale gas reservoir and can be effectively used to calculate the NPV distributions. The peak of these distributions corresponds to the value with most probable occurrence, or mode. These values, in the order they appear in the figure, are 380 nd, 87 m (286 ft), 13.7 m (45 ft), 37.2 md-m (122 md-ft), 77%, 0.14 m (0.46 ft), 5.3%, 528, 69.2 m (227 ft), and 1.4 md-m (4.6 md-ft), respectively. The modes of some of the listed uncertainty parameters are worth pointing out. For example, the mode of matrix permeability was around the average shale value (shale permeability ranges from 10 to 1000 nd). Then, the mode of the fracture half length was lying in the empirical engineering range, which is from 80 m to 120 m (262 ft to 394 ft). Since the reservoir thickness was only 20 m, and due to the reason that the layering effect was very prominent in this shale formation, the mode for fracture height was reasonable. The number of natural fractures was also considered as reasonable, because the experience from the field data show that the natural fractures system was very developed in this geological formation.
The same flowing bottomhole pressure (BHP) was applied for each well as the reservoir simulation constraint, as shown in Figure 5. During the first 511 days, the flowing BHP was designed to match the provided data. After 511 days of history production, the flowing bottomhole pressure was dropped to 6.9 MPa (1000 psi) up to 20 years to simulate the maximum drawdown. The drawdown was gradually performed instead of a sudden drop, since the BHP is scheduled to follow a linear decline from the last day of the history matching to about 2 years when the BHP finally becomes the targeted 6.9 MPa. Because the rate of this linear decline is still very drastic when compared to the BHP decline encountered in the history period, it was expected that there would be more pressure differential between the wellbore and the reservoir to release more gas and water. Therefore, the flow rates for both phases would increase slightly (if a very gradual linear decline slope was adopted for the BHP schedule, there would not be a significant increase in the flow rates for both phases). The flow rates for both phases reached their relative maximum values at 2 years and then dropped afterwards, because the minimum BHP was reached and the pressure differential would gradually reduce. The two-phase flow of gas and water was simulated using the relative permeability curves as shown in Figure 6b, and it was kept unchanged for all the scenario cases. These curves were based on the history matching solution combination that gives the least global error ( Figure 6a shows the relative permeability results of all the history matching solutions that satisfy the error threshold). The reason to only utilize the best-match results for the relative permeability was that the effect of the relative permeability curve on the EUR predictions was much less pronounced than the hydraulic fracture geometry (such as the fracture height and fracture half length). In addition, the matrix parameters, such as permeability and porosity, also had a bigger influence on the 20-year EUR predictions. We also assumed that the fractures and matrix shared the same Corey type relative permeability curve.   The same flowing bottomhole pressure (BHP) was applied for each well as the reservoir simulation constraint, as shown in Figure 5. During the first 511 days, the flowing BHP was designed to match the provided data. After 511 days of history production, the flowing bottomhole pressure was dropped to 6.9 MPa (1000 psi) up to 20 years to simulate the maximum drawdown. The drawdown was gradually performed instead of a sudden drop, since the BHP is scheduled to follow a linear decline from the last day of the history matching to about 2 years when the BHP finally becomes the targeted 6.9 MPa. Because the rate of this linear decline is still very drastic when compared to the BHP decline encountered in the history period, it was expected that there would be more pressure differential between the wellbore and the reservoir to release more gas and water. Therefore, the flow rates for both phases would increase slightly (if a very gradual linear decline slope was adopted for the BHP schedule, there would not be a significant increase in the flow rates for both phases). The flow rates for both phases reached their relative maximum values at 2 years and then dropped afterwards, because the minimum BHP was reached and the pressure differential would gradually reduce. The two-phase flow of gas and water was simulated using the relative permeability curves as shown in Figure 6b, and it was kept unchanged for all the scenario cases. These curves were based on the history matching solution combination that gives the least global error ( Figure 6a shows the relative permeability results of all the history matching solutions that satisfy the error threshold). The reason to only utilize the best-match results for the relative permeability was that the effect of the relative permeability curve on the EUR predictions was much less pronounced than the hydraulic fracture geometry (such as the fracture height and fracture half length). In addition, the matrix parameters, such as permeability and porosity, also had a bigger influence on the 20-year EUR predictions. We also assumed that the fractures and matrix shared the same Corey type relative permeability curve.

Simulation results
The gas and water flow rates and cumulative gas and water production of all 500 simulation cases for the five different well spacing scenarios are shown in Figure 7. It can be clearly seen that the gas and water flow rates reached their peak values very early in the history. Then, both fluid flow rates reached a maximum at around 2 years, which implies that the drawdown to 6.9 MPa (1000 psi) was becoming effective to increase both the gas flow rate and water flow rate. The range of 20 year cumulative production with complex natural fractures is 212-637 Mm 3 : million cubic meters (7. 5-22.5 Bcf: billion cubic feet) for gas and 4.32-12.8 Km 3 : thousand cubic meters (27-80 MSTB: thousand stock tank barrels) for water. The maximum fluid flow rates for gas and water at about 2 years range from 0.127 to 0.396 Mm 3 /day: million cubic meters per day (4.5-14 MMscf/day: million standard cubic feet per day) for gas and from 2.4 to 11.2 m 3 /day (0.015-0.07 MSTB/day). These ranges are rather wide, providing a comprehensive quantification for the uncertainties in the well spacing optimization problem.

Simulation Results
The gas and water flow rates and cumulative gas and water production of all 500 simulation cases for the five different well spacing scenarios are shown in Figure 7. It can be clearly seen that the gas and water flow rates reached their peak values very early in the history. Then, both fluid flow rates reached a maximum at around 2 years, which implies that the drawdown to 6.9 MPa (1000 psi) was becoming effective to increase both the gas flow rate and water flow rate. The range of 20 year cumulative production with complex natural fractures is 212-637 Mm 3 : million cubic meters (7.5-22.5 Bcf: billion cubic feet) for gas and 4.32-12.8 Km 3 : thousand cubic meters (27-80 MSTB: thousand stock tank barrels) for water. The maximum fluid flow rates for gas and water at about 2 years range from 0.127 to 0.396 Mm 3 /day: million cubic meters per day (4.5-14 MMscf/day: million standard cubic feet per day) for gas and from 2.4 to 11.2 m 3 /day (0.015-0.07 MSTB/day). These ranges are rather wide, providing a comprehensive quantification for the uncertainties in the well spacing optimization problem.
After gathering the 20 year cumulative production data, it is better to further visualize the EUR of gas and water by plotting the lognormal probability plot, as shown in Figures 8 and 9. The lognormal probability plot is a powerful data analysis presentation to visualize distribution that spans a large range, such as the fluid EUR in this study. It combines both the advantages of the logarithm scale and normal distribution scale to provide insights into the P10, P50 and P90 of the distribution. Rough estimates of these values are tabulated in Table 3. The x axis is the interested EUR in the logarithm scale, and the y axis is the cumulative probability of the corresponding observation points, meaning that the percentage of the data that is below this specific observation value. There is usually a linear fit operation performed on the compiled data, to give a rough but continuous extrapolated value that is between the observed values. When comparing the value of both the gas and water EUR to study the effect of the well spacing on recovery, it would be best to use the P50 value (or the median) to make any judgements, since the median value is considered to be more statistically representative because it would be less sensitive to the extreme outliers. The P10 and P90 values would also be less reliable measures, since these situations are statistically less likely and thus less representative scenarios. Table 4 compares the relative change of the both gas and water P50 EUR values when an additional well is added (thus tighter spacing). It is clearly seen that the P50 EUR values for both phases decrease when more wells are present, but the degradation of the EUR values for both phases is less severe when going from a two-well scenario to three-well scenario. This can be explained by the fact that the well interference (production interference) was significantly less severe than when there are more than three wells. Another observation from Table 4 is that the degradation of the EUR values for both phases was also less severe when going from a five-well scenario to a six-well scenario.
This phenomenon implies that the negative impact of well interference was maximized for both of these two scenarios. Looking at (d) and (e) from Figure 3, it is obvious that the spacing between the fracture tips of the neighboring wells is very similar (the fractures from two neighboring wells are almost connected). However, from Figure 3c, there is still some distance between the fracture tips of the two neighboring wells. Therefore, it would also be reasonable to conclude that the effect of well interference on the 20 year two-phase EUR can be optimized. Furthermore, it is never practical to model a well spacing above 472 m (1550 ft), according to field experience, because a large drainage area is not well exploited. This is due to the fact that there is a limit for the lateral growth of the fracture (fracture half length, in this reservoir it rarely goes above 137 m or 450 ft), and well interference (both production and pressure) is therefore negligible for spacing above 472m (1550 ft), disregarding the phase. After gathering the 20 year cumulative production data, it is better to further visualize the EUR of gas and water by plotting the lognormal probability plot, as shown in Figures 8 and 9. The lognormal probability plot is a powerful data analysis presentation to visualize distribution that spans a large range, such as the fluid EUR in this study. It combines both the advantages of the logarithm scale and normal distribution scale to provide insights into the P10, P50 and P90 of the distribution. Rough estimates of these values are tabulated in Table 3. The x axis is the interested EUR in the logarithm scale, and the y axis is the cumulative probability of the corresponding observation points, meaning that the percentage of the data that is below this specific observation value. There is usually a linear fit operation performed on the compiled data, to give a rough but continuous extrapolated value that is between the observed values. When comparing the value of both the gas and water EUR to study the effect of the well spacing on recovery, it would be best to use the P50 Furthermore, it is never practical to model a well spacing above 472 m (1550 ft), according to field experience, because a large drainage area is not well exploited. This is due to the fact that there is a limit for the lateral growth of the fracture (fracture half length, in this reservoir it rarely goes above 137 m or 450 ft), and well interference (both production and pressure) is therefore negligible for spacing above 472m (1550 ft), disregarding the phase.

Economic Analysis
Then, the NPVs of all the runs were calculated based on a well cost of USD 8.0 million, as shown in Figure 10. The gas price was modeled using 1.8 dollars per Mscf. The gas tax rate was set at 6% and the miscellaneous tax rate was set at 5%. Water disposal cost was 0.6 dollars per barrel and the discount rate was 12%. All the NPVs (calculated by the equations in the methodology section) were collected for each scenario and the distributions of this economic performance was presented as box-plots in Figure 10. The box-plot is a classical method to directly and effectively represent any distributions in the data analysis field. The x axis is usually categorical and discrete, representing the scenarios or the intervals that independent variables fit. The y axis is the interested response value and is usually continuous. There are specific meanings to the geometry of the box-plot. The lower bar with a horizontal line is the minimum observed value, the lower edge of the box is the P25 (first quartile) of the observed values, the middle bar in the box is the P50 (median), the upper edge of the box is the P75 (third quartile), and the upper bar with a horizontal line is the maximum observed value. Sometimes, there are points outside the described geometry, which are classified as the outliers. The detailed information regarding the P25, P50, and P75 values is tabulated in Table 5.
From Figure 10, it can be observed that a maximum value of NPV is situated at approximately 236 m (775 ft) of well spacing, by directly using the history matching solutions. The NPV after 20 years ranges from USD 15 million to USD 35 million for this optimal well spacing scenario. Another observation that can be made is that the economic performance is not good for low wells per section scenario, since the NPV for a well spacing of 472 m (1550 ft) was mainly distributed in the interval with lower values (from USD 10 million to USD 23.5 million), with no values above USD 25 million as opposed to all other cases. This observation implies that it is better for the operator to drill at tighter spacing rather than separating wells further away. The random generation of natural fractures might explain this observation, since some of the natural fractures are not effectively connected and thus would not contribute to the drainage. value. Sometimes, there are points outside the described geometry, which are classified as the outliers. The detailed information regarding the P25, P50, and P75 values is tabulated in Table 5.   Figure 10, it can be observed that a maximum value of NPV is situated at approximately 236 m (775 ft) of well spacing, by directly using the history matching solutions. The NPV after 20 years ranges from USD 15 million to USD 35 million for this optimal well spacing scenario. Another observation that can be made is that the economic performance is not good for low wells per section scenario, since the NPV for a well spacing of 472 m (1550 ft) was mainly distributed in the interval with lower values (from USD 10 million to USD 23.5 million), with no values above USD 25 million as opposed to all other cases. This observation implies that it is better for the operator to drill at tighter spacing rather than separating wells further away. The random generation of natural fractures might explain this observation, since some of the natural fractures are not effectively connected and thus would not contribute to the drainage.

Pressure Visualizations
After determining the optimal well spacing of four wells per section with a spacing of 236 m (775 ft), it would be insightful to visualize both the short-term (1 year and 5 years) and long-term (10 years and 20 years) pressure distributions in the matrix and fracture, respectively, as shown in Figures 11 and 12, in which the drainage efficiency can be clearly observed. As can be seen, the natural fractures further away from the wellbore do not contribute to the drainage in the early production period (less than 10 years). However, the effect of natural fractures on the pressure drawdown becomes more pronounced after 20 years of production. In addition, pressure drawdown in the hydraulic fractures is very fast, since it reached bottomhole pressure in less than 5 years, implying

Pressure Visualizations
After determining the optimal well spacing of four wells per section with a spacing of 236 m (775 ft), it would be insightful to visualize both the short-term (1 year and 5 years) and long-term (10 years and 20 years) pressure distributions in the matrix and fracture, respectively, as shown in Figures 11 and 12, in which the drainage efficiency can be clearly observed. As can be seen, the natural fractures further away from the wellbore do not contribute to the drainage in the early production period (less than 10 years). However, the effect of natural fractures on the pressure drawdown becomes more pronounced after 20 years of production. In addition, pressure drawdown in the hydraulic fractures is very fast, since it reached bottomhole pressure in less than 5 years, implying that the hydraulic fractures play a more important role in contributing to production than natural fractures, except the ones that are intersected by the hydraulic fractures near the wellbores.
Energies 2020, 13, x FOR PEER REVIEW 18 of 23 that the hydraulic fractures play a more important role in contributing to production than natural fractures, except the ones that are intersected by the hydraulic fractures near the wellbores.

Conclusions
An innovative and streamlined well spacing optimization workflow was proposed in conjunction with the EDFM coupled with a compositional reservoir simulator to simulate the shale gas production with a two-phase flow, to perform economic calculations, to determine the optimal well spacing from box-plots and to clearly visualize the pressure distributions at the desired time steps. The effects of well interference due to complex natural fractures is well embedded. In addition, the effects of uncertainties of the reservoir is greatly captured quantitatively. The proposed workflow is seamlessly applied to an actual well in the Sichuan Basin of southwest China. Most of the findings in this study are only restricted to the shale gas well located in Sichuan Basin, and might not be applicable to other shale gas formations. The following conclusions can be summarized from this study: (1) Well spacing can be optimized by designing well placement scenarios using multiple history matching solutions and by maximizing the calculated NPV. There exists a range of solutions that can improve the long-term project NPV for the well spacing optimization problem, so landing outside this solution range will jeopardize the economic performance.

Conclusions
An innovative and streamlined well spacing optimization workflow was proposed in conjunction with the EDFM coupled with a compositional reservoir simulator to simulate the shale gas production with a two-phase flow, to perform economic calculations, to determine the optimal well spacing from box-plots and to clearly visualize the pressure distributions at the desired time steps. The effects of well interference due to complex natural fractures is well embedded. In addition, the effects of uncertainties of the reservoir is greatly captured quantitatively. The proposed workflow is seamlessly applied to an actual well in the Sichuan Basin of southwest China. Most of the findings in this study are only restricted to the shale gas well located in Sichuan Basin, and might not be applicable to other shale gas formations. The following conclusions can be summarized from this study: (1) Well spacing can be optimized by designing well placement scenarios using multiple history matching solutions and by maximizing the calculated NPV. There exists a range of solutions that can improve the long-term project NPV for the well spacing optimization problem, so landing outside this solution range will jeopardize the economic performance. (2) The better the reservoir and fractures conditions (permeability, porosity, fracture half length, fracture height, etc.), then the higher the predictive NPV is. It would be substantially more beneficial if additional data can be gathered to build more realistic geological models (such as permeability/porosity map) of the reservoir. (3) The natural fractures that are near the wellbores contribute more to the total EUR and the economic performance than those that are not intersected with hydraulic fractures. However, the 2D modeling of the natural fractures are over simplified solutions, and more data (such as wellbore image log) should be acquired to accurately capture the spatial distributions of the network and to provide more insights regarding the well spacing problem. (4) The well spacing that maximizes the project NPV for this study is found to be approximately 236 m (775 ft), with 20-year gas EUR ranging from 101.9 to 135.9 Mm 3 per well (3.6 to 4.8 Bcf per well) and 20 year water EUR ranging from 2.0 to 2.8 Km 3 per well (12.5 to 17.5 MSTB per well). Nevertheless, more information regarding the hydraulic fractures' geometry (using microseismic event data) might significantly affect the results obtained above.