The Evaluation and Sensitivity of Decline Curve Modelling

The development of prediction tools for production performance and the lifespan of shale gas reservoirs has been a focus for petroleum engineers. Several decline curve models have been developed and compared with data from shale gas production. To accurately forecast the estimated ultimate recovery for shale gas reservoirs, consistent and accurate decline curve modelling is required. In this paper, the current decline curve models are evaluated using the goodness of fit as a measure of accuracy with field data. The evaluation found that there are advantages in using the current DCA models; however, they also have limitations associated with them that have to be addressed. Based on the accuracy assessment conducted on the different models, it appears that the Stretched Exponential Decline Model (SEDM) and Logistic Growth Model (LGM), followed by the Extended Exponential Decline Model (EEDM), the Power Law Exponential Model (PLE), the Doung’s Model, and lastly, the Arps Hyperbolic Decline Model, provide the best fit with production data.


Introduction
In recent years, shale gas reservoirs (SGR) or unconventional reservoirs have steadily become the main bases of natural gas production around the world [1]. Wang [2] notes that shales and sediments are the richest sedimentary rocks in the Earth's crust and, according to recent activities, shale gas will constitute the largest component in gas production globally, as conventional reservoirs continually decrease. It is further mentioned by Wang [2] that SGR, unlike conventional reservoirs, tend to be more costly to develop and require special tools to enable the gas to be produced at a cost-effective rate due to their extremely low matrix permeability and porosity [3]. Accordingly, the modelling of shale gas production and its decline is essential to predict how fast the gas can be produced and turned into revenue from each well, as well as the feasibility of producing natural gas from operated shale plays from a cost perspective [2].
Currently, the oldest and most commonly used tool for the modelling of shale gas production is the rate versus time decline curve estimation due to its ease. Current efforts in decline curve analysis (DCA) have been concentrating on a computer statistical approach, the basic objective being to arrive at a distinctive "unbiased" interpretation [2]. In recent years, several DCA models have been suggested and compared with previous shale gas production figures, prior to being used on more reservoirs [4]. This paper focusses on the evaluation and sensitivity of the current DCA models and proposes a new hybrid model to be investigated in SGR decline analysis. The main ideas are (a) to characterise and evaluate the current decline curve models used to explain shale gas reservoir forecasting and (b) use the goodness-of-fit regression test to assess the sensitivity of the decline curve models in (a).

•
Fracture/Early Linear Flow: A transient flow regime that occurs when the production flow is linear to the single fractures. This flow regime governs the known life of most shale wells. A negative half slope on a log-log plot of rate versus time can be used to differentiate this linear flow. • Fracture Boundary Flow: Follows after a certain period of production when an interference occurs i.e., from linear to simulated reservoir volume (SRV). Many of the existing horizontal shale wells have not experienced this regime, but some of the newer wells with huge fracture treatments have been observing this regime early. This can be observed on a log-log plot by deviation from a -1/2 slope line on a log-log plot of rate versus time. • Matrix Linear Flow: When production from the matrix, beyond the SRV, starts to govern the production, a linear type flow will be seen. This regime is most likely will not be observed in the economic life of the well. Comparable to fracture linear flow, this regime can be observed using a negative half slope line on a log-log plot of rate versus time. • Matrix Boundary Flow: After the outer matrix transient has reached the drainage boundaries of the well, a deviation from the negative half slope, corresponding to matrix linear flow, will be observed. This deviation is equivalent to matrix boundary flow. Similar to the matrix, linear flow will most likely not be observed.

Overview of Decline Curve Models
Consistently forecasting the long-term production performance of shale (unconventional) reservoirs has been a challenge [11]. The petroleum industry requires simple, useful, and speedy means of predicting production and assessing reserves; hence, DCA has been an attractive alternative in contrast to other methods [11]. Due to the relative ease of DCA, it is considered the most used Energies 2020, 13, 2765 3 of 16 method in the industry [11]. The current DCA models will be evaluated based on their characteristics, strengths, weaknesses, and sensitivity to production data.

Arps Decline Curve and the Modified Hyperbolic Decline Model (MHD)
Arps decline curve analysis is the most commonly used method of estimating ultimate recoverable reserves and future performance [12]. Paryani et al. [13] reasons this to reliable history match (even with b > 1) and its simplicity. The model process is based on vital assumptions: that past operating conditions will remain unaffected, a well is produced at or near capacity, and the well's drainage remains constant and is produced at a constant bottom-hole pressure [14]. Notably, the Arps model is only applicable in pseudo-steady flows when the flow regime transfers from linear flows to boundary-dominated flows (BDF) [15]. This indicates the Arps equations are not applicable to the production forecasting of the entire decline process of horizontal wells in low-permeability reservoirs [16]. The Arps decline curve analysis can be summarised into three types: exponential Equation (1), hyperbolic Equation (2), and harmonic Equation (3) [17,18].
where q is the flow rate in STB/day or Mscf/day, q i is the initial flow rate in STB/day or Mscf/day, D is the decline constant while D i is the initial decline constant, which are both measured in days − 1, and b is the decline exponent. The most commonly employed hyperbolic form of Arps decline Equation (2) is used for shale reservoirs. The hyperbolic decline equation is suitable to use due to the "best fit" that it provides for the long transient linear-flow regime observed in shale gas wells with b values greater than unity [18]. The model results in post-production overestimation due to the decrease in the decline rate with production time. Due to the overestimation, Robertson et al. [19] suggested a revised version of the hyperbolic decline model for shale gas production decline. The equation is given as: where q is the production rate in m 3 /d or STB/day, D lim is the decline rate in d −1 , and n is the time exponent. They suggested that the hyperbolic decline model sometimes yields unrealistically high reserve estimates. They made an assumption that the rate of decline starts at 30% of flow and usually declines in a hyperbolic way [19]. This modified model considers when the hyperbolic decline in the early life of a well transfers to exponential decline in the late life [19]. The switching process can be determined by applying computer programs. The switching point is when the decline rate is smaller than a certain limit (usually 5%) [19]. The MHD model addresses the overestimation limitation of EUR; however, it is still unable to determine D lim for production data [15].
To test the behavior of the Arps hyperbolic model and the modified version shown in Figure 1, a semi log plot (log q versus t) illustrates the sensitivity of the models to various estimated field data. The R 2 values denote the goodness of fit or the degree of linear correlation, which is a measure of the level of association of a group of actual observations to the model's forecasts [20]. As observed from the regression lines for the various data, the resulting fit appears to capture the trend in the data well. Arps fits Data 1 and 2 fairly, similarly for the MHD. However, the methods matches the other cases poorly, because it cannot model multiple flow regimes. In the instance of the MHD model, there is a shift in the curves downward, which results in a change in the R 2 value. Upon closer inspection of the Arps fits Data 1 and 2 fairly, similarly for the MHD. However, the methods matches the other cases poorly, because it cannot model multiple flow regimes. In the instance of the MHD model, there is a shift in the curves downward, which results in a change in the R 2 value. Upon closer inspection of the EUR values for both models, which are shown in Table 1, it is evident that the MHD model corrects for the overestimation of the Arps model.

Power Law Exponential Model (PLE)
Ilk et al. [22] presented the PLE, which is an extension of the exponential Arps formula for the decline degree in shale reservoirs. This model was developed precisely for SGR and approximates the rate of decline with a power law decline. The PLE model matches production data in both the transient and boundary-dominated regions without being hypersensitive to remaining reserve estimates [23]. Seshadri and Mattar [24] presented that the PLE model can model transient radial and linear flows, while Kanfar and Wattenbarger [25] proved that the model is reliable for linear flow, bilinear flow followed by linear flow, and linear flow followed by BDF, or bilinear flow followed by linear flow and finished with BDF flow. Vanorsdale [26] deduced that when the flow regime changes throughout the initial 10 years of the well, the PLE model will yield a very optimistic recovery. The model characterizes the decline rate by infinite time, D∞ which is defined as a "loss ratio" (which is assumed to be constant from Arp) [16]. The production rate is derived as follows: where / is the slope, D∞ is the decline rate over a long-term period, and is the time exponent.
By substituting the above equations, the production rate is obtained:

Power Law Exponential Model (PLE)
Ilk et al. [22] presented the PLE, which is an extension of the exponential Arps formula for the decline degree in shale reservoirs. This model was developed precisely for SGR and approximates the rate of decline with a power law decline. The PLE model matches production data in both the transient and boundary-dominated regions without being hypersensitive to remaining reserve estimates [23]. Seshadri and Mattar [24] presented that the PLE model can model transient radial and linear flows, while Kanfar and Wattenbarger [25] proved that the model is reliable for linear flow, bilinear flow followed by linear flow, and linear flow followed by BDF, or bilinear flow followed by linear flow and finished with BDF flow. Vanorsdale [26] deduced that when the flow regime changes throughout the initial 10 years of the well, the PLE model will yield a very optimistic recovery. The model characterizes the decline rate by infinite time, D∞ which is defined as a "loss ratio" (which is assumed to be constant from Arp) [16]. The production rate is derived as follows: where dq/dt is the slope, D ∞ is the decline rate over a long-term period, andn is the time exponent. By substituting the above equations, the production rate is obtained: In this model, there are four unknown variables:q i ,D i, D ∞ andn, which result in several degrees of freedom and may be clumsy to use or solve [27]. According to Johnson et al. [28], the D∞ parameter is difficult to determine. However, there are advantages to this model in that the extra variables permit for both transient and boundary flow, and the equation for production rate seems comparable to the Arps exponential equation [13]. With the PLE model ( Figure 2), which uses a log-log plot (log q versus log t) to test the sensitivity of the data, the resulting fit appears to capture the trend in the data better compared to the Arps Hyperbolic Model. This model fits Data 1, 2, 3, and 4 fairly accurately. This can be attributed to the PLE model, matching production data in both transient and boundary-dominated regions.
In this model, there are four unknown variables: , , , which result in several degrees of freedom and may be clumsy to use or solve [27]. According to Johnson et al. [28], the D∞ parameter is difficult to determine. However, there are advantages to this model in that the extra variables permit for both transient and boundary flow, and the equation for production rate seems comparable to the Arps exponential equation [13]. With the PLE model ( Figure 2), which uses a log-log plot (log q versus log t) to test the sensitivity of the data, the resulting fit appears to capture the trend in the data better compared to the Arps Hyperbolic Model. This model fits Data 1, 2, 3, and 4 fairly accurately. This can be attributed to the PLE model, matching production data in both transient and boundary-dominated regions.

Stretched Exponential Decline Model (SEDM)
Valkó, Valkó, and Lee [29,30] applied the SEDM in shale wells, which is an empirical method different from Arps equations, as it describes the decline trend of production data obtained from unconventional reservoirs. It was developed to fit transient flow regimes [10,25]. The significant advantages of the model are the bounded nature of estimated ultimate recovery (EUR) without limits on time or rate, and the straight-line behavior of a recovery potential expression [30]. The model differs from other models, since it does have a basis in physics and is directed by a major differential equation [14]. It is used to model aftershock decay rates [31]. The production rate declines with time according to the following equations: This method defines a characteristic number of periods, τ, and a dimensionless exponent, n, of the ratio of time, t. It also uses observed cumulative production along with theoretical cumulative production derived from the integral of the rate-time equation to estimate remaining technically recoverable volumes. Equation

Stretched Exponential Decline Model (SEDM)
Valkó, Valkó, and Lee [29,30] applied the SEDM in shale wells, which is an empirical method different from Arps equations, as it describes the decline trend of production data obtained from unconventional reservoirs. It was developed to fit transient flow regimes [10,25]. The significant advantages of the model are the bounded nature of estimated ultimate recovery (EUR) without limits on time or rate, and the straight-line behavior of a recovery potential expression [30]. The model differs from other models, since it does have a basis in physics and is directed by a major differential equation [14]. It is used to model aftershock decay rates [31]. The production rate declines with time according to the following equations: This method defines a characteristic number of periods, τ, and a dimensionless exponent, n, of the ratio of time, t. It also uses observed cumulative production along with theoretical cumulative production derived from the integral of the rate-time equation to estimate remaining technically recoverable volumes. Equation (10) appears similar to the PLE model; however, it differs, as it does not rely on a single interpretation of parameters. Instead, it uses two-parameter gamma functions [29]. In addition, there is no single τ and n parameters, but instead, a sum of multiple exponential declines, which follows the fat tail distribution [30]. Stretched Exponential Decline Model (SEDM) requires an iterative process to determine the value of the parameter, n. The model can only estimate the recoverable volumes with an abandonment rate of zero as opposed to commercial volumes with economic cut-off rates and has not been widely used [32]. However, Can et al. [32] showed that in tight formations where transient flow period is extremely long, SEDM has been successful in modeling the rate-time behavior and provides more realistic reserve estimates compared to Arps decline relations.
Testing the behavior of the SEDM, Figure 3, which is a plot of production rate versus the cumulative production (q versus Q) to test the sensitivity of the data, the resulting fit appears to capture the trend in the data poorly. The SEDM method fits all cases inaccurately (lower R 2 values). This is due to the SEDM model's transient flow rather than boundary-dominated flow and requirement for a sufficiently long production time (usually >36 months) to accurately estimate the parameters τ and n [33].
which follows the fat tail distribution [30]. Stretched Exponential Decline Model (SEDM) requires an iterative process to determine the value of the parameter, n. The model can only estimate the recoverable volumes with an abandonment rate of zero as opposed to commercial volumes with economic cut-off rates and has not been widely used [32]. However, Can et al. [32] showed that in tight formations where transient flow period is extremely long, SEDM has been successful in modeling the rate-time behavior and provides more realistic reserve estimates compared to Arps decline relations.
Testing the behavior of the SEDM, Figure 3, which is a plot of production rate versus the cumulative production (q versus Q) to test the sensitivity of the data, the resulting fit appears to capture the trend in the data poorly. The SEDM method fits all cases inaccurately (lower R 2 values). This is due to the SEDM model's transient flow rather than boundary-dominated flow and requirement for a sufficiently long production time (usually >36 months) to accurately estimate the parameters τ and n [33].

The Extended Exponential Model (EEDM)
Zang et al. [11] presented a renewed experimental method, the EEDM, as a simple formula to forecast shale oil and gas well performance. They proposed a mechanism of "growing drainage volume" to conceptualize and model the performance of shale wells. This model combines the exponential decline equation proposed by Fetkovich et al. [34] Equation (13) with the derived empirical Equation (14). The EEDM includes both transient and BDF flow in a single equation, and it can match the historical data with a smooth curve throughout the transition period from transient to BDF flow regimes. Furthermore, the model is simple and can easily be applied [11]. It is also able to project the future production by fitting all of the historical production data from the beginning of the production decline.
Paryani et al. [13] stated that the model contains two decline constants and a decline exponent; particularly noteworthy, the production data fits using a smooth curve through the whole flow systems [16]. The advantage of the model is that both early and late production profiles can be captured once have been calibrated using the production data [11]. However, as parameter has an incomplete influence on the curve fitting, it is therefore fixed.

The Extended Exponential Model (EEDM)
Zang et al. [11] presented a renewed experimental method, the EEDM, as a simple formula to forecast shale oil and gas well performance. They proposed a mechanism of "growing drainage volume" to conceptualize and model the performance of shale wells. This model combines the exponential decline equation proposed by Fetkovich et al. [34] Equation (13) with the derived empirical Equation (14). The EEDM includes both transient and BDF flow in a single equation, and it can match the historical data with a smooth curve throughout the transition period from transient to BDF flow regimes. Furthermore, the model is simple and can easily be applied [11]. It is also able to project the future production by fitting all of the historical production data from the beginning of the production decline.
Paryani et al. [13] stated that the model contains two decline constants and a decline exponent; particularly noteworthy, the production data fits using a smooth curve through the whole flow systems [16]. The advantage of the model is that both early and late production profiles can be captured once β e and β l have been calibrated using the production data [11]. However, as parameter β l has an incomplete influence on the curve fitting, it is therefore fixed.
Energies 2020, 13, 2765 where a is the nominal decline rate, β l is the late-life period constant, and β e is the early period constant. Combining Equations (13) and (14) and taking the logarithm of each side, the equation below (the exponential decline equation) is obtained.
where q o is the initial production rate in m 3 /s. Using the EEDM (Figure 4), which is a plot of −ln q qo /t versus t to test the sensitivity of the data, the resulting fit appears to also capture the trend in the data poorly. The method fits all cases inaccurately (lower R 2 values). This type of method is best to forecast short-term trends in the absence of recurring variations. Hence, the EEDM would only be accurate when a realistic amount of stability between the past and future is assumed.
where a is the nominal decline rate, βl is the late-life period constant, and βe is the early period constant. Combining Equations (13) and (14) and taking the logarithm of each side, the equation below (the exponential decline equation) is obtained.
where qo is the initial production rate in m 3 /s. Using the EEDM (Figure 4), which is a plot of − / versus t to test the sensitivity of the data, the resulting fit appears to also capture the trend in the data poorly. The method fits all cases inaccurately (lower R 2 values). This type of method is best to forecast short-term trends in the absence of recurring variations. Hence, the EEDM would only be accurate when a realistic amount of stability between the past and future is assumed.

Doung's Decline Model
Doung [35] presented an unconventional rate decline method to evaluate the performance of shale gas wells that does not depend on the fracture types. The model assumes linear or near-linear flow, as indicated by a log-log plot of rate over cumulative production versus time, which yielded a straight-line tendency [36]. The rate is calculated in the model using the following equation [27]: where t (a,m) is the time constant in 1/s, and q∞ is the production rate at infinite time in m 3 /s. The cumulative production and time constant is calculated as: where Gp is the cumulative gas production in Bcf and m is the slope. Paryani et al. [13] indicated the key restrictions of the model: if the well is closed for extended periods, a proper rate initialization against pressure is required to obtain precise values of parameters

Doung's Decline Model
Doung [35] presented an unconventional rate decline method to evaluate the performance of shale gas wells that does not depend on the fracture types. The model assumes linear or near-linear flow, as indicated by a log-log plot of rate over cumulative production versus time, which yielded a straight-line tendency [36]. The rate is calculated in the model using the following equation [27]: where t (a,m) is the time constant in 1/s, and q ∞ is the production rate at infinite time in m 3 /s. The cumulative production and time constant is calculated as: where G p is the cumulative gas production in Bcf and m is the slope. Paryani et al. [13] indicated the key restrictions of the model: if the well is closed for extended periods, a proper rate initialization against pressure is required to obtain precise values of parameters a and m and, secondly, that in the event of water breakthrough, there is a sudden decrease in the Energies 2020, 13, 2765 8 of 16 decline rate, and this causes an increase in the values of the a and m parameters. Vanorsdale [26], similar to in the case of the PLE model, also indicated that the Doung model will yield a very optimistic recovery when the flow regime changes throughout the initial 10 years. He went on to indicate that the model may provide conservative recovery estimate in vertical, non-hydraulic fractured classical shale wells [26]. However, Lee et al. [36] has indicated that the Doung's model appears to fit field data from various shale plays quite well and provides an effective alternative to Arps hyperbolic model. With the Duong's model (Figure 5), which uses a log-log linear plot (log q versus log t) to test the sensitivity of the data, the resulting fit appears to capture the trend in the data well. The method fits Data 1, 2, and 4 fairly accurately. For Case 3, the method fits the data poorly with a lower R 2 value of 0.8371. The model probably provides a good fit because it was specifically developed for unconventional reservoirs with very low permeability. a and m and, secondly, that in the event of water breakthrough, there is a sudden decrease in the decline rate, and this causes an increase in the values of the a and m parameters. Vanorsdale [26], similar to in the case of the PLE model, also indicated that the Doung model will yield a very optimistic recovery when the flow regime changes throughout the initial 10 years. He went on to indicate that the model may provide conservative recovery estimate in vertical, non-hydraulic fractured classical shale wells [26]. However, Lee et al. [36] has indicated that the Doung's model appears to fit field data from various shale plays quite well and provides an effective alternative to Arps hyperbolic model. With the Duong's model ( Figure 5), which uses a log-log linear plot (log q versus log t) to test the sensitivity of the data, the resulting fit appears to capture the trend in the data well. The method fits Data 1, 2, and 4 fairly accurately. For Case 3, the method fits the data poorly with a lower R 2 value of 0.8371. The model probably provides a good fit because it was specifically developed for unconventional reservoirs with very low permeability.

Logistic Growth Model (LGM)
Logistic Growth Models developed belong to a group of mathematical models used to forecast growth in numerous applications [36] and were previously used to model population growth [37,38]. It was developed to forecast reservoirs with extremely low permeability [27]. LGM is very flexible and confident in modelling long transient boundary-dominated performances of unconventional reservoirs [16]. The model incorporates known physical volumetric quantities of oil and gas into the forecast to constrain the reserve estimate to a reasonable quantity.
LGM is capable of trending existing production data and providing reasonable forecasts of future production. The logistic growth model does not extrapolate to non-physical values [38]. Tsoularis and Wallace [39] discussed a development in this regard by Verhulst [40], who considered that for the population model, a steady population would consequently possess a saturation level characteristic, typically called the carrying capacity, K, which forms a numerical upper bound on the growth size. In order to include this limiting characteristic, they introduced the logistic growth equation as an extension to the exponential model [39]. Zhang et al. [1] adopted this model for SGR with very low permeability and developed the LGM as an empirical method to forecast the gas production. The LGM can be represented as follows:

Logistic Growth Model (LGM)
Logistic Growth Models developed belong to a group of mathematical models used to forecast growth in numerous applications [36] and were previously used to model population growth [37,38]. It was developed to forecast reservoirs with extremely low permeability [27]. LGM is very flexible and confident in modelling long transient boundary-dominated performances of unconventional reservoirs [16]. The model incorporates known physical volumetric quantities of oil and gas into the forecast to constrain the reserve estimate to a reasonable quantity.
LGM is capable of trending existing production data and providing reasonable forecasts of future production. The logistic growth model does not extrapolate to non-physical values [38]. Tsoularis and Wallace [39] discussed a development in this regard by Verhulst [40], who considered that for the population model, a steady population would consequently possess a saturation level characteristic, typically called the carrying capacity, K, which forms a numerical upper bound on the growth size. In order to include this limiting characteristic, they introduced the logistic growth equation as an extension to the exponential model [39]. Zhang et al. [1] adopted this model for SGR with very low permeability and developed the LGM as an empirical method to forecast the gas production. The LGM can be represented as follows: where K is the carrying capacity. The main benefit of the LGM is that the reserve estimate is inhibited by the parameter K as well as the production rate, which terminates at infinite time [1]. The main assumption in this model is that the whole reservoir can be drained by a single well over a suitably long period and requires the approximation of at least two parameters, or parameters as per the available well information [4,7]. Figure 6, a plot of production rate versus time (q versus t), illustrates the sensitivity of the model to various estimated field data. As observed from the regression lines for the various data, the resulting fit appears to capture the trend in the data well. The LGM fits Data 1 and 2 fairly. However, the method matches the other cases poorly, as indicated by the lower R 2 values. This could be attributed to the data size, which is too small to yield an accurate fit, since the underlying principle of this model is population growth, which stipulates that growth is only possible up to a certain size.
Energies 2020, 13, x FOR PEER REVIEW 9 of 17 where K is the carrying capacity. The main benefit of the LGM is that the reserve estimate is inhibited by the parameter K as well as the production rate, which terminates at infinite time [1]. The main assumption in this model is that the whole reservoir can be drained by a single well over a suitably long period and requires the approximation of at least two parameters, or parameters as per the available well information [4,7]. Figure 6, a plot of production rate versus time (q versus t), illustrates the sensitivity of the model to various estimated field data. As observed from the regression lines for the various data, the resulting fit appears to capture the trend in the data well. The LGM fits Data 1 and 2 fairly. However, the method matches the other cases poorly, as indicated by the lower R 2 values. This could be attributed to the data size, which is too small to yield an accurate fit, since the underlying principle of this model is population growth, which stipulates that growth is only possible up to a certain size.

Autoregressive Intergrated Moving Average (ARIMA) and Neutral Network Models (NNM) (Hybrid Model)
The accuracy of time series forecasting is challenging for scientists [41]. Time series data often comprise linear as well as non-linear components [42]. In some cases, linear-based approaches might be more suitable than non-linear ones due to the data characteristics. The hybrid method is a combination of ARIMA and the neural network method. According to Faruk [42], hybrid methods have a higher degree of accuracy than neural networks. ARIMA can recognize time-series patterns well but not non-linear data patterns. On the other hand, neural networks only handle non-linear data. Therefore, hybrid models combine the advantages of ARIMA with respect to linear modelling and neural networks in terms of non-linear edge modelling [43]. Notwithstanding, in some circumstances, the single model approach can outperform hybrid models [41] Mathematically, time-series data can be expressed as a combination of linear and non-linear components [44]: where Yt shows the time-series data, Lt indicates the linear components, and the non-linear components are represented by Nt.
Mathematically, the neural network model for residual of n input nodes can be expressed as the following:

Autoregressive Intergrated Moving Average (ARIMA) and Neutral Network Models (NNM) (Hybrid Model)
The accuracy of time series forecasting is challenging for scientists [41]. Time series data often comprise linear as well as non-linear components [42]. In some cases, linear-based approaches might be more suitable than non-linear ones due to the data characteristics. The hybrid method is a combination of ARIMA and the neural network method. According to Faruk [42], hybrid methods have a higher degree of accuracy than neural networks. ARIMA can recognize time-series patterns well but not non-linear data patterns. On the other hand, neural networks only handle non-linear data. Therefore, hybrid models combine the advantages of ARIMA with respect to linear modelling and neural networks in terms of non-linear edge modelling [43]. Notwithstanding, in some circumstances, the single model approach can outperform hybrid models [41].
Mathematically, time-series data can be expressed as a combination of linear and non-linear components [44]: where Y t shows the time-series data, L t indicates the linear components, and the non-linear components are represented by N t . Mathematically, the neural network model for residual of n input nodes can be expressed as the following: e t = f (e t−1 + e t−2 , . . . , e t−n ) Energies 2020, 13, 2765 10 of 16 where f is a non-linear function that is specified by the neural network. With regard to the results of the prediction error of N t , the combination forecast using the hybrid method can be expressed as: There has been limited work conducted using this model for shale gas reservoirs. Hence, the next step would be to investigate this model for shale gas reservoirs.
To summarize all eight DCA models for an easy reference of readers, Table 2 lists the name of each model, its DCA equation, the characteristic, strength, weakness, and lastly the related references. linear to BDF flow reliable and simple to use post-production overestimation [12][13][14][15][16][17][18] 2 Modified Hyperbolic Curve linear and non-linear high degree of accuracy approach can be found to not be fit all types of data [40][41][42][43]

Accuracy of Current Decline Curve Models with Field Data
Yuhu et al. [15] discussed comparisons of EURs with five types of decline models from single-well production data. They explained that according to the prediction results, the highest predicted EUR was gained by the hyperbolic decline model, followed by the Modified Hyperbolic Model (MHD), Doung's Model, PLE and, lastly, the EDM. Hu et al. [27] conferred production data for wells with a production time greater than 10 years, for which the PLE decline model was recommended for multiple flows. It was also pointed out that the hyperbolic decline model predicted higher estimates of reserves than the PLE decline model. Another study that they reviewed recommended the MHD rather than the PLE decline model, which in their view was complicated.
It is noted that the differences in EURs with different decline models decrease with an increase in production time [45]. On the other hand, prediction consistency increases with an increase in production time. Based on this distinctive production data, the order of predicted EURs from high to low were through the hyperbolic decline model, the MHD, the PLE decline model, and the EDM respectively [45]. The predicted EURs decreased with an increase of production time for the hyperbolic decline and the modified hyperbolic decline model. The predicted EURs increase with an increase of production time for the PLE decline and the EDM model [45]. Currently, the applicability of these different decline models is uncertain. The general trend found in their paper was that the hyperbolic decline model overestimates the production and that the other decline models will still have to be investigated for reliability and accuracy [45].
In their study, Guo et al. [46] investigated shale gas wells in the Barnett shale play, where they found that from the results of goodness of fit, the hyperbolic curve fits well for both the aggregate and individual shale gas wells. On the other hand, Kenomore et al. [47] in their production decline study of the Barnett shale found that either the Arps hyperbolic or Doung's model can be used only if the historical data exceeds 10 months. They used root mean square error (RMSE) analysis and the results indicated that the Arps hyperbolic model showed better forecasting compared to the Doung's model for the top three longest production histories. Zhang et al. [1] concurred with the findings of the Doung's model in their paper, noting that it is more accurate for linear flow and bilinear flow; however, if the production history is shorter than 18 months, this model provides unreliable results for EUR. In most circumstances, the Doung's model overestimates the total EUR. Harris [48], in his research study of the Elm Coulee field production data, found that the Duong method will produce the most optimistic forecasts followed by the Arps model with 5% minimum decline, and then the SEPD model. Shah [49] in his research developed new methods of combining the SEPD and Arps hyperbolic equation, the Doung's with the Arps hyperbolic equation, and the Arps super hyperbolic combined with the Arps hyperbolic decline equation. He found that the SEPD and Arps hyperbolic equation gave the most conservative results of all the methods in the study, even if there was insufficient data available. This equation can also work without enough boundary-dominated flow (BDF) data available.
Hu et al. [27] studied DCA techniques for the Eagle Ford and Austin Chalk reservoirs. They found that in the case of the Eagle Ford reservoir, the MHD and the Doung's model provided the highest EUR estimations and the two lowest matching errors, while the PLE decline model with D ∞ 0 produced the lowest EUR estimates with the highest matching errors in all cases. In another study, according to the results of goodness of fit (R 2 and N-RMSE), the hyperbolic model fits well with aggregated well data and with individual wells [1]. Further, in their study, Hu et al. [27] explain that the LGM and PLE model with D ∞ = 0 gave production projections neither too positive nor too traditional with modest matching errors. Therefore, they recommend both the MHD and Doung's model for this reservoir. However, Zhang et al. [11] developed the EEDM and verified their model using field data from the Eagle Ford. They found this model to be more rigorous in that it will include the effects of interference among adjacent fractures, variable permeability, and discontinuous pressure distribution, which are difficult to capture and model with other DCA methods [11]. In the case of the Austin Chalk reservoirs, all DCA methods resulted in fairly similar EUR forecasts and matching errors; hence, any method can be used [27]. Figure 7, which uses estimated production data versus time values, indicates that using the R 2 values as a goodness of fit to determine the accuracy of the different decline models, the SEDM, followed by the LGM, EEDM, PLE, Doung's decline model and, lastly, the hyperbolic decline model would predict the EUR accurately.  During their case study analysis, Paryani et al. [13] found that the LGM, PLE, and Doung's models overcame Arps limitations to a certain degree. The PLE model always predicted the lowest forecasts of all the models with the most conservative production forecasting and reserve estimation. Doung's model performed the best for longer when less noisy production data was available; however, erratic EUR was observed, which indicates that this model requires further improvements [13]. The LGM gave reasonable EUR estimates when compared to the Arps model. There was an 81% fit of the wells' past production rate and cumulative production. The LGM also appears most effective at historically matching past production and predicting finite reasonable EUR. However, Tan et al. [4] found that due to the constraints of K and the vanishing production rate at infinity time, the LGM provides a finite estimate of EUR. They also determined by using normalized and logarithmic rate-time residuals that the limitations of the Arps model are overcome and accuracy improves in cases of unconventional reservoirs.

Conclusions
Shale gas reservoirs have become an essential source for providing natural gas globally and the process of hydraulic fracking has been used in the extraction of shale gas. During the fracking process, there are different flow regimes, which occur during the life cycle of SGRs, these being fracture linear flow, fracture boundary flow, matrix linear flow, and matrix boundary flow. They are significant because they impact both the production and decline behavior of SGRs.
Based on previous studies conducted, it was found that the Arps hyperbolic decline, the MHD and Doung's models provided the best fit with production data. However, contrary to the reviewed studies when estimated production data was used in the evaluation process for the basis of this paper, using the goodness-of-fit technique, the PLE and Doung's decline models aligned the best with the production data compared to the other models.
It is evident from the accuracy assessment decline curve modelling impacts the EUR of SGRs, and it was observed that all decline models yield a different EUR result, which is either over or underestimated. Studies have revealed that the production time significantly impacts the EUR depending on which decline model is being used. When each model was assessed for accuracy once again using the goodness-of-fit technique, the results indicated the SEDM, followed by the LGM, EEDM, PLE, Doung's decline model and, lastly, the hyperbolic decline model align with the production data.
It is evident from the decline curve evaluation there are advantages in using the current DCA models; however, they also have limitations associated with them, which have to be addressed. Therefore, the next step will be to evaluate the use of the hybrid model in evaluating the decline of SGR.