Development and Application of a Production Data Analysis Model for a Shale
Gas Production Well

This paper presents the development and application of a production data analysis software that can analyze and forecast the production performance and reservoir properties of shale gas wells. The theories used in the study were based on the analytical and empirical approaches. Its reliability has been confirmed through comparisons with a commercial software. Using transient data relating to multi-stage hydraulic fractured horizontal wells, it was confirmed that the accuracy of the modified hyperbolic method showed an error of approximately 4% compared to the actual estimated ultimate recovery (EUR). On the basis of the developed model, reliable productivity forecasts have been obtained by analyzing field production data relating to wells in Canada. The EUR was computed as 9.6 Bcf using the modified hyperbolic method. Employing the Pow Law Exponential method, the EUR would be 9.4 Bcf. The models developed in this study will allow in the future integration of new analytical and empirical theories in a relatively readily than commercial models.


Introduction
The shale reservoirs are known to have an enormous amount of natural gas. However, shale reservoirs have such low permeability that it releases gas very slowly, which is why shale is the last major source of natural gas to be developed. The most prolific shales are relatively flat, thick, and predictable, and the formations are so large that their wells will continue producing gas at a steady rate for decades [1,2].
In conventional oil and gas fields, transient flow generally appears during the short time period early in the production stages, while boundary dominated flow (BDF) occurs thereafter. During the transient flow period, reservoir properties are estimated by well testing analysis utilizing measured pressure. In addition, production forecasting is performed utilizing the decline curve analysis (DCA) proposed by Arps [3]. However, due to the low permeability of shale gas reservoirs, the transient flow period is relatively longer than that of conventional reservoirs. Therefore, Fetkovich [4] suggested that rate transient analysis (RTA) could be used to predict reservoir properties during the transient flow period.
Anderson et al. [5] suggested that to forecast productivity for shale gas wells, it is essential to consider the impact of fracture caused by the characteristics of matrix permeability and hydraulic fracturing. Production data analysis for shale gas wells uses production data and reservoir properties to calculate parameters for analyzing future productivity. And then it proposed that production decline rate and decline exponent would be calculated, and that derived the estimated ultimate recovery (EUR) [6,7]. Shale gas wells with multi-stage hydraulic fractured horizontal well (MHFW) are now the primary method for exploiting unconventional gas reservoirs. Analytical and empirical approaches for forecasting these wells were in their mature study. Analytical approaches, such as numerical and analytical simulation have been discussed in the literature [8][9][10].
Many studies [11][12][13][14][15][16][17] using the empirical DCA and flow characteristic have been conducted over the past years on shale gas wells with MHFW. Production data analysis that study a shale gas fields flow characteristics are actively being researched worldwide, and some commercial software has been developed to support this research. However, to estimate highly reliable results using commercial software, the analysis method should be understood beforehand. Moreover, additional time will be required for the latest analysis theories to be applied to commercial software. Furthermore, maintenance cost that are incurred each year must be considered. To facilitate efficient and effective shale gas production, methods to solve these problems will be required.
During this study, methods employed to analyze production data analysis of shale gas wells with MHFW were researched. In addition, an in-house software was developed applying the most current analysis theories. The developed software was verified utilizing a comparison with a similar commercial software package. This study provides research to enhance well productivity forecasting by utilizing the developed software for field production data analysis.
We have developed the spread-sheet based model which had applied the most current analysis theories. The developed model could reduce the additional time to require for the latest analysis theories to be applied to a commercial software, and maintenance costs that are incurred each year. Also, we have proposed the methodology of determining the parameters (b, D ∞ , etc.,) with the case studies in order to apply the theories (the Modified hyperbolic, the PLE method) to the field production data in Canada.

Production Forecasting Theory in Shale Reservoirs
Productivity is forecast by analyzing the decline trend with a DCA, employing production data acquired from conventional gas wells. However, the production data acquired from shale gas wells for one to two years displays only a transient flow pattern. Therefore, productivity forecast using a decline exponent produces an infinite acting scenario, and thus is not appropriate as an analysis method. Consequently, production forecasting should be performed by analyzing production performance while considering flow regime.
Given the production data of the transient flow region, which is demonstrated as a half slope on log-log plot (normalized rate vs. time), presented on the square root time (SRT) plot, the trend line's slope, m, is computed. Thus, a variable for productivity forecasting, t elf (end time of linear flow or transient flow), can be calculated as follows: where A is the area of stimulated reservoir volume (SRV) (multiplication of x and y in Fig. 1), h indicates the reservoir thickness, and φ refers to the reservoir porosity. μ g indicates gas viscosity, c t indicates total compressibility, m refers to the trend line's slope displayed in the SRT plot, p pi indicates initial pseudopressure, p pwf indicates pseudo-flowing bottom hole pressure, and T is the reservoir temperature.
Based on the formula above, the SRV area can be calculated if the original gas in place (OGIP) is computed through a linear regression (straight line extrapolation) on the BDF data utilizing a flowing material balance (FMB) plot. This yields the following equation: where A SRV defines the SRV area, OGIP refers to original gas in place, B gi refers to initial gas formation volume factor and S g refers to gas saturation.
The value of t elf cannot be obtained via linear regression from the FMB plot of transient flow period data, and thus, t elf should be obtained using SRV from the data acquired by microseismic analysis [18]. Fig. 2 presents the record of SRV on a microseismic signal. Given t elf , EUR can be calculated by applying the decline exponents of the hyperbolic decline curve equation before and after t elf using Eq. (3).
(3) Figure 1: Schematic diagram for multi-stage hydraulic fractured horizontal wells Figure 2: Example of microseismic data analysis [18] where q refers to the production rate, q i is the initial production rate, b is the decline exponent, and D i is the initial production decline rate. Nobakht et al. [9] analyzed by numerical simulation and found that b = 0.5 for homogeneous hydraulic fracturing reservoirs.
In general, the application of a modified hyperbolic DCA is limited owing to the uncertainty of microseismic data. Consequently, Ilk et al. [11] forecast that transient flow data significantly affects the decline rate; thus, they included a decline exponent to account for the variations in the hyperbolic decline curve at initial production periods of tight gas. Therefore, this study calculated the decline rate by analyzing the flow characteristics of shale and tight gas fields, and determined a future production rate utilizing Eq. (4).
where b q i refers to the production rate when t = 0. D ∞ is the decline rate at infinite time, b D i is the output computed by D i = D 1 /n and n refers to the time exponent. When the production rate was charted on a graph according to time change, and transient flow data was analyzed by setting a power law, specifically D ∞ = 0, the decline rate exhibited exponential decline. Ilk et al. [11] defined this as a power law exponential (PLE) DCA.

Development and Verification of a Production Data Analysis Model
The production data analysis model, which was designed with high standards for accessibility and convenience, was developed utilizing Microsoft Excel (Fig. 3). The model formula was converted into Visual Basic, and the arithmetic operations were designed to be automatically performed utilizing macro functions. The model was implemented to display the results graphically. The macro buttons for analysis were added to the ribbon menu at the top of the Excel window for a more intuitive use. Production data is stored in a database, and it can be imported into the Excel model when needed.
The developed model utilizes input data to analyze production performance and estimate productivity forecasting. When reservoir properties, oil well completion data, and fluid properties are entered in the input section, a macro is invoked to automatically calculate variables. The production performance analysis section generates the log-log plot, SRT plot and FMB plot, and computes the reservoir properties required for productivity forecasting (Fig. 4). Using the log-log plot, flow regime can be examined, and the linear flow parameters and pseudo-skin factor can be computed utilizing straight line fitting. The developed model can compute the SRV area and fracture half-length, if OGIP is estimated from the FMB The productivity forecasting section is divided into the Arps DCA, modified hyperbolic method, and PLE method. Furthermore, it was designed to allow an analysis between production data and production performance to be organically feasible. To verify the validity of the developed model, we compared our production data analysis with the results Fekete Harmony commercial software. A mutual comparison of the results verified the validity of the model and field applicability.
By analyzing the straight-line trend of either model, the two models' results were observed to be essentially equal (Fig. 5). Graphical analysis showed that the data deviated from the straight line in some sections; we believe these deviations are due to the difference in the method employed to calculate the pseudo-pressure function and the filtering function of the commercial model. Slight differences in the linear flow parameters computed from the SRT analysis using the developed model and the Fekete Harmony model were observed as 41,900 md 1/2 ft 2 and 37,071 md 1/2 ft 2 , respectively. These differences are believed to be the result of a difference in data selection required to perform regression analysis from the respective models. Furthermore, this produced a slight difference in the y-intercept value, b. However, an FMB plot analysis of the two models showed that their trendlines virtually matched (Fig. 6). Similar to the SRT plot analysis, some differences were observed, but this determined due to differences between the pseudo-pressure calculation methods used between the models. As a result of observing the productivity forecasting, the EUR was calculated as 3,652 MMscf for the developed model and 3,667 MMscf for commercial software, with only 0.4% of the error.

Synthetic Data Analysis
Although analysis methods that consider the flow characteristics of shale gas fields were presented, the EUR might differ according to analysis methods that only employ transient flow data. To present a valid analysis method, this study generated production data by building an MHFW model, reflecting the characteristics of a shale gas reservoir (Fig. 7). Tab. 1 lists the main properties of the reservoir model used in this study. Fifty years of production data were generated from January 2013 to 2062 utilizing Schlumberger's Eclipse simulator, until the production rate reached the economic limit of 80 MMscf/d (Fig. 8). Furthermore, cumulative production of gas was computed as approximately 9 Bcf. This study selected production data beginning 1.6 years after production when transient flow occurred. Different productivity forecasting techniques were used for each of the three cases. First, case A is productivity was forecast by applying Arps DCA method. Second, case B calculated the t elf through SRV and applied different productivity forecasting techniques depending on the flow regime. Finally, in case C maximum cumulative production was determined utilizing a production decline rate and time exponent, in line with the RTA methods described above, by applying the PLE method. A valid analysis method was presented by forecasting productivity according to each analysis method, as well as the EUR comparison. In Case of A, a decline exponent of Arps DCA method, which is generally used in conventional oil and gas production processes, was estimated using 0.5.  EUR was calculated as 5 Bcf with a relative error of approximately 44% from the actual EUR (Fig. 9). Therefore, the first analysis method was confirmed to be inappropriate. The model failed due to the long transient flow period of a shale gas reservoir, which spans several years. To correct for this, the productivity during this period is defined as reservoir BDF and a decline exponent was applied.
In Case of B, t elf was calculated as 6.63 years utilizing the trend line slope, which was derived from plotting production data on the SRT plot and the SRV area calculated from the reservoir system (Fig. 10). Decline exponents of 1.5 and 0.5 were applied to the transient flow period and BDF period, respectively. In this model, EUR was 9.4 Bcf, which was approximate to the actual EUR with a relative error rate of roughly 4%. Based on these observations, it was determined that this method had been valid for productivity forecasting. Because this analysis employed simulated data, which can be significantly more reliable than field data, it was easier to obtain more precise estimates of SRV and t elf . However, when utilizing field data, microseismic data should be used to increase the reliability of SRV calculations, and subsequently, EUR estimates.
In Case of C (Fig. 11), The initial production decline rate calculated from the production data and time exponent was 1.91. As the production period becomes infinitely long, the production decline rate approaches 0; 1.0 Â 10 −6 Production period with economic rate (year) 50 Figure 8: Production performance of simulated reservoir using synthetic data thus, D ∞ = 0. A time exponent of 0.11 was applied using a sensitivity analysis of the production decline trend. As a result, maximum cumulative gas production was 9.4 Bcf with a relative error of approximately 4%; therefore, the results were not significantly different from the actual EUR.  The PLE method does not define t elf , but instead forecasts productivity applying only a time exponent; therefore, this method is judged to be robust in cases where, either microseismic data does not exist, or its analysis is difficult.

Field Data Analysis
In this study, productivity forecasting methods were evaluated utilizing production data from a shale gas field (Horn River basins) located in British Columbia (Fig. 12). A productivity forecasting model was investigated using the modified hyperbolic method and PLE method. Analysis was based on a production period of one year and three months, which began after the start of production in production Well A. Tab. 2 shows the reservoir and fluid properties for field data. Based on a summary of deviation survey for the field data of shale gas well A, B (Tabs. 3 and 4). The horizontal well length was 15,446 ft. Additionally, 17-stage hydraulic fracturing was applied, and production was performed from June 2012 to August 2013 (Fig. 13).
In case of shale gas well A, The SRT plot's trendline slope for pseudo-pressure was estimated, as shown in Fig. 14. Subsequently, t elf was estimated utilizing the SRT and SRV slope obtained from the microseismic data. t elf was forecast as June 2015, which is three years after starting production. Based on the results, productivity was forecast by applying the decline exponent to the hyperbolic decline curve equation (Fig. 15). The EUR was computed as 9.6 Bcf. In addition, the initial decline rate computed from the  production data and time exponent was 2.21 when applying the PLE method. Because the decline rate approaches 0 as the production period extends indefinitely, D ∞ = 0. The sensitivity analysis results showed that the time exponent was 0.13. the maximum cumulative production estimated by this was 9.4 Bcf (Fig. 16).   The reservoir property of shale gas well B, the initial reservoir pressure was 3,190 psia, reservoir temperature and thickness were 149°F and 147 ft, respectively, and the horizontal well length was 14,388 ft. In addition, 17-stage hydraulic fracturing was applied. This production well was utilized from June 2012 to August 2013 (Fig. 17). In the case of production well B, regression analysis produced an irregular scatterplot of normalized pressure, as was observed on the SRT plot. To correct this problem, the time term of the SRT plot was comparatively analyzed by utilizing a superposition time function plot (Fig. 18).
Production forecasting results are shown in Fig. 19. Initially, t elf was estimated from the SRV using the microseismic data, and the trend line slope was computed from the SRT plot. As a result, t elf was forecast to be March 2019, which was six years and nine months from the start of production. Accordingly, the EUR was computed as 7.7 Bcf based on the forecasted productivity results. For the subsequent analysis method, t elf was calculated using a regression analysis based on the superposition time SRT graph and was forecast to be March 2015. Based on this, the decline exponent according to flow regime was applied to the hyperbolic decline curve equation, future productivity was forecast, and the EUR was calculated as 5.5 Bcf.
To establish the most appropriate analysis method among the above results, they were each compared with t elf of well A. Because well B was operated in the same layer as well A, reservoir properties may be similar. When t elf was compared between estimation methods, the application of superposition time  produced a similar result of June 2015. This date is when the transient flow ends, which was calculated by utilizing the production data acquired from well A. In the interim, the initial production decline rate was computed as 0.31 and time exponent as 0.29 by applying the PLE method. The maximum cumulative gas production calculated by the analysis method above was computed as 5.2 Bcf (Fig. 20).

Conclusion
This study developed a production data analysis software that can analyze and forecast the production behavior of shale gas wells and confirmed the software's validity by comparing its analysis results with a commercially available program. Utilizing this model, this study provided research to support well productivity forecasting by analyzing field production data from wells in British Columbia. The developed model can improve the results and rapidly reflect new analytical and empirical methods compared to commercial models.
1. The software used to analyze shale gas well production data was developed utilizing MS Excel; the software supports well production performance analysis and productivity forecasting. The developed model was calculated to be robust, as it produced results similar to that of a commercial program and were statistically significant within an error range of 5%. 2. In this study, productivity forecasting methods were evaluated using initial production data produced from an MHFW model of a simulated shale gas reservoir. The modified hyperbolic method was validated by an EUR result with an error rate of approximately 4%, compared with the actual EUR. If t elf is computed utilizing production performance analysis and microseismic data, the reliability of the results can be increased. If the t elf calculation becomes difficult, the PLE method becomes more appropriate. 3. This study analyzed future productivity by applying the modified hyperbolic method on one year and three months of production data from shale gas production well A. t elf was forecast to be June 2015, three years after beginning of production, and the EUR was computed as 9.6 Bcf. Based on the results of productivity analysis employing the PLE method, the EUR was 9.4 Bcf. 4. Assuming production well B's reservoir properties were similar to well A, we applied a superposition time plot to estimate t elf . The difference in results between models was statistically insignificant; thus, this method was judged to be valid. Our model produced a t elf of March 2015, and the EUR was computed as 5.5 Bcf.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.