Optimization of Multiple Hydraulically Fractured Horizontal Wells in Unconventional Gas Reservoirs

Accurate placement of multiple horizontal wells drilled from the same well pad plays a critical role in the successful economical production from unconventional gas reservoirs. However, there are high cost and uncertainty due to many inestimable and uncertain parameters such as reservoir permeability, porosity, fracture spacing, fracture half-length, fracture conductivity, gas desorption, and well spacing. In this paper, we employ response surface methodology to optimize multiple horizontal well placement to maximize Net Present Value (NPV) with numerically modeling multistage hydraulic fractures in combination with economic analysis.This paper demonstrates the accuracy of numericalmodeling ofmultistage hydraulic fractures for actual Barnett Shale production data by considering the gas desorption effect. Six uncertain parameters, such as permeability, porosity, fracture spacing, fracture half-length, fracture conductivity, and distance between two neighboring wells with a reasonable range based on Barnett Shale information, are used to fit a response surface of NPV as the objective function and to finally identify the optimum design under conditions of different gas prices based on NPV maximization. This integrated approach can contribute to obtaining the optimal drainage area around the wells by optimizing well placement and hydraulic fracturing treatment design and provide insight into hydraulic fracture interference between single well and neighboring wells.


Introduction
The combination of horizontal drilling and multistage hydraulic fracturing technology has made possible the current flourishing gas production from shale gas reservoirs in the United States, as well as the global fast growing investment in shale gas exploration and development.Multiple transverse hydraulic fractures are generated when all wellbores are drilled in the direction of the minimum horizontal stress.Maximizing the total stimulated reservoir volume (SRV) plays a major role in successful economic gas production.The unprecedented growth of shale reservoirs has brought a new perspective and focus to the optimization of multiwell placement in the same pad.Drilling multiple horizontal wells from a single pad has increasingly become a common approach for developing shale reservoirs due to significant cost, time, and environmental savings.The surface footprint is reduced greatly by drilling multi-well from the same pad due to minimizing the number of surface locations required while increasing the bottom hole contact of the shale resource [1].Zipper fracturing (zipper-frac) and simultaneous fracturing (simul-frac) [2], where two adjacent horizontal wells are hydraulically fractured alternatingly and simultaneously, respectively, are two commonly used hydraulic fracturing techniques to stimulate multi-well from the same pad.Although hydraulic fractures improve gas production from shale gas wells, the cost of operation is expensive.Long laterals require greater volume of liquids and proppants which contribute to greater cost [3].The well economics is also sensitive to well performance and natural gas price due to higher drilling and completion costs.Therefore, optimizing well parameters such as well number and well distance in conjunction with hydraulic fracture parameters, such as fracture spacing and fracture half-length based on economic analysis, are very important, especially in the current environment of low natural gas prices.
Optimization of multi-well placement is primarily valuable for overall project economic viability and minimizing  the risks of well collision in shale gas reservoirs.Closer well placement will result in stimulated reservoir volume intersects, leading to well competition and penalizing overall production [5].However, only recently there have been limited studies of optimizing fracturing design together with multiwell placement simultaneously.Esmaili et al. [6] defined three types of horizontal well based on the number of neighboring wells for sharing drainage area when drilling multiple wells from a pad.Rafiee et al. [7] proposed a new design for two horizontal wells by the modification of the traditional zipperfrac and demonstrated this new design maximized reservoir contact and improved well performance when compared to the original zipper-frac design from both rock mechanics and fluid production aspects.However, such a design assumed a larger fracture spacing of 500 ft and did not optimize well placement and fracture spacing simultaneously to obtain the optimal design for economic gas production.Díaz de Souza et al. [8] did sensitivity studies of three different wells placement with 2, 3, and 4 horizontal wells within the same stimulation volume in the Haynesville Shale to obtain the optimal well spacing.They stated that four horizontal wells with a well distance of 660 ft was a near-optimal solution for this reservoir.However, critical parameters for developing a play economically, such as fracture spacing and fracture half-length, have not been considered for optimization based on optimal well spacing.Harpel et al. [9] reported that well spacing becomes tighter in parts of the Fayetteville Shale from 600 ft to 400 ft or 300 ft, leading to fluid and proppant volume reductions, while further optimization of stimulation design, especially fracture half-length, is very much required for future development.Ramakrishnan et al. [10] suggested that it is particularly challenging to optimize the stimulation treatment in a multi-well design drilled with few hundreds feet of spacing between one another to maximize coverage around each horizontal well, because the optimization process involves perforation and stage placements, sequential stimulation of these wells, fluid and proppant schedules, treatment rates, and application of diversion technology in an effort to achieve effective stimulation along these wells and between the wells.Hence, a detailed study and a comprehensive approach for optimization of fracturing design and multi-well placement are still significantly necessary.
In this paper, we employed response surface methodology to build the response surface in terms of NPV with six parameters such as reservoir porosity, permeability, fracture half-length, fracture conductivity, fracture spacing, and well distance from Barnett Shale, to obtain the best economic scenario for a given range of these influential parameters.The effect of gas desorption is integrated in the numerical modeling of multistage hydraulic fractures.The impact of different gas prices is also taken into account for the optimization process.The goal of this work is to provide insights into the effective exploitation of shale gas reservoirs via optimization of the fracturing design and multiple wells placement simultaneously.

Shale Gas Reservoir Modeling
Given the complex nature of hydraulic fracture growth and the very low permeability of the matrix rock in shale gas reservoirs, coupled with the predominance of horizontal completions, reservoir simulation is the preferred method to predict and evaluate well performance [11][12][13].Local grid

3.6E + 06
Field data Gas rate without gas desorption Gas rate with gas desorption Gas production without gas desorption Gas production with gas desorption (b) Gas production for a 30-year period Figure 4: History matching of Barnett Shale with and without gas desorption effect.

Well 1
Well 2 refinement with logarithmic cell spacing is used in the simulation to accurately model flow from the shale to the fracture, that is, properly incorporate the transient flow behavior from the matrix to the fracture.In a block, the hydraulic fracture is explicitly modeled; moreover, the matrix is described as some subcells whose size increases logarithmically, while moving away from the hydraulic fracture to properly simulate the large pressure drop between the matrix and the fracture.In addition, a dual permeability grid is used to allow simultaneous matrix-to-matrix and fracture-to-fracture flows.This method can accurately and efficiently model transient gas production from hydraulic fractures of the horizontal wells in shale gas reservoirs [14,15].The reservoir is assumed to be homogeneous and the fractures evenly spaced, with stressindependent porosity and permeability.It is assumed that there is no water flow in the reservoir modeling of shale gas.In our simulation, gas is only flowing into the wellbore through the hydraulic fractures, that is, no matrix-wellbore communication.The turbulent gas flow due to high gas flow rate in hydraulic fractures is modeled as non-Darcy flow.The non-Darcy Beta factor, used in the Forchheimer number, is determined using a correlation proposed by Evans and Civan [16] as follows: where the unit of  is md and the unit of  is ft −1 .The  () correlation was obtained using over 180 data points including those for propped fractures and was found to match the data very well with the correlation coefficient of 0.974 [14].This equation is implemented into the numerical model and used for accounting for non-Darcy flow in hydraulic fractures.Figure 1 is a diagram of typical shale gas completion design with a multistage hydraulic fracture treatment, which illustrates several important geometric fracture parameters, such as outer fractures, inner fractures, fracture spacing, and fracture half-length.

Economic Model
NPV is one of the most common methods used to evaluate the economic viability of investing a project.It is referred to   the sum of all cash flows discounted to a specific point in time at the investor's minimum discount rate.The correlation between the present value  and the future value  is where  is the currency escalation rate or interest rate;  is the number of periods.
The NPV is calculated using the following expression: where   is future value of production revenue for a fracture reservoir;   is future value of production revenue for an unfractured reservoir; FC is the total fixed cost;  well is the cost of one horizontal well;  fracture is the cost of hydraulic fracture in a horizontal well;  is the number of horizontal wells.The costs of well and fracture used in the economic analysis are based on the work of Schweitzer and Bilgesu [17], as shown in Table 1.

Langmuir Isotherm
Gas shales are organic-rich formations.Gas storage in the shale is mainly divided into free gas in natural fractures and matrix pore structure and adsorbed gas in organic materials.Langmuir isotherm is widely used to describe the gas adsorption phenomenon.The amount of gas stored in shale is often described by Langmuir equation: where   is the gas content in scf/ton,   is the Langmuir volume in scf/ton,   is the Langmuir pressure in psi, and  is pressure in psi.The bulk density of shale (  ) is needed to convert the typical gas content in scf/ft 3 to scf/ton.Langmuir pressure and Langmuir volume are two key parameters.Langmuir volume is referred to as the gas volume at the infinite pressure representing the maximum storage capacity for gas; Langmuir pressure is referred to as the pressure corresponding to one-half Langmuir volume.As the reservoir pressure is decreased, gas is desorbed from the surface of the matrix.Figure 2 shows a graph of gas content with pressure for the adsorbed gas and total gas used for Barnett Shale [4].Both free gas and adsorbed gas add together to generate the total gas content.In Barnett Shale, the adsorbed gas is approximately 46% of the total gas.Contrary to conventional gas reservoirs, the amount of gas desorption in the matrix is commonly described by the Langmuir isotherm in a range of reservoir pressures.The Langmuir isotherm of the Barnett Shale used in this study is illustrated in Figure 3.It is clearly shown that higher Langmuir pressure releases more adsorbed gas and results in higher gas production.Generally, in early stage of production, when reservoir pressure is high, the gas desorption contribution to the gas production is insignificant; however, for long-term production, it is necessary to account for gas desorption, based on a laboratory measured isotherm due to the more substantial pressure depletion, resulting in more gas desorption.CMG [18] was used to model the effect of gas desorption from a shale gas reservoir in a black oil model with a technique developed by Seidle and Arri [19].A Langmuir isotherm is replicated by a black oil model's solution gas ratio to include the effect of gas desorption in shale.

History Matching for Barnett Shale
Published average reservoir data for a Barnett Shale well were used for history matching [20].In this case, the well was stimulated by a multistage fracturing with a single, perforated interval for each stage.In this simulation study, we set up a reservoir with a volume of 3000 ft × 1500 ft × 300 ft.The fracture spacing and half-length are set at 100 ft and 150 ft, respectively, and the number of fractures is 28.Detailed reservoir information about this section of the Barnett Shale is listed in Table 2.The reservoir is assumed to be homogeneous and the fractures evenly spaced, with stress-independent porosity and permeability.Only gas is flowing in the reservoir, which is assumed to behave as non-Darcy flow.The history matching of field data is presented in Figure 4(a).It shows a more reasonable match between the numerical simulation results and the actual field gas flow data, considering the effect of gas desorption, contributing to 15.6% of total gas production at around 4.5 years of gas production.In addition, Figure 4(b) shows the forecasting of gas production for a 30year period with and without considering gas desorption.As shown, with the gas production, the gas desorption contributes more due to substantial pressure depletion and larger gas drainage area and finally contributes to 20.7% of the total gas production at 30 years of gas production.Thus, the impact of gas desorption cannot be ignored when performing history matching and assessing production forecast of gas production in Barnett Shale formation.Hence, this study takes into account gas desorption effect for the subsequent optimization of multiwell placement in Barnett Shale.

Multiwell Modeling
Two scenarios describing multiple horizontal well placement were studied, as illustrated in Figure 5. Scenario 1 is referred to as aligning fracturing, where hydraulic fracturing is between two wells in an aligned pattern, and Scenario 2 is referred to as alternating fracturing, where hydraulic fracturing is between two wells in a staggered pattern.We investigate the effect of fracture spacing towards comparison of gas production between these two scenarios.We set up a shale gas reservoir model with a volume of 5000 ft × 1600 ft × 200 ft.Detailed reservoir information used for Barnett Shale is listed in Table 3.Comparison of cumulative gas production is shown in Figure 6.It shows that there is almost no difference of gas production for these two scenarios when fracture spacing is below 400 ft; however, Scenario 2 yields higher gas production than Scenario 1 when the fracture spacing is 400 ft or above.Figure 7 presents the comparison of average reservoir pressure with fracture spacing of 600 ft.It can be seen that for Scenario 2, it has a larger average reservoir pressure drawdown, leading to higher cumulative gas production.Pressure distribution at    30 years of gas production for these two scenarios is shown in Figure 8.It shows that more contact reservoir area is drained effectively by Scenario 2, compared to Scenario 1.In addition, Rafiee et al. [7] reported that Scenario 2 design can increase the stress interference between fractures and create more effective stimulated reservoir volume to improve gas production.Therefore, Scenario 2 is used for optimization of multi-well placement in the subsequent study of this paper.

Multiwell Optimization
Response surface methodology (RSM) approach is applied to the optimization of two horizontal well placement based on the Barnett Shale reservoir information.RSM is utilized to approximate a response, in terms of maximum NPV in this work, over a region of interest specified by the range of variability of input factors based on the least squares criterion.The RSM model can be linear or fully quadratic.It can offer a cost-effective and efficient way to manage the uncertainties for shale gas reservoir development.More detailed mathematical and statistical theories of RSM can be found in the work by Myers and Montgomery [21].We set up a shale gas reservoir model with a volume of 5000 ft × 2000 ft × 200 ft.Detailed reservoir information used for the Barnett Shale is listed in Table 4. Six uncertain parameters such as porosity (), permeability (), fracture half-length (), fracture conductivity (), fracture spacing (), and well distance () are given a reasonable range with the actual maximum and minimum values or coded symbol of "+1" and "−1, " respectively, as shown in Table 5.The number of hydraulic fractures is 87 and 35, corresponding to fracture spacing of 40 ft and 100 ft, respectively.These uncertainty ranges come from field data, analogues, and history matching of the Barnett Shale reservoirs.According to 6 variables, 38 cases were required, based on the approach of D-Optimal design, which was originated from the optimal design theory [22].Table 6 shows the 38 combinations of these uncertain parameters generated by the D-Optimal design.After numerical simulation of each case, cumulative gas production and gas flow rate are obtained and shown in Figures 9 and 10 for selected cases with various fracture spacing and fracture halflength.It clearly shows that the cumulative gas production at 30 years of gas production is in the range of 3934-6529 MMSCF and gas flow rate has a large uncertainty.This means further optimization is needed.
Once the cumulative gas production of 38 run cases obtained, the economic Excel spreadsheet is used to calculate the corresponding NPVs.In this work, the impact of different gas price is considered for calculation of NPVs and optimization of fracturing design in multi-well placement to determine the optimal stimulation design.The price of natural gas plays a critical role in economic success of shale gas development.The activity in shale gas exploitation increases with the increasing gas price.Negative NPV of the project suggests that an unoptimized shale gas production can result in not only no income but also loss of money.For gas price of $3/MSCF, $4/MSCF, and $5/MSCF, the positive NPV is in the range of $0.36-13.39 million, $2.10-18.88million, and $3.42-24.38 million, respectively, as shown in Figure 11.
Once NPVs of 38 run cases obtained, the Design-Expert Software is used to build the NPV response surface model in this study.To select the appropriate model, the statistical approach was used to determine which polynomial fits the equation among linear model, two-factor interaction model (2FI), quadratic model, and cubic model, as shown in Tables 7, 8, and 9, for different gas prices.The criterion for selecting the appropriate model is choosing the highest polynomial model, where the additional terms are significant and the model is not aliased.Although the cubic model is the highest polynomial model, it is not selected because it is aliased.Aliasing is a result of reducing the number of experimental runs.When it occurs, several groups of effects are combined into one group and the most significant effect in the group is used to represent the effect of the group.Essentially, it is important that the model is not aliased.In addition, other criteria are to select the model that has the maximum "Adjusted -Squared" and "Predicted -Squared".Thus, the fully quadratic model is selected to build the NPV response surface in the subsequent optimization process.
The equation fitted to the NPV response surface with the coded symbol for different gas prices is presented below.
For gas price of $  ( where  is formation porosity;  is formation permeability, md;  is fracture half-length, ft;  is fracture conductivity, md-ft;  is fracture spacing, ft;  is well distance, ft.The normal plot of residuals, reflecting the distribution of the residuals, for different gas prices is shown in Figure 12.All the points in the "Normal Plot of Residuals" fall on the straight line, meaning the residuals are normally distributed.Figure 13 shows the plot of "Predicted versus Actual" for different gas prices, illustrating whether the generated equation of NPV response surface accurately predicts the actual NPV values.It can be seen that generated NPV response surface models provide such reliable predicted values of NPV, as compared with the actual values of NPV.This means that the generated NPV response surface models are reliable.
Figure 14 shows the 3D surface of the well distance and fracture spacing with gas price of $4/MSCF.It shows that there exists an optimal combination between well distance and fracture spacing.For fracture half-length of 200 ft, the NPV first increases and then decrease with increasing well distance.For fracture half-length of 400 ft, the NPV increases with increasing well distance within the range of this study.With the increasing fracture half-length, the optimal point moves to larger well distance.Figure 15 presents the 3D surface of the well distance and fracture half-length with gas price of $4/MSCF.It can be seen that larger well distance and fracture half-length will lead to higher NPV.Similarly, Figure 16 presents the 3D surface of the fracture spacing and fracture half-length with gas price of $4/MSCF.It shows that there exists an optimal value for fracture spacing.With the increasing well distance, the optimal point moves to larger fracture half-length and smaller fracture spacing.Therefore, it can provide some insights into optimization of stimulation designs and completion strategies to obtain the maximum economic viability of the field.
The numerical optimization option selects the set of variables that leads to the maximum NPV value.In this study for the Barnett Shale, the optimal designs with porosity of 0.06 and permeability of 0.0001 md and different gas prices are listed in Table 10.The optimal NPV value is $8.70 MM, $12.40 MM, and $16.34 MM, for gas price of $3/MSCF, $4/MSCF, and $5/MSCF, respectively.What is more, the optimal fracture spacing is 80 ft for gas price of $3/MSCF, 70 ft for gas price of $4/MSCF, and 60 ft for gas price of $5/MSCF, as shown in Figure 17.It is extremely vital to validate the optimal results.Verification is performed by running the simulation with the best design condition.Figure 18 shows the cumulative gas production and gas flow rate for these three optimum cases with different gas prices, respectively.The NPV is calculated as $7.91 MM, $11.81 MM, and $15.75 MM, for gas price of $3/MSCF, $4/MSCF, and $5/MSCF, respectively.The absolute NPV difference between the NPV from the response surface and the real NPV is small.As indicated by the small relative error, all three solutions show a very good agreement between the calculated NPV and the real NPV.

Summary and Conclusions
The economic success of shale gas reservoirs depends on optimization of the number of treatment stages and number of fractures and horizontal wells.In this paper, response surface methodology was used to obtain the optimal design for shale gas production by optimizing reservoir and fracture uncertainty parameters in two horizontal wells.We applied this method to optimize 6 uncertain parameters, such as permeability, porosity, fracture spacing, fracture half-length, fracture conductivity, and well distance for the Barnett Shale development.Also, the gas desorption effect is considered for modeling shale gas production because of its large contribution to the estimated ultimate recovery for the Barnett Shale.The impact of varying natural gas price is taken into account during the optimization process.The following conclusions can be drawn from this study.
(1) With a porosity of 0.06, a permeability of 0.0001 md, and a fracture conductivity of 26 md-ft, the optimal design combinations for Barnett Shale is fracture half-length of 400 ft, well distance of 1000 ft, and fracture spacing of 80 ft, 70 ft, and 60 ft for gas price of $3/MSCF, $4/MSCF, and $5/MSCF, respectively.
(2) The gas desorption contributes to 20.7% of EUR at 30 years of gas production for an actual Barnett Shale horizontal well.
(3) The proposed approach is practical and efficient for the design and optimization of hydraulic fracturing for multiple horizontal wells in shale gas reservoirs.Gas content, scf/ton   :

Nomenclature
Bulk density of shale, g/cm 3 .

FractureFigure 1 :
Figure 1: A sketch of multiple hydraulic fractured horizontal shale gas well.

Figure 5 :
Figure 5: Two scenarios of multiple horizontal well placement.

Figure 7 :
Figure 7: Comparison of average reservoir pressure.

Case 7 :Figure 10 :
Figure 10: Gas flow rate of 30 run cases for a 30-year period.

Figure 11 :
Figure 11: NPVs of 38 simulation cases at 30 years of gas production with different gas prices.

CEFigure 14 :FFFigure 15 :
Figure 14: 3D surface of NPV at varied values of well distance and fracture spacing with gas price of $4/MSCF.

CDesignCFigure 16 :
Figure 16: 3D surface of NPV at varied values of fracture spacing and fracture half-length with gas price of $4/MSCF.

Figure 17 :
Figure 17: The optimum fracture spacing with different gas prices.

Figure 18 :
Figure 18: Three optimum cases with different fracture spacing for different gas prices.

Table 1 :
Economic data for NPV calculation.

Table 2 :
Parameters used in history matching.

Table 3 :
Parameters used in multiwell modeling.

Table 4 :
Parameters used in multiwell optimization.

Table 5 :
Uncertainty parameters in this study.

Table 7 :
Statistical approach to select the RSM model with gas price of $3/MSCF.

Table 8 :
Statistical approach to select the RSM model with gas price of $4/MSCF.

Table 9 :
Statistical approach to select the RSM model with gas price of $5/MSCF.

Table 10 :
Optimal combinations and optimization validation.