Development of Production-Forecasting Model Based on the Characteristics of Production Decline Analysis Using the Reservoir and Hydraulic Fracture Parameters in Montney Shale Gas Reservoir, Canada

This study developed a production-forecasting model to replace the numerical simulation and the decline curve analysis using reservoir and hydraulic fracture data in Montney shale gas reservoir, Canada. A shale-gas production curve can be generated if some of the decline parameters such as a peak rate, a decline rate, and a decline exponent are properly estimated based on reservoir and hydraulic fracturing parameters. The production-forecasting model was developed to estimate five decline parameters of a modified hyperbolic decline by using significant reservoir and hydraulic fracture parameters which are derived through the simulation experiments designed by design of experiments and statistical analysis: (1) initial peak rate (Phyp), (2) hyperbolic decline rate (Dhyp), (3) hyperbolic decline exponent (bhyp), (4) transition time (T transition), and (5) exponential decline rate (Dexp). Total eight reservoir and hydraulic fracture parameters were selected as significant parameters on five decline parameters from the results of multivariate analysis of variance among 11 reservoir and hydraulic fracture parameters. The models based on the significant parameters had high predicted R2 values on the cumulative production. The validation results on the 1-, 5-, 10-, and 30-year cumulative production data obtained by the simulation showed a good agreement: R2 > 0:89. The developed production-forecasting model can be also applied for the history matching. The mean absolute percentage error on history matching was 5.28% and 6.23% for the forecasting model and numerical simulator, respectively. Therefore, the results from this study can be applied to substitute numerical simulations for the shale reservoirs which have similar properties with the Montney shale gas reservoir.


Introduction
The prediction of shale gas production is important for evaluating the feasibility of shale gas developments. Some methods have been proposed to predict shale gas production which typically shows an early peak production and sharp decline production trend. The Arps method [1] is the most widely used decline curve analysis technique for estimating the ultimate recovery of oil and gas. Table 1 shows the three forms of the Arps equation, in which "q," "q i ," "b," "D i ," and "t" are the mean production rate, initial production rate, decline exponent, initial decline rate, and production time, respectively. Robertson [2] introduced the modified hyperbolic decline formula to avoid overestimating the estimated ultimate recovery. Ilk et al. [3] also combined hyperbolic and exponential decline formulas to predict shale gas production with different flow regimes over time. In addition, Duong [4] introduced an empirical decline model for fracture-dominated shale reservoirs. These methods, however, cannot forecast production without production history data. Numerical reservoir simulation is an effective technique to predict oil and gas production. However, the simulation is time-consuming, particularly for simulating a huge and complicated reservoir model. Therefore, it is necessary to develop a prediction model which can predict shale gas production without running a reservoir simulation.
Some studies proposed prediction models for the shale gas production and net present value (NPV) [5][6][7][8][9]. The prediction models were developed using various analysis techniques and reservoir simulation data. Yu and Sepehrnoori [5] performed a study on the response surface methodology (RSM) to optimize eight parameters (porosity, permeability, reservoir thickness, reservoir pressure, bottom-hole pressure, fracture spacing, fracture half-length, and fracture conductivity) based on the Barnett shale reservoir for predicting the NPV. Kim et al. [6] developed the neural network-based proxy model for predicting shale gas production using the reservoir and fracture design parameters. Nguyen-Le and Shin [7] presented a framework for the development of an economic indicator of shale gas projects based on the reservoir parameters, hydraulic fracturing parameters, and gas price scenarios. Furthermore, Nguyen-Le et al. [8,9] developed long-term shale gas production models using early production data in the Barnett shale reservoir. On the other hand, few studies have examined the effect of the reservoir and hydraulic fracture parameters on the decline characteristics of shale gas production. The decline parameters of shale gas production need to be evaluated because they can be used to characterize the shale-gas production curve.
A general shale-gas production curve shows a high peak production rate followed by a steep initial decline trend [9,10]. The decline curve analysis for shale gas production using the Arps model showed that the decline exponent starts at 4.0 in the first few days, followed by a moderate value in several weeks or months (b = 2:0) and finally becomes zero which follows the exponential decline model [11,12]. Therefore, the use of a single decline curve model with a constant decline exponent may not be suitable for shale-gas production decline analysis.
In this study, a modified hyperbolic decline curve analysis technique was employed to model shale-gas production decline, in which a single hyperbolic decline model was used to describe the early production period and a single exponential model was used to describe the latter production decline trend. Prediction models were then developed to predict five decline parameters including (1) initial peak rate (P hyp ), (2) hyperbolic decline rate (D hyp ), (3) hyperbolic decline exponent (b hyp ), (4) transition time (T transition ), and (5) exponential decline rate (D exp ), from the simulation results based on the Montney reservoir. Among the independent parameters evaluated in the statistical analysis, the reservoir parameters included the following: reservoir thickness (H), reservoir pressure (P i ), matrix porosity (Φ m ), matrix permeability (K m ), matrix initial water saturation (S wim ), Langmuir pressure (P L ), and Langmuir volume (V L ). The hydraulic fracture parameters included the main fracture half-length (X mf ), main fracture conductivity (C mf ), fracture initial water saturation (S wif ), and the number of clusters per stage (N cluster ).

Methodology
This study is aimed at forecasting the shale-gas production curve by estimating five decline parameters of the modified hyperbolic decline through the regression models based on the reservoir and hydraulic fracture parameters. Therefore, it is important to select the significant parameters on five decline parameters and develop regression models for predicting the decline parameters. In multivariate analysis of variance (MANOVA), parameters with p value <0.05 are considered significant parameters.
This study consists of three parts: selection of the significant parameter, development of the forecasting model, and model validation ( Figure 1). The first part concentrated on selecting the significant parameters on the shale-gas production decline. The simulation was performed based on Plackett-Burman design, which is a two-level fractional factorial design [13]. In the second part, the productionforecasting model was developed based on the significant parameters selected. Finally, the developed model was validated by comparing it with the simulation results. In both the second and third parts, Latin hypercube sampling was used to design the simulation experiments.
The Western Canada Sedimentary Basin (WCSB) accounted for about 87 percent of crude oil and 94 percent of natural gas production in Canada [14,15]. There are 11 major formations in the WCSB, and the Montney formation was first developed [14]. The average depth of the Montney shale is 10,000 ft (4500~12500 ft), and the thickness of the formation is between 170 ft and 680 ft [16,17]. Dry gas is mainly produced in the deep buried southwest, with wet gas, gas condensate, and oil produced more eastward [18]. Table 2 lists the range of the input parameters, which are the properties of the Montney reservoir, and Table 3 presents the five decline parameters that are the output parameters. The decline parameters were obtained from the results of the modified hyperbolic decline analysis on the reservoir Table 1: Three forms of the Arps equation [1].  2 Geofluids simulation. P hyp and T transition were the first to be derived. P hyp is the initial peak point of the shale-gas production curve. T transition is the time when the hyperbolic decline changes to the exponential decline, which is the straight line in the semilog graph [3]. The values of D hyp , b hyp , and D exp were then derived from the random numbers that minimize the mean absolute percentage error (MAPE) between the reservoir simulation results and the modified hyperbolic decline equation.

Reservoir Simulation Model.
A 3D fractured reservoir model was constructed using the GEM simulator of the Computer Modeling Group [27]. A dual porositypermeability model was employed: a matrix system and a fracture system. Furthermore, the simulation of a shale gas reservoir needs hydraulic fracture models based on a local grid refinement. A local grid refinement is a method to calculate the flow and matrix behavior near the hydraulic fracture zone accurately. The average data for the base reservoir model, completion design, and operating condition were selected from the literature data on Montney shale reservoir (Table 4). Figure 2 shows the 3D reservoir and hydraulic fracture model. The symmetry models were used for the reservoir simulation, which decrease the simulation run time considerably (Figures 3 and 4). The length, width, and thickness of the symmetry model were 5500 ft (I direction), 525 ft (J direction), and 170 ft (K direction).

Selection of Significant Parameters.
For the eleven parameters in Table 2, a total of 2048 (2 11 ) simulations are required to investigate the full factorial design response, including the main effects and the interaction effects of the parameters. Therefore, the design of experiments needs to reduce the number of total simulations runs. Only 32 simulations were required to evaluate the effects of the input   3 Geofluids parameters by applying the Plackett-Burman design [13]. In addition, it is an efficient method for screening the less influential parameters from the input parameters on the output parameters and simplifying the setting for further simulations. The significant parameters were selected from the results of multivariate analysis of variance (MANOVA) [5,[28][29][30]. Table 5 lists the p values from MANOVA. The p value is an indicator to decide the relative significance between parameters. The criterion for selecting the significant parameters was a significance level of 0.05 [5,[28][29][30]; the significant parameters are highlighted in bold in Table 5.
The significant parameters were used in the simulation experiments designed by Latin hypercube sampling to develop regression models for predicting the decline parameters.   [31]. The system can evaluate all the response space without replications. The Latin hypercube sampling designed in CMOST guarantees two important characteristics. One is the approximate orthogonality of the input parameters, and the other is space filling. That means that all selected samples are statistically independent and evenly distributed within the range [32]. To develop the regression models, we simulated 100 cases designed by CMOST based on the Latin hypercube sampling. Five regression models were employed to find the fittest estimation model based on the predicted R 2 , which means how well the models can predict the responses on the new values of the input parameters (Table 6). In Table 6, the values of all the parameters, which are β 1 , β 2 , ⋯, β k , were calculated to find the best fit models. To develop effective regression models, the backward elimination method removed the term with the highest p value continuously until only statistically significant variables existed in the models 1, 2, 3, 4, and 5. The results of the regression analysis are presented, and the models with the highest predicted R 2 are highlighted in bold Table 7. Model 5 was the most suitable model to predict the decline Table 6: Multiple regression formulas applied for the regression models of production decline analysis.

Multiple regression formulas
Linear ln y = β 0 + β 1 ln x 1 + β 2 ln x 2 + ⋯+β k ln x k +∑∑β ij ln x i ln x j ⟶y = e β 0 +β 1 ln x 1 +β 2 ln x 2 +⋯+β k ln x k +∑∑β ij ln x i ln x j Formula 5 5 Geofluids parameters based on the reservoir and hydraulic fracture parameters except for b hyp . The models were not overfitted because there was a little difference between the actual R 2 and predicted R 2 . Furthermore, the highlighted models also had the highest adjusted R 2 . This means that the models are remarkable compared to others, regardless of the number of terms; Table 8 lists the equations of five regression models for five decline parameters. Figure 5 presents the predicted production rate and cumulative production, which were characterized using the five regression models of decline parameters (P hyp , D hyp , b hyp , T transition , and D exp ) debased on the significant reservoir and hydraulic fracture parameters of the Montney shale gas reservoir with the average values of the reservoir and hydraulic fracture parameters (Table 4). P hyp and T transition determine the initial peak rate and transition time. D hyp , b hyp , and D exp determine the overall decline trend of production curve between P hyp and T transition and after T transition .
The accuracy of the production-forecasting model was examined by validating the developed models; the validation results are presented in this section. We performed 23 simulations to make the validation data on the productionforecasting model. In particular, the accuracy of cumulative production within five years is important because the economic feasibility of horizontal wells to produce the shale gas is strongly related to the net present value of the first five years rather than the ultimate recovery of the well [33]. Figure 6 presents the validation results of five regression models for predicting the decline parameters. The R 2 values for the regression models were between 0.48 and 0.99. Figure 7 shows the prediction results of cumulative production in 1, 5, 10, and 30 years.

Applications of Production-Forecasting Models.
The developed production-forecasting model can be applied to predict shale gas production without a numerical simulation which requires a considerable effort and time to build reservoir model and reservoir simulation if the average reservoir properties are known.
To evaluate the applicability of production-forecasting models, it was validated by comparing the results of history matching on the field data between the proposed model and numerical simulator (CMOST), and the validation results are presented in this section.
The field data was obtained from the horizontal well of Sunset field, Montney, Canada. Table 9 lists the reservoir and completion data of the well. The field production data was adjusted to remove the noise and calibrate the horizontal well length because the production-forecasting models were developed based on the 5500 ft horizontal well with 10 stages.
The history matching process was conducted to find the best combination of matching parameters (H, X mf , S wif , and C mf ) to achieve the minimum error between gas production of field data and numerical simulation. The history matching algorithms used in this study are the designed exploration  7 Geofluids and controlled evolution (DECE) in CMOST, and the distribution of matching parameters (H, X mf , S wif , and C mf ) was assumed to have uniform distribution. The total number of numerical simulations for history matching was 500. Figure 8 shows the history matching results (optimal case) using the CMOST.
A random data set with 10,000 cases was made for history matching using the forecasting model, and the ranges of the random data correspond to the properties of the Montney reservoir. The fittest case was selected from 10,000 cases based on the MAPE, and MAPE of the fittest matching result was 5.28% (forecasting model). On the other hand, MAPE in the numerical simulator was 6.23% (CMOST) as shown in Figure 9.
Numerical simulation consumed three days and 12 hours for competing 500 numerical simulation runs (CMOST), but the production-forecasting model does not require any simulation time. Therefore, the productionforecasting model could be applied to replace the numerical simulation for shale gas production cases.

Discussion
The significant parameters were presented to develop a production-forecasting model for the Montney shale gas reservoir. The developed production-forecasting model showed good agreement with the cumulative production of shale gas.   Several studies have examined the effects of the reservoir and hydraulic fracture parameters on shale gas production. Warpinski et al. [34] reported the effects of the reservoir permeability, fracture spacing, fracture conductivity, stimulated reservoir volume, and cleanup on the production performance and the response of the reservoir pressure. In this study, the initial water saturation of a fracture was also evaluated to consider the flowback effects on the characteristics of shale-gas production decline. The initial water saturation of a fracture in the shale gas reservoir had a significant effect on the early decline rate of the shale gas production. Kim et al. [6] evaluated 11 key parameters for shale gas production. But the significant parameters were not selected among the eleven key parameters. Therefore, the relatively insignificant parameters on shale gas production could be used as independent parameters in the prediction model.
The present study not only investigated the effects of the reservoir and hydraulic fracture parameters but also  Figure 8: Comparison of (a) cumulative gas production and (b) gas rate between gas production of field data and numerical simulation. Gas rate (MMscf/d) Prediction model Field data CMOST Figure 9: Results of the fittest history matching on the field data from the numerical simulator (CMOST) and prediction model. 9 Geofluids identified the relative significance of key parameters. The prediction accuracy of cumulative production is strongly related to the prediction accuracy of the hyperbolic decline characteristics because the production performance of shale gas is dominated by the period of the transient flow regime. This can be proven from the results of sensitivity analysis on the characteristics of production decline analysis. Figure 10 presents the results of sensitivity analysis on the characteristics of production decline analysis based on the mean absolute percentage error on the cumulative  production between the production-forecasting models and numerical simulation. The values estimated from each decline parameter model were marked as a "base" in each graph. The base values are changed to both above and below up to 50%. Sensitivity analysis shows that the transition time and exponential decline rate had little or no effect on the mean absolute percentage error. Therefore, forecasting the transition time and exponential decline rate accurately is not critical for a prediction of shale gas production.

Conclusion
This study proposed production-forecasting model for the Montney shale reservoir. The study included procedures to select the significant reservoir and hydraulic fracture parameters and develop the production forecasting model. The following conclusions were made: (1) The relative significance of the reservoir and hydraulic fracture parameters on the production decline parameters was evaluated based on the p values derived from MANOVA (2) The reservoir thickness, porosity, permeability, initial reservoir pressure, main fracture half-length and conductivity, fracture initial water saturation, and the number of clusters per stage were selected as significant parameters in the fractured shale reservoir (4) The production-forecasting models showed a good match on the cumulative production despite the low accuracy on transition time forecasting: R 2 > 0:89 The production-forecasting model is very helpful for forecasting shale gas production without the numerical simulation if the reservoir properties are known Nomenclature P i : Reservoir pressure, psi H: Reservoir thickness, ft K m : Matrix permeability, md Φ m : Matrix porosity, fraction D nf : Natural fracture spacing, ft P L : Langmuir pressure, psi V L : Langmuir volume, scf/ton S wim : Matrix initial water saturation, fraction S wif : Fracture initial water saturation, fraction X mf : Main fracture half-length, ft C nf : Natural fracture conductivity, md-ft C mf : Main fracture conductivity, md-ft N cluster : Number of clusters per stage P hyp : Initial peak rate, scf/day b hyp : Hyperbolic decline exponent D hyp : Hyperbolic decline rate, /month D exp : Exponential decline rate, /month T transition : Transition time, month.

Data Availability
Data are available on request.

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