Methods of Decline Curve Analysis for Shale Gas Reservoirs

: With help from horizontal wells and hydraulic fracturing, shale gas has made a signiﬁcant contribution to the energy supply. However, due to complex fracture networks and complicated mechanisms such as gas desorption and gas slippage in shale, forecasting shale gas production is a challenging task. Despite the versatility of many simulation methods including analytical models, semi-analytical models, and numerical simulation, Decline Curve Analysis has the advantages of simplicity and efﬁciency for hydrocarbon production forecasting. In this article, the eight most popular deterministic decline curve methods are reviewed: Arps, Logistic Growth Model, Power Law Exponential Model, Stretched Exponential Model, Duong Model, Extended Exponential Decline Model, and Fractural Decline Curve model. This review article is dedicated to summarizing the origins, derivations, assumptions, and limitations of these eight decline curve models. This review article also describes the current status of decline curve analysis methods, which provides a comprehensive and up-to-date list of Decline Curve Analysis models for petroleum engineers in analysis of shale gas reservoirs. This work could serve as a guideline for petroleum engineers to determine which Decline Curve models should be applied to different shale gas ﬁelds and production periods.


Introduction
Shale gas has become an important energy production component due to the availability of horizontal well drilling and hydraulic fracturing techniques since the late 1980s. Figure 1a shows that currently natural gas accounts for one-quarter of the world energy supply , and Figure 1b illustrates that shale gas accounts for about half of the natural gas produced in the USA and could increase to 70% in 2040, as predicted by the Energy Information Agency of the United States.
Shale formations are characterized by natural fractures (Figure 2a), low-permeability rocks, and nanopores associated with organic matter, clay minerals, carbonate, and various silts (Figure 2b). In shale, gas content is mainly composed of free gas in the nanopores and adsorbed gas in the rocks and other components such as carbon dioxide, nitrogen, and hydrogen sulfide. The co-existence of these components with different size and properties has posed a major challenge to investigating shale. (b) High-resolution images of Montery shale, California, USA [4]. A, B, C/D are the scan images under different resolutions for the same shale core sample. A, B, and C are for the same core sample, under different resolutions, 10μm, 3 μm and 2 μm respectively. D is for another core sample under 2 μm resolution. Arrows in B and C mark pores less than 80 nm. Arrows in D mark carbonate micrite.
To simulate gas production in shale, complex gas transport mechanisms and the presence of natural and hydraulic fractures need to be considered, which makes the simulation of shale gas  (b) High-resolution images of Montery shale, California, USA [4]. A, B, C/D are the scan images under different resolutions for the same shale core sample. A, B, and C are for the same core sample, under different resolutions, 10μm, 3 μm and 2 μm respectively. D is for another core sample under 2 μm resolution. Arrows in B and C mark pores less than 80 nm. Arrows in D mark carbonate micrite.
To simulate gas production in shale, complex gas transport mechanisms and the presence of natural and hydraulic fractures need to be considered, which makes the simulation of shale gas  [4]. A, B, C/D are the scan images under different resolutions for the same shale core sample. A, B, and C are for the same core sample, under different resolutions, 10µm, 3 µm and 2 µm respectively. D is for another core sample under 2 µm resolution. Arrows in B and C mark pores less than 80 nm. Arrows in D mark carbonate micrite.
To simulate gas production in shale, complex gas transport mechanisms and the presence of natural and hydraulic fractures need to be considered, which makes the simulation of shale gas reservoirs a challenging task. There are many gas transport mechanisms coexisting during the gas production process, such as gas diffusion, gas slippage, and gas desorption. The complexity of these mechanisms makes classical free gas diffusivity equations inadequate to model the gas transport in shale. Existing simulation models for gas transport in shale have not fully captured the complex mechanism, and a new comprehensive and reliable gas transport model still needs to be constructed and studied [5][6][7].
With meticulously measured petrophysical information [8][9][10][11], the analytical, semi-analytical, or numerical solutions based on an accurate reservoir model could lead to reasonable production forecast. However, it is expensive and time-consuming to perform petrophysical experiments in order to build such a sophisticated reservoir model, because a lot of geological and petrophysics data need to be gathered and analyzed before being implemented into the reservoir model. For instance, test wells need to be drilled, and multiple well logging processes need to be performed in order to obtain the necessary input (e.g., formation location, depth, resistivity, permeability, fracture properties, gas properties) for the reservoir model. These processes are expensive and time-consuming, and some logging techniques are not sensitive enough to measure the multi-scale properties in shale. In contrast, DCA only requires the production data, and little reservoir information is needed. DCA is efficient with adequate accuracy to meet industrial needs. In the past few decades, several DCA models have been proposed and benchmarked with commercial reservoir simulators or shale gas production data before being applied to more shale gas reservoirs.
In this work, eight most popular DCA models for shale gas reservoirs will be reviewed. These models start from the classical Arps model and end with the new Fractional Decline Curve model. There have been previous review works or comparison works for different DCA models, such as Kanfar and Wattenbarger [12], Joshi and Lee [13], Lee [14], and Hu et al. [15], but there are a few new models proposed during the past several years and there is a need to include new applications of these DCA models in an updated review. This work will provide the most current theory and applications of DCA models in shale gas reservoirs.

Origins of DCA
In essence, DCA models are regressions for historical production data. When the production is controlled, the application of DCA is limited because the production rates are mostly constant. As the exploration and production techniques are further developed, there have been more capacity wells drilled; the production rates of these wells usually show a rapid drop at the early stage, then the decreasing rate of the production slows down at late production time. At late production time, these wells require more sophisticated DCA models to simulate the production decline trend.
Ralph Arnold is one of the pioneers who first applied data analysis techniques to estimate the oil reserves. The earliest reference to DCA was made by Arnold and Anderson [16], who assumed that the production rate at any time is a constant fraction of the rates at the preceding time. In Arnold [17], the author states that "Twenty years ago the estimate of oil reserves was computed generally by calculating the contents of the supposed reservoir rock from data regarding thickness, extent, etc., guessing at the saturation and percentage of recoverable oil, and finally arriving at a very rough approximation of the desired information. Nowadays, thanks to the great mass of data available and to the perfection of methods of computation, more accurate results are obtained. If the production of a given well over a reasonable period (even a period of days in some instances) is known, the future production of the well by years and in totality can be computed with remarkable accuracy, and thus arrive at the recoverable reserve for this well and its surrounding area." Arnold did not specifically use the terminology "Decline Curve Analysis," but "great mass of data" and "methods of computation" are the exact input data for DCA. The main idea introduced in Arnold [17] is exactly what DCA does to forecast the production. Since then, there have been more applications of DCA models for production forecasting [18,19]. For a complete review of early DCA work, refer to [20], where the author reviewed all the major DCA models developed before 1945.

Arps Decline Model
The classical DCA model was proposed in Arps [20], where the author proposed a hyperbolic function with three parameters to simulate the decline of flow rate. In the Arps model, bottomhole pressure is fixed, the skin factor is constant, and the flow regime is boundary dominated flow. To derive the DCA model, the concept of loss ratio was first introduced. The loss ratio is defined as the ratio between the production drop of the current time step and that of the previous time step. Based on observations, Arps proposed two different scenarios of loss ratio. The first scenario is to assume that the loss ratio is a constant, where b is a positive constant. Integrating Equation (1), we get the exponential decline functions as follows: where q i is the initial rate, in bbl/day. The second scenario is to assume that "first differences of the loss ratios are approximately constant", i.e., d( The double integration of Equation (3) allows us to obtain the rate-time relationship for hyperbolic decline: where a i is the initial loss ratio.
Based on this decline model, the curve has a slope − 1 b on a log-log paper. During the fitting process, these parameters can be determined by calculating the derivatives of production data with respect to time.
Although Arps model is simple and fast, it often fails to accurately fit the decline curve of unconventional reservoirs and predict the estimated ultimate recovery (EUR) [21]. The Arps model often tends to overestimate the EUR for shale gas wells because it assumes that a boundary-dominated flow regime prevails [22,23]. Since most shale gas wells rarely reach the boundary-dominated flow regime, the Arps model cannot be applied directly to shale gas reservoirs without significant modifications [24]. Following the Arps decline curve, various new models were proposed to better model the boundary-dominated flow phase, and will be presented below.

Modified Hyperbolic Decline Model
Arps model was originally designed for pseudo-radial boundary-dominated flows (BDF), which occur in medium-high-permeability formations. However, for tight-gas shale reservoirs, the fracture-dominated flow regimes dominate and the late-time-flow regime is rarely reached. Meanwhile, the Arps decline model leads to an unbounded and unrealistic estimate of EUR for b ≥ 1 [25]. To resolve this issue, Equations (2) and (4) need to be applied in a piecewise fashion, which means using different DCA curves for different production time intervals. The hyperbolic decline function, Equation (4), is applied in the early stage; after the decline rate reaches a certain value, the exponential decline function, Equation (2), is used. This process could be achieved by applying computer programs to determine the switch point, which is the point when the decline rate is smaller than a certain limit (5% is often used). To generalize Equations (2) and (4) in one expression so that there is no need to change formulas for different time intervals, Robertson [26] proposed the following Modified Hyperbolic Decline (MHD) model: where a, β, and N are all positive constant. Note that when β → 1 , Equation (5) is equivalent to Equation (4); when N → ∞ , Equation (5) is equivalent to Equation (2) ( [26]). Another equivalent way to express the modified hyperbolic model was proposed in Seshadri and Mattar [27], where hyperbolic decline in the early life of a well is shifted to exponential decline in the late life by imposing a predetermined decline rate, D * . Once the decline rate reaches D * , the decline curve switches from a hyperbolic function to an exponential function. The equivalent MHD model is written as where q i is the initial rate, in bbl/day, and D 1 and D 2 are decline rates. Meyet et al. [28] pointed out that the selection of the decline rate D * can be empirically determined and has no physical basis. The authors also showed that the EUR forecasts by MHD are often greater than those of other DCA methods. Paryani et al. [29] also showed that the MHD model provides over-optimistic estimates of reserves and longer remaining life of shale oil wells.

Power Law Exponential Decline Model
As mentioned in Section 2.2, Arps [20] introduced the so-called "loss-ratio" as in Equation (1). When b is a constant, the exponential decline, Equation (2), is obtained. In Ilk et al. [30], the authors proposed another approach to formulate parameter b: where D ∞ , D 1 are the decline constant at infinite time and initial time, andn is the time exponent. The physical interpretation of Equation (7) is that the loss ratio can be approximated by a decaying power function with constant behavior at late production time. By substituting Equation (7) into Equation (1) and integrating, the following production rate is derived: whereq i is the rate "intercept",D i is the initial decline constant, D ∞ is the decline constant at infinite time, andn is the time exponent. In addition, the parameters D andD i are defined as follows: The model defined by Equation (8) is called the Power Law Exponential (PLE) Decline Model. It is based on the Arps decline curves and uses the power law to approximate the production rate. This model is developed specifically for shale gas wells.
The PLE model has extra variables to account for both transient and boundary-dominated flows. However, four unknowns in this model,q i ,D i , D ∞ , andn, cause many degrees of freedom while fitting real-world data, which can be impractical. There are many methods to solve the PLE model, and multiple solutions are acceptable. For example, the interpolation algorithms could be written using regression techniques like a classical least squares method.
However, Equation (10) only holds true during transient flows [31], which could be determined by using the derivate plot of the flow rate where the slope is −1. Therefore, attention is needed to select appropriate data when using Equation (10).
McNeil et al. [32] showed that the PLE model matches production data in transient and boundary-dominated regions without being hypersensitive to remaining reserve estimates, and the authors prefer the PLE model over the traditional Arps method with insufficient boundary-dominated data or a very long transient period. Also, the authors showed that the PLE model provides consistently reasonable results for the EUR. Seshadri and Mattar [27] pointed out that the PLE model can model transient radial and linear flows. Their results show that the PLE model does not appear to provide any significant tangible advantage over the MHD model. Kanfar and Wattenbarger [12] proved that the PLE 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. Meyet et al. [28] illustrated that the PLE model varies the least with respect to length of production history for all wells, and the PLE and LGA models usually return similar results. Paryani et al. [29] showed that the PLE model fits the historical production data 67% of the time, and provided the lowest forecasts among all the DCA models in their study and thus the most conservative forecast.

Stretched Exponential Decline Model
To avoid the disadvantages of the Arps decline model and use the relatively easier access of large dataset of well productions, Valko [33] and Valkó and Lee [34] proposed the Stretched Exponential Decline Model (SEPD), in which they assume that the product rate satisfies the stretched exponential decay: Integrating Equation (11) yields where τ is the characteristic time constant and n is the exponent. Valkó and Lee [34] mentioned that a natural interpretation of this model is that the actual production decline is determined by a great number of contributing volumes. All these volumes have exponential decay rates, but with a specific distribution of characteristic time constants (τ). To construct Equation (11), nonlinear Equation (13) need to be solved in order to find the undetermined parameters, which is time-consuming: where r 21 and r 31 are the cumulative production ratio for two and three years, respectively. Γ(t) and Γ(s, x) are the gamma function and incomplete gamma function, respectively. They are defined as follows: Akbarnejad-Nesheli et al. [21] showed that the SEPD model is advantageous for combining the concave and convex portions of decline curves without increasing the number of model parameters, and could provide a finite (bounded) value of EUR without cutoffs in time or rate. Zuo et al. [35] also illustrated that the SEPD model provides a bounded EUR rather than an infinite value; moreover, the authors pointed out that the SEPD model models transient flow rather than boundary-dominated flow, and requires a sufficiently long production time (usually >36 months) to accurately estimate the parameters τ and n. Also, the construction of the SEPD model requires solving complicated nonlinear equations to determine unknown parameters (Equation (13)).

Duong Model
The Duong model [36] is introduced based on one empirically derived rule, that is, the log-log plot of q/G p (G p is the cumulative gas production) vs. t forms a straight line. The author conjectured that Based on Equation (16), the author derived the formula for well production rate and cumulative production as follows: where a is the intercept coefficient from Equation (16) and m is the slope in the log-log plot. Kanfar and Wattenbarger [12] showed that the Duong model is more accurate for linear flows and bilinear-linear flows. Meyet et al. [28] showed that the EURs determined with PLE and Duong model vary the least with respect to the length of production history for all wells among all of the DCA methods in their study, and other DCA methods tend to converge towards the modified Duong model and PLE model. Furthermore, the Duong model tends to provide the most conservative results. Zuo et al. [35] pointed out that if m and a in Equation (17) are within certain ranges, the gas flow rate should decrease monotonically, and q i determination in the model may lead to unreliable results if the production history is shorter than 18 months. Paryani et al. [29] fitted well with 51% of the historical production data, and the Duong model fits better with longer and less noisy historical production data. In Wang et al. [37], the authors proposed a new empirical method and compared it with the SEPD model and Duong model, and concluded that the SEPD underestimates EUR and the Duong model overestimates the ultimate recoveries.

Logistic Growth Model
Logistic Growth Model was once used by Spencer and Coulombe [38] to model the regrowth of livers. The rat livers in their experiment were reduced to one-third of their original sizes, and appeared to regenerate hyperbolically. Clark et al. [39] adopted this model for shale gas reservoirs with extremely low permeability, and developed the Logistic Growth Model (LGM) as an empirical method to forecast the gas production. The cumulative production in this model is up to a maximum carrying capacity and there is no further growth after reaching this maximum value. The cumulative production for forecast is calculated using the following equation: where Q is cumulative production, in bbl, K is carrying capacity, a is a constant of t n , n is the hyperbolic exponent, and t is the time, in days. K can be determined by volumetric methods or curve fitting methods. n controls the curvature of the decline curve, and has values between 0 and 1. The constant a is the time to the power n at which half of the carrying capacity has been produced, which is easy to see when you take the limit of t n going to a in Equation (19). In fact, when t n goes to a, Equation (19) goes to K/2. The production rate of this model can be obtained by differentiating Equation (19), which is shown below: The relevant parameters (K, n, andâ) can be obtained by either optimizing the parameters with a numerical scheme such as the least squares method or linearizing the equations and plotting the data.
The major advantage of the LGM is that the reserve estimate is constrained by the parameter K as well as the production rate, which terminates at infinite time [39]. However, an upward inflection in the curve would occur if n > 1 [28]. The main assumption in this model is that the entire reservoir can be drained by a single well over a sufficiently long period [29].

Extended Exponential DCA Model
Another nice piece of work to avoid using piecewise functions was done by Zhang et al. [40]. The authors combined the exponential form of the decline equation proposed by Fetkovich [41]: with a newly proposed empirical formula to calculate index a: where β l , β e are constants, combining Equations (21) and (22) to obtain: This new method is called Extended Exponential Decline Curve Analysis (EEDCA). In Zhang et al. [40], this model is validated by comparison with the Arps model and SEPD model, as well as numerical simulations. The comparison showed that EEDCA projections matched the production data. The authors also tested the model with hindcasts of 85 shale wells, before applying this model to fit seven shale reservoir wells in the Eagle Ford. Note that if we substitute Equation (22) into Equation (21), we can get q = q 0 e −(β l +β e e −t n )t = q 0 e −β l t+β e te −t n , which is different to the PLE model in Equation One of the advantages of the EEDCA model is that these two new parameters, β e and β l , once calibrated using known production data, can capture both the early production profile and the late production profile. Early on, the term β e e −t n has the dominant impact, but as time goes on, the impact of β e e −t n reduces, and the other term, β l , becomes dominant. Thus, there is no need to solve at least two sets of equations as MHD (Equation (6)).

Fractional Decline Curve Model
For unconventional reservoirs, the flow decay rate is slower than exponential decay so that the decline curve has a long tail, which is statistically caused by anomalous diffusion phenomena. Ordinary diffusion is an important process described by a Gaussian distribution. In one dimensional diffusion, assuming a particle is initially at x = 0 when t = 0, the probability density function P(x, t) at time t is A feature of this process is the linear relationship between the mean square displacement and the time, which is x 2 = 2t. However, for anomalous diffusion, we see that x 2 ∝ t γ , γ = 1 or x 2 might be a divergent integral for t = 0. The latter process is called Lévy flight. Chaves [42] proposed a generalized form of Fick's Law wherein he replaced the gradient with an α-th order (1 ≤ α ≤ 2) Riemann-Liouville fractional derivative (whose definition will be given shortly) and obtained a fractional order diffusion equation. Then he rewrote this fractional diffusion equation in an anisotropic medium, in which case it generates asymmetric Levy statistics instead of normal Gaussian distribution. It describes Lévy flight very well. From the point of view of statistical physics, classical diffusion is caused by the Brownian motion of the particles. The spatial probability density function evolving over time governs Brownian motion, and is a Gaussian distribution whose variance is proportional to the first power of time. In contrast, there have been many experiments [43,44] showing that the anomalous diffusion is characterized by its variance behaving like a non-integer power of time.
For example, Nakagawa et al. [45] found a significant difference between the actual diffusion profile of contaminants underground and the theoretical profile predicted by a conventional diffusion equation, which was pointed out by Adams and Gelhar [43]. The corresponding equation predicting the anomalous diffusion is called a fractional diffusion equation. There have been many studies on anomalous diffusions and the corresponding fractional diffusion equations. For example, Zhou and Selim [46] showed that subdiffusion in porous media can be characterized by a long-tailed profile in the spatial distribution of densities; Metzler and Klafter [47] illustrated a fractional diffusion equation with respect to a non-Markovian diffusion process; Chaves [42] derived the same equation to describe Levy flights; Nigmatullin [48] first considered the fractional initial boundary value problem and derived several theoretical and numerical results; Luchko et al. [49,50] obtained the maximum principle and the unique existence of the generalized solution; Sakamoto and Yamamoto [44] studied the uniqueness of weak solutions and the asymptotic behavior; Luchko and Zuo [51] studied the explicit form of the forward solutions for subdiffusion equations with Dirichlet or Neumann boundary conditions; Chen and Raghavan [52] constructed a one-dimensional, fractional-order transient diffusion equation to model diffusion in complex geological media; The pressure distribution was derived in terms of the Laplace transformation and the Mittag-Leffler function; Holy and Ozkan [53] presented a fractional diffusion-based approach for modeling two-phase flow in fractured porous media, such as that encountered in unconventional reservoirs. These fractional diffusion models have shown better accuracy in fluid flow mechanisms in a complex porous medium, such as groundwater contamination and shale gas production.
Based on anomalous diffusions, Zuo et al. [35] introduced a new Fractional Decline Curve (FDC) model based on anomalous diffusions to match the flow rate and to predict the future production rate. In petroleum engineering, the following fractional diffusion equation is the most popular to model the fluid dynamics in porous media [54]: where ∂ α t is the Caputo fractional derivative, which is defined as follows [55]: where p xx is the 1-D Laplace operator defined as p xx = ∂ 2 p ∂x 2 , and n − 1 ≤ α ≤ n. The classical Mittag-Leffler function [56] is a special case (β = 1) of the following definition of the general Mittag-Leffler function: where z ∈ C, α, β ∈ R. It is an entire function in z with order 1 α and type one. Based on Equation (27), Zuo et al. [35] proposed the following DCA model: Energies 2018, 11, 552 10 of 18 In Figure 3, the decline profile for the Mittag-Leffler function for λ = π 2 is illustrated. The long-tail behavior is observable for anomalous diffusion profile. This behavior is very different from Gaussian diffusion. For anomalous diffusions, the flow rate decays quickly in the beginning but declines more slowly at later stages. When α gets smaller, the long-tail profiles become more obvious. Note that α = 1 corresponds to the Gaussian diffusion.
where ∈ , , ∈ . It is an entire function in with order and type one.
Based on Equation (27), Zuo et al. [35] proposed the following DCA model: In Figure 3, the decline profile for the Mittag-Leffler function for 2 λ π = is illustrated. The longtail behavior is observable for anomalous diffusion profile. This behavior is very different from Gaussian diffusion. For anomalous diffusions, the flow rate decays quickly in the beginning but declines more slowly at later stages. When α gets smaller, the long-tail profiles become more obvious. Note that 1 α = corresponds to the Gaussian diffusion. There is the following asymptotic behavior of Mittag-Leffler functions [55]: If 0 2 α < < and μ is an arbitrary real number such that min{ , } 2 πα μ π π α < < , then for an integer 1 p ≥ , we have the following expansions, where ( ) Γ ⋅ is the gamma function and arg( ) z is the argument of complex number z . Equation (29) will be used to estimate the possible ranges of parameters λ and α . Four steps are needed to reconstruct the FDC model: Step 1: Data Preconditioning. During the production period, some wells might be shut off so the production data will drop significantly on the following day. In this step, the production data are filtered out during the time when the well was shut off.
Step 2: α ,λ, and m Determinations. In Equation (29), for long enough t , the behavior of function so that the log-log plot of 1/ q vs. t form a straight line, with slope α and an intercept . Solving for α and λ to obtain a close approximation of these two parameters. m There is the following asymptotic behavior of Mittag-Leffler functions [55]: If 0 < α < 2 and µ is an arbitrary real number such that πα 2 < µ < min{π, πα}, then for an integer p ≥ 1, we have the following expansions, where Γ(·) is the gamma function and arg(z) is the argument of complex number z. Equation (29) will be used to estimate the possible ranges of parameters λ and α. Four steps are needed to reconstruct the FDC model: Step 1: Data Preconditioning. During the production period, some wells might be shut off so the production data will drop significantly on the following day. In this step, the production data are filtered out during the time when the well was shut off.
Step 2: α, λ, and m Determinations. In Equation (29), for long enough t, the behavior of function so that the log-log plot of 1/q vs. t form a straight line, with slope α and an intercept log(λΓ(1 − α)). Solving for α and λ to obtain a close approximation of these two parameters. m will be taken as the largest production rate of this group of data. The iteration and trial and error m methods are used to obtain all the possible α, λ values and calculate the l 2 -norm of the difference between the actual production data and the approximation value. The parameter triple (α, λ, m) corresponding to the least l 2 -norm will be our solution.
Step 3: Validation of these parameters with the actual data.
Step 4: EUR forecast comparison. The EUR is calculated using the following formula: In Zuo et al. [35], the FDC model was validated with a numerical shale gas model and was compared with the Arps model. The comparison results showed that the FDC model well matches the production rate produced by the commercial software CMG and has a much smaller error than the Arps model. Zuo et al. [35] also applied the new model for five wells in Fayetteville shale and obtained results that matched the measurement data well.
Wang et al. [37] stated that the FDC model requires iterative programming to optimize the original parameters obtained by production, and EUR is calculated by daily production accumulation, which makes the FDC rather complicated and inconvenient.
To summarize all eight DCA models for a quick reference of readers, Table 1 lists the name of each method, its DCA expressions, and the related references.   [6,35]

Probabilistic Decline Curve Model
Despite the simplicity and easy application of DCA methods as compared to other production forecast tools such as analytical methods, semi-analytical methods, and numerical methods, DCA methods have some shortcomings due to their nature of the deterministic inversion procedures. These shortcomings include the strong dependency on the smoothness of the production data. If there are many outliers in the production data, the determination of the unknown parameters in DCA models could become difficult. This difficulty becomes worse in shale gas reservoirs, as the production data of shale gas reservoirs are more complex and more random than those of conventional reservoirs. Common disadvantages of deterministic inversion methods include the following [57]: (1) The nonlinear minimization requires explicit information about model derivatives; (2) no explicit information is provided about the uncertainty of the interpretation results; (3) inversion-based interpretation results are deeply biased by the initial search point in model space. To minimize these disadvantages and prevent large errors in deterministic DCA models, various probabilistic DCA methods can be used to forecast the production in shale gas reservoirs based on the limited production data available. Although this review focuses on deterministic DCA models, probabilistic DCA methods are still briefly reviewed below so that interested readers could use this as a reference list.
Bayes's Theorem [58] has been widely used for uncertainty quantification [59][60][61][62]. Gong et al. [59] introduced a Bayesian probabilistic methodology using Markov chain Monte Carlo (MCMC) coupled with the standard Metropolis-Hasting (MH) algorithm [63,64]. The Arps decline curve model was applied to quantify the uncertainty in production forecasts and reserve estimates. Applying the new method to 197 Barnett shale gas wells, the authors concluded that their methodology can quantify the uncertainty in hindcasted production with a narrower P90-P10 interval and higher computational efficiency than the modified bootstrap method [65], which was originally proposed by Jochen and Spivey [66]. A Bayesian machine learning method with the MCMC coupled with the MH algorithm was proposed in Fulford et al. [60], and was used to forecast production in shale wells by utilizing the transient hyperbolic model introduced by Fulford and Blasingame [67]. However, the behavior of b and D parameters as a function of time is only valid for homogeneous reservoirs with equal fracture half-length and spacing, not suitable for heterogeneous reservoirs with varying fracture half-length and spacing [68]. In addition, the standard MH algorithm has a lower acceptance ratio [69] than other MCMC algorithms such as delayed rejection [70], adaptive Metropolis [71], and the combination of delayed rejection and adaptive Metropolis [72]. A good summary and discussion of the advanced MCMC algorithms can be found in Goodwin [73]. To improve the MH algorithm, Yu et al. [74] proposed a new probabilistic approach with the Bayesian methodology combined with MCMC sampling and a FDC decline curve model to increase the efficiency and reliability of the uncertainty quantification. The authors used an adaptive Metropolis algorithm to replace the MH algorithm and applied their model to multiple gas wells in Fayetteville shale. Recently, de Holanda et al. [75] constructed a physics-based DCA model accounting for linear flow and material balance in horizontal multi-stage hydraulically fractured wells. In that work, the model was applied to a large dataset in a workflow that incorporates heuristic knowledge into the history matching and uncertainty quantification by assigning weights to rate measurements. The uncertainty quantification was performed via a Bayesian approach with hindcasts, and applied to 992 gas wells from the Barnett shale.

Comparisons of DCA Models with Field Data
The above reviewed eight DCA models are popular candidates to be used for shale gas production forecasting. However, different models provide different predictions for production data from various shale formations. In this section, several comparison studies are illustrated to show the differences among these DCA models.
Kanfar and Wattenbarger [12] compared five DCA methods explained in our study: Arps, PLE, SEPD, Duong, and LGM models. These models were compared by matching field data from Barnett shale and Eagle Ford shale. The field data are Barnett shale gas Well 314 and Eagle Ford shale gas Well 204. Figure 4 shows the matching results for these five models, which shows that for Well 314 in Figure 4a, the PLE and LGM methods have close results while Arps returns the most optimistic estimate and SEPD is the lowest estimate. For Well 204 in Figure 4b, the SEPD method has the lowest estimate and Arps has the highest estimate. From these applications, the authors concluded that: if a well is within a linear flow regime, the most accurate EUR result could be obtained by using Arps, PLE, or Duong; for a well within the bilinear flow, the best DCA models are the PLE or Duong model; if the well production shows boundary-dominated flow, the best DCA model is the PLE method.  Joshi and Lee [13] compared three DCA models-Arps, SEPD, and Modified Duong-by matching field data from Barnett Shale. Figure 5 illustrated that, for two cases in Figure 5a,b, the Arps model overestimated the production compared to historical data, while SEPD and Modified Duong provided acceptable results. Joshi and Lee [13] compared three DCA models-Arps, SEPD, and Modified Duong-by matching field data from Barnett Shale. Figure 5 illustrated that, for two cases in Figure 5a,b, the Arps model overestimated the production compared to historical data, while SEPD and Modified Duong provided acceptable results.  Joshi and Lee [13] compared three DCA models-Arps, SEPD, and Modified Duong-by matching field data from Barnett Shale. Figure 5 illustrated that, for two cases in Figure 5a,b, the Arps model overestimated the production compared to historical data, while SEPD and Modified Duong provided acceptable results.
(a) (b) Figure 5. Comparison results for two well groups from Barnett shale [13], where EOH stands for the end of history used in the matching, and EOP is the end of production data. In both cases 36 months of historical data were used for matching. Here we also studied the differences of Arps, Duong, FDC, and SEPD by matching the Fayetteville shale gas production data with these four models. The comparison results of the flow rate and EUR are plotted in Figure 6. They show that Arps returns the most optimistic estimate, and SEPD returns the most pessimistic one, while FDC and Duong return similar intermediate matching results. Here we also studied the differences of Arps, Duong, FDC, and SEPD by matching the Fayetteville shale gas production data with these four models. The comparison results of the flow rate and EUR are plotted in Figure 6. They show that Arps returns the most optimistic estimate, and SEPD returns the most pessimistic one, while FDC and Duong return similar intermediate matching results.  Joshi and Lee [13] compared three DCA models-Arps, SEPD, and Modified Duong-by matching field data from Barnett Shale. Figure 5 illustrated that, for two cases in Figure 5a,b, the Arps model overestimated the production compared to historical data, while SEPD and Modified Duong provided acceptable results.
(a) (b) Figure 5. Comparison results for two well groups from Barnett shale [13], where EOH stands for the end of history used in the matching, and EOP is the end of production data. In both cases 36 months of historical data were used for matching. Here we also studied the differences of Arps, Duong, FDC, and SEPD by matching the Fayetteville shale gas production data with these four models. The comparison results of the flow rate and EUR are plotted in Figure 6. They show that Arps returns the most optimistic estimate, and SEPD returns the most pessimistic one, while FDC and Duong return similar intermediate matching results.

Discussion
Despite challenges to shale gas reservoir modeling due to complicated transport mechanisms and the existence of fracture networks, oil companies have not slowed down on shale gas investment and production using horizontal well drilling and hydraulic fracturing techniques. Because building a reservoir model usually requires drilling test wells and performing well logging measurements, many small oil companies may not have the budget to do so. Even for large oil companies, building a reservoir model is not worthwhile for evaluation of small-scale oil fields. In those cases, comprehensive numerical simulation methods are likely impractical. In order to forecast the production of these reservoirs, DCA is one of the most convenient and practical techniques.
With the rapid increase in shale gas production over the past 30 years, there have been numerous production data for shale gas reservoirs. Many different DCA models have been constructed to model the shale gas production rate, from the classical Arps to the latest FDC model; each has its advantages and shortcomings. In practice and in all existing commercial DCA software, most of these eight DCA models are implemented and open to be used. They are also extensively compared by clients, shown in this review. Most of the deterministic DCA models are empirical and lack a physical background so that they cannot be used for history matching of the reservoir properties. Up until now there have been only a few studies that probed the relationships between DCA parameters and reservoir properties, such as permeability, porosity, fracture length, and fracture widths. Zuo et al. [35] discovered that, for shale gas wells in the Fayetteville formation, the DCA parameters are all within a certain small range. For example, for Fayetteville shale, parameter α defined in Equation (28) lies in the range from 0.65 to 0.75, which indicates that this parameter might have some relationship with the physical parameters. Zhang et al. [40] showed that some recommended empirical coefficients can be determined for specific shale formations such as Eagle Ford, which may help, particularly in the early stages of development when little production data are available. In the future, once more information about the fractures and reservoirs becomes available, more detailed investigations about the links between DCA parameters and reservoir properties can be performed.
Nowadays, with the advancement of machine learning and data analysis techniques, the combination of deterministic DCA models with machine learning techniques could also be an interesting future trend of DCA model applications.

Conclusions
Because of their simple mathematical expressions, fast computation speed, and low dependence on the geophysical data of the reservoirs, DCA models have been widely applied to to forecast shale gas production using data. As the simplest technique for shale gas production forecasting, DCA models do not need specific reservoir or fracture data input. DCA models do not require the construction of a complicated reservoir model. Using only the production data, DCA models can return accurate history matching with the production data by different type curves and forecast the EUR using interpolated curves. In this study, eight popular DCA models for shale gas reservoirs are reviewed, including their origins, formulations, and the types of reservoirs they fit. Their advantages and disadvantages have also been clearly presented.
As one of the earliest DCA models, Arps is designed for boundary-dominated flows but most shale gas reservoirs rarely reach the boundary-dominated flow regime, so the original Arps curve is not appropriate to simulate the shale gas reservoir production. If used directly on shale gas production, Arps tends to overestimate EUR.
To be compatible with two different flow regimes (i.e., the transient flows and the boundarydominated flows), two decline rates can be used for different time periods. By doing so, the modified Arps curve and MHD could return a good match for shale gas production by fitting the unknown parameters. However, the determination of the critical decline rate D * relies on empirical methods.
The SEPD model has a better performance for transient flows than for boundary-dominated flows. Because its expression has a finite limit, it has the advantage of providing a bounded value of EUR. The disadvantage is that SEPD requires a sufficiently large set of production data in order to obtain good determination of the unknown parameters. With few production data, SEPD usually return the lowest EUR. The modified Duong model is more accurate for linear flow and bilinear flow than other DCA models proposed before. However, if the production history is shorter than 18 months, the Duong model could return unreliable results for the EUR forecast. Most of the time, the Duong model overestimates the total EUR.
For shale gas reservoirs with a single well produced over a sufficiently long time, the SEPD model can be applied. As the SEPD model, LGM returns a finite estimate of EUR because of the constraints of carrying capacity K as well as the vanishing production rate at infinity time. SEPD is suitable for extremely low-permeability shale gas reservoirs and could provide less EUR than the Arps model. It usually yields similar results to the PLE model.
By considering both the early and late production profile, the EEDCA model does not require a switch from a transient model to a boundary-dominated flow model. It could result in smooth decline profiles. EEDCA model can be applied for shale gas reservoirs with both early and late production data.
Originating from the anomalous diffusion phenomena with a rigorous physical background, the FDC model is constructed based on the anomalous diffusion kernel function. It has historical dependence effects and could match the long tail behavior of the shale gas production well. For shale gas wells close to each other, the FDC model could have very similar parameter values.
To decrease the error caused by noisy data and better use the new production data, many probabilistic DCA models can be constructed to increase the accuracy of the deterministic DCA methods. When sufficient data are available, probabilistic DCA models can be applied to return reliable probability ranges.
The biggest disadvantage of all eight DCA models is that most of them are empirical without physical background with exception being FDC for some physical background. There are potential relationships between DCA parameters and reservoir properties, which need to be further investigated when more accurate field measurement of reservoir and fracture properties are available.