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 a ^ ) 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

which is different to the PLE model in Equation (8)

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^{−tn} 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 ∈ ℂ , α , β ∈ ℝ . 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:

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.

**Figure 3.** Decay rates comparison between Gaussian and anomalous diffusion for λ = π^{2} .

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 mE_{α,1}(−λt^{α}) is dominated by the second term

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.