Production decline analysis of shale gas based on a probability density distribution function

Production decline analysis is a simple and efficient method to forecast production dynamics of shale gas. The traditional Arps decline model is also widely used in the production decline analysis of shale gas, but an obvious error is often generated. Based on the Weibull and 𝜒 2 probability density distribution function, the monotonic decreasing production prediction equations of shale gas are established. It is deduced that recently, the widely used Duong model is essentially a Weibull probability density distribution model. Decline analysis results of production data from actual shale gas well and numerical simulations indicate that the fitting results of the Weibull (Duong) model and 𝜒 2 distribution model are better than the Arps model whose deviation of early data is large. For a shale gas reservoir with very low permeability, pressure conformance area is small and it is obviously influenced by fractures. Early shale gas production rate mainly contributed to by fractures declines quickly and the later production rate mainly contributed to by the matrix declines slowly over time. The production decline curve has obvious long-tail distribution characteristics and it is a better fit to the data with a 𝜒 2 distribution model. As for the increase of permeability, the fitting accuracy of the Weibull (Duong) model gradually becomes better than the 𝜒 2 distribution model. Research results provide theoretical guidance for choosing a reasonable production decline model of a shale gas reservoir with a different permeability.


Introduction
Empirical method, type-curve analysis and numerical simulation are the main methods of production decline analysis (Clarkson 2013;Zhang et al. 2015;Zheng et al. 2016;Li et al. 2018). Since the Arps decline model has been proposed, it has been widely used in conventional and unconventional oil and gas reservoirs (Ilk et al. 2008;Ikewun & Ahmadi 2012;Mishra et al. 2014). In most cases, exponential decline, harmonic decline and hyperbolic decline of the Arps model can fit most of the production data. However, the Arps model is mainly suitable for a flow regime influenced by a reservoir boundary. In an unconventional shale gas reservoir developed with a fractured horizontal well, pressure is difficult to propagate to the boundary because of ultra-low permeability. The Arps model used in the production decline analysis of a shale gas reservoir can easily generate some errors. An empirical method such as the Arps model is mainly based on statistical analysis of production data. Type-curve analysis and numerical simulation methods are based on seepage mechanism of shale gas and they can be used to verify the results of empirical model analysis.
Different medium including organic kerogen, inorganic matrix, natural micro fractures and hydraulic fractures exist in shale gas reservoirs (Curtis et al. 2010;Akkutlu & Fathi 2012). Early shale gas production rate is mainly contributed to by fractures and it declines quickly. After that, the shale gas production rate declines slowly because of the long-term contribution of gas in the matrix. It is important to look for a more accurate production decline model of shale gas. Besides the Arps model (Arps 1945), a stretched exponential decline model (Valko & Lee 2010), Duong model (Duong 2011) and fractional decline curve (FDC) (Zuo et al. 2016) have been proposed. To a certain extent, they improve the accuracy of model prediction. The Duong model has a good application effect in a reservoir with a fracture flow regime. It assumes that there is a log-log straight relationship between q/G p and time, which has been observed in a tight or shale gas reservoir. The FDC can describe long-tail behavior of a shale gas production profile where gas production rate declines rapidly in the early stage and then slowly over time. Miao et al. (2018) proposed declining models of shale gas in an early production period with self-diffusion characteristics and a late production period with surface diffusion characteristics. Zuo et al. (2016) and Miao et al. (2018) consider that there is a log-log straight relationship between 1/q and time after a certain time, yet sometimes the correlation coefficient of the straight line is not good.
In the process of shale gas production, gas rate and pressure are constantly changing. A conventional method is only suitable for conditions with constant rate or pressure, but not suitable for variable rate or pressure production. For variable rate and variable pressure problems, Blasingame's model (Blasingame et al. 1991) and the Agarwal-Gardner model (Agarwal et al. 1998) have been proposed and widely used (Sun 2015). They can effectively identify fluid flow state as well as calculate formation flow coefficient, and they can control radius, skin factor and control reserves. In recent years, these methods have been used in dynamic analysis of shale gas reservoirs (Wei et al. 2016). These models reflect the seepage mechanism of shale gas, but they are not convenient because of mathematical complexity. It is still of great significance to seek relatively simple analytical methods to analyse the production dynamics of shale gas.
Probability density function reflects the distribution law of probability f(x) with variable x. The relationship between production rate q(t) and time t also embodies some distribution law. Many mathematical methods involving probability distribution functions have been used to describe growth or decline laws of production rate and estimate uncertain reserves (Gong et al. 2011;Gonzalez et al. 2012;Mishra et al. 2014). Some of these models include the Gompertz, Richards, logistic and Weibull models (Zeide 1993;Ricklefs & Scheuerlein 2002;Tabatabai et al. 2005;Clark et al. 2011). Through analogy between probability density distribution function and production change function of shale gas, the monotonic decreasing production prediction equations of shale gas are established and the adaptabilities are researched in this paper. The new models get good fitting results when analysing production decline data of shale gas from actual gas wells and numerical simulations. Traditional Arps and other models cannot be completely abandoned in the production decline analysis of shale gas, yet in some cases the model in this paper is an alternative method for the shale gas reservoir with different permeability or long-tail characteristics of production profile.

Establishment of production decline equation
The Weibull probability density distribution function is If < 1, Weibull distribution belongs to heavy-tailed distribution.
In the interval of x from 0 to ∞, the integral of Weibull distribution function is F(x) = ∫ ∞ 0 f (x)dx = 1. Equation (1) can also be written in the following more concise form: where x is a distribution variable, whose distribution interval is 0 ∼ ∞, is the shape parameter that controls the distribution form and is the scale parameter that controls peak position and peak value. In the interval of x from 0 to ∞, the integral of f(x) also equals 1.
It can be seen that equations (1) and (2) are equivalent. The probability density curve of Weibull distribution is shown in figure 1. It shows that the curve shape is influenced by the value of . When ≤ 1, the curve monotonously decreases. Particularly when = 1, the curve is exponential. When >1, the curve is a single-peak distribution that first increases and then decreases. For the decline curve of shale gas production without a peak, this can be analysed using the Weibull model with a value of ≤ 1. For production of shale gas where q is daily production rate of shale gas, 10 4 m 3 /d; G R is recoverable reserve of shale gas, 10 4 m 3 . Comparing equation (3) with equation (4), it is known that f (x) = q∕G R , then Assuming Equation (6) is production forecasting model of shale gas based on a Weibull probability density distribution function.
Assuming (7) is converted into where t ≥1, whose initial time is t 1 . Considering initial time is t = 1, cumulative production rate of shale gas is calculated with an integral of equation (8), which is Cumulative shale gas production rate G p increases as time increases. When the time t reaches a certain value, G p can be approximately expressed as Equation (8) is the Duong model that has been widely used in recent years. In the process of model deduction, Duong states that the ratio of shale gas instantaneous production rate to cumulative production rate satisfies the following equation.
Shale gas instantaneous production rate formula equation (8) and cumulative production rate formula equation (10) can be deduced according to equation (11). Based on the Weibull probability density distribution function, Duong model can also be deduced in this paper. Therefore, the Duong model is essentially a Weibull probability density distribution model. Actually, equation (10) is an approximate expression of equation (9). Only the time t is high enough; equation (10) and equation (11) are satisfied. That is to say, the loglog curve of q/G p ∼ t is a straight line only after a period of production, which is shown in equation (11). According to statistical analysis of a large amount of production data in tight and shale gas reservoir, Duong found that the production data with linear flow characteristics (including linear flow and bilinear flow) can be described by equation (11), which provides a scientific basis for cumulative shale gas production rate formula of equation (10). In order to conveniently solve equation coefficients in equation (8), equation (10) rather than equation (9) is used in the following study.

Solution of production decline equation
Equation (6) can be written as follows The logarithm of the equation is Given a different value of , a linear trial-and-error method is used to solve the equation. The value of with the highest correlation coefficient in the linear regression is the correct value. According to the slope and intercept of the straight line, G R and can be solved.
In the process of equation solution, an initial value of is difficult to determine, therefore G R and are not easy to be calculated. However, according to equation (8), (10) and (11) of the Duong model, the corresponding parameters can be solved quickly. The solution procedures are as follows. (i) A representative straight line with the highest correlation coefficient square R 2 in the log-log curve of q/G p ∼ t is selected to regress the slope and intercept, then the parameters of a and m are determined. (ii) It is known from equation (8) that the curve between q and t −m e a 1−m (t 1−m −1) is a straight line passing through the origin. Initial gas production rate q 1 can be obtained according to the slope. (iii) According to the economic limit production rate of shale gas and the corresponding time, the recoverable reserve of shale gas can be calculated from equation (10).

Establishment of production decline equation
The term 2 is a distribution derived from the normal distribution, whose probability density function is where Γ(n/2) is the Gamma function, whose expression is x is a distribution variable whose range is from 0 to ∞, and n is a parameter that controls the function distribution shape. The probability density curve of 2 distribution is shown in figure 2. It shows that the curve shape is influenced by the value of n. When n ≤ 2, the curve monotonously decreases. Particularly when n = 2, the curve is exponential. When n > 2, the curve is a single-peak distribution that first increases and then decreases. A decline curve of shale gas production without a peak can be analysed using the 2 distribution model with a value of n ≤ 2. Compared with the Weibull model, the 2 distribution model can better describe the production decline law of shale gas with a rapid decline in the early stage, which will be verified in the following example analysis.
Assuming that = n/2 and = 2, then In the interval of x from 0 to ∞, the integral of 2 distribution function equals 1. That is ∫ ∞ 0 f (x)dx = 1. Similar to the Weibull distribution, the production rate calculation formula of shale gas based on 2 distribution is Assuming Assuming a = G R c b+1 Γ(b + 1) , then Equation (22) is a generalised Weng's prediction model (Weng 1991). When b equals 0, it is an exponential decline model without a peak. Zuo et al (2016) proposed the FDC to forecast the production rate of shale gas.
where m is an coefficient corresponding to well index; E (− t ) is the Mittag-Leffler function proposed by Mittag-Leffler (1903) and and are fitting parameters. After a certain time, equation (23) can be simplified to Equation (24) can be derived through mathematical method, which is omitted here. It can also be seen intuitively in figure 3, where m = 1, and are from the fitting of actual data. In a case where later data are the same, the early data of equation (24) are bigger than equation (23). If the early data decline too quickly, equation (23) cannot fit accurately.
Actually, equation (24) can be written as a form of the power function equation where a ′ = m Γ(1− ) , b ′ = − . When the value of c in equation (21) or equation (22) is big enough, equation (21) or equation (22) can approach equation (25). It can be included that the 2 distribution model covers the characteristics of the FDC in equation (23), which is actually the power function equation in a later stage.

Solution of production decline equation
Getting the logarithm in both ends of equation (22), then Assuming a 0 = ln a, a 1 = b , a 2 = −1∕c, y = ln q , x 1 = ln t, x 2 = t, then Parameters a 0 , a 1 , a 2 can be obtained through binary linear regression of equation (27). Then a, b and c can be calculated as follows.
The solution process of the 2 distribution model is easier than for the Weibull distribution model, whose three parameters cannot be obtained at once through a linear regression method.   In summary, the Arps model can only describe the production decline process, while the Weibull (Duong) and 2 distribution models can not only describe the production decline process, but also the entire life cycle including production increase and decline processes. In the Arps model, exponential decline rate is constant, while the decline rate of hyperbolic decline and harmonic decline gradually decreases. When the initial decline rate is the same, the exponential decline is the fastest. The second is hyperbolic decline and harmonic decline is the slowest. The Arps model is widely used in a reservoir with boundary flow regime, while the Duong model is suitable for a linear flow regime. Although the Duong model can be derived from the Weibull model, the Duong model is more adaptable than the Weibull model, in which negative value of G R may appear. Because of the complexity of porous media and seepage mechanisms in shale gas reservoirs, there may be deviation in the fitting of production data with the Arps and Weibull (Duong) models. The 2 distribution model is essentially the product of the power function and exponential function of time.
It is suitable to describe the production characteristics of long-tail distribution with a quick early decline and gradually slow decline in the later stage, which will be demonstrated in following examples.

Field case analysis
The production data from three wells of the Fayetteville shale gas reservoir in the study by Zuo et al. (2016) are analysed with production decline formulas of the Arps model, Weibull (or Duong) model and 2 distribution model. The production time of each shale gas well has been more than 5 years, which shows an obvious production decline law.
In the process of data fitting with the Arps model, a different decline index n from 0 to 1 with steps of 0.001 are set to regress the equation. The Arps decline type can be determined by the value of n with a maximum R 2 (square  Analysis results are shown in Table 1, figures 4-6. It is indicated that the overall fitting accuracy of the Arps model, Weibull (Duong) model and 2 model are good, whose average relative error and square of correlation coefficient R 2 are 1.73% and 0.9397, respectively. It is known that Arps model is suitable for boundary dominated flow. It is also possible that the flow is not affected by the boundary, but the fitting of production data, as well 4 and well 5 are good only because of the wide adaptability of the Arps model. The Duong model is suitable for the fitting of a linear or bilinear flow regime, which widely exists in fractures of shale gas reservoirs. For well 2, the R 2 of the Weibull (Duong) model is the highest. The curve of the 2 decline model has the characteristics of long-tail distribution with a quick early decline and slow late decline. Therefore, for well 2 whose early production data declines slowly, the fitting error of the 2 distribution model in the early stage is relatively big, but the fitting in the late stage is good. It is also concluded that the Arps model cannot be completely abandoned in production decline analysis of a shale gas reservoir, but other models can be used as supplements.

Decline analysis of production data from numerical simulations
Because of many interference factors, such as accidental well shutdown and frequent adjustments of working systems in the actual production process of shale gas reservoirs, the true production decline law may be concealed. Therefore, production data from numerical simulations of fractured horizontal wells in shale gas reservoirs are further analysed. A shale gas reservoir numerical simulation model with dimensions of 1600 (length) × 800 (width) × 80 m (thickness) is established. Reservoir porosity is 2.5%. Permeability values are set to be 0.00025, 0.0025 and 0.025 mD in different cases. Initial reservoir pressure and temperature are 30 MPa and 60°C. The minimum bottom-hole pressure constraint during production is 5.0 MPa. Langmuir volume and Langmuir pressure are 2.91 m 3 /t and 5.15 MPa. Reserves of adsorption gas and free gas are 6.57 × 10 8 m 3 and 6.43 × 10 8 m 3 , with a total of 13.0 × 10 8 m 3 . The length of the horizontal well is 1200 m. There are 10 hydraulic fractures with fracture spacings of 120 m, fracture half-length of 100 m and fracture width of 0.005 m. Initial shale gas production rate is 10 × 10 4 m 3 /d and production time is 30 years. Pressure distributions of shale gas reservoirs with different permeabilities are shown in figures 7-9. It is shown that the pressure conformance area of a shale gas reservoir increases with the increase of permeability under the same production time. When permeability is 0.00025 mD, the pressure conformance area is small and it is obviously influenced by fractures. As with the increase of permeability, pressure conformance area enlarges and the influence of fractures on pressure distribution decreases. When permeability is 0.025 mD, the pressure conformance area is elliptical, where the influence of a horizontal well on gas drainage radius is dominant but the influence of fractures is subtle.  Shale gas production curves under different matrix permeabilities are shown in figure 10. It shows that the stable production duration of shale gas increases with an increase in permeability. When permeability is 0.00025 mD, the gas production curve has obvious long-tail distribution characteristics, where early decline rates are large and then decline rates decrease over time because the gradual desorption of shale gas in the process of pressure decreases. When permeability reaches 0.025 mD, the characteristics of long-tail distribution are not obvious and total decline rate is obviously bigger than that of 0.00025 mD.
Production decline analysis results of numerical simulation data are shown in Table 2, figures 11-13. The main characteristic of simulation data is a long-tail distribution whose early data decline quickly and late data decline slowly. On the whole, fitting results of the Weibull (Duong) model and 2 distribution model are better than the Arps model whose deviation of early data is large. In figures 11-13, the Arps model is in harmonic decline, which can better describe a quick early decline and slow decline in the later stage than exponential and hyperbolic decline. If the early production declines too quickly, it is difficult to fit the data with the Arps model. In this case, the Weibull (Duong) and 2 distribution model can be considered. For the case in figure 7 with a permeability of 0.00025 mD, the influence of fractures and the character-istics of long-tail distribution are more obvious. As a result, fitting accuracy of Arps model is the lowest (R 2 = 0.9552). Fitting accuracy of a 2 distribution model (R 2 = 0.9988) is better than the Weibull (Duong) model (R 2 = 0.9942). As for the increase in permeability, the fitting accuracy of the Weibull (Duong) model gradually becomes better than the 2 distribution model. Although the Arps model can also describe the decreasing decline rate over time in a shale gas reservoir, sometimes the results are not good because they does not satisfy boundary dominated flow. Linear flow is a main flow mode of a fractured horizontal well in a shale gas reservoir. If the fracture spacing is large enough, an early radial flow will appear after early linear flow perpendicular to the fractures. When pressure disturbance fronts between adjacent fractures converge, the fracture boundary flow appears. After that, compound linear flow from the unfractured area to fractured area occurs. Therefore, production decline analysis of Weibull (Duong) model may generate some errors because the flow is not pure linear flow. Comparative analysis results demonstrated that 2 distribution models are more suitable for fitting production data with long-tail distribution characteristics than the Arps model and Weibull (Duong) model, which has also been verified in previous field data analysis. 10