# Performance Analysis of Gas Wells Based on the Conventional Decline Parameters and the Flow Integral Equation

Authors: Lixia Zhang; Yong Li; Xinmin Song; Mingxian Wang; Yang Yu; Yingxu He; Zeqi Zhao
Journal: Geofluids
 (2021)

---

## Abstract

The estimation of reserves and performance prediction are two vital tasks for the development of gas reservoirs where the evaluation of gas in place or well-controlled reserves, as the foundation of the performance analysis of gas wells, turns to be exceedingly significant. Advanced production data analysis or modern rate transient analysis (RTA) methods mainly depend on the iterative calculations of material balance quasitime (tca) and type curve fitting, the essence of which is to update the average reservoir pressure data time and again. The traditional Arps’ decline models are of empirical nature despite the convenience and applicability to the constant bottomhole pressure (BHP) condition. In order to avoid the implicit iteration, this paper develops an explicit method for estimating the average reservoir pressure on the basis of dynamic material balance equation (DMBE), termed “flow integral method,” which can be applied to various gas production systems under boundary-dominated flow (BDF). Based on the flow integral method and the decline parameter evaluation, we employ the hyperbolic decline model to model the gas well performance at a constant BHP. The analytical formulations of decline rate and decline exponent are deduced from the DMBE and the static material balance equation (SMBE) considering the elastic compressibilities of rock pore and bound water. The resulting decline parameter method for explicit estimation of gas reserves boasts a solid and rigorous theory foundation that production rate, decline rate, and average reservoir pressure profiles have reference to each other, and its implementation steps are explained in the paper. The SMBE can, combined with the estimated pressure profile by the flow integral method, also be used to determine gas reserves which is not limited to the constant-BHP condition and can calibrate the estimates of the decline parameter method. The proposed methods are proven effective and reliable with several numerical cases at different BHPs and a field example.


---

## Body

## 1. Introduction

The production performance analysis of gas wells runs through the gas reservoir development. The variations of production rate, bottomhole pressure (BHP), cumulative gas production, and average reservoir pressure with time are the basic data for reservoir engineers to implement the analysis and prediction of gas well performance, among which measurements of average reservoir pressure, nonetheless, are far from easy. Therefore, the static material balance relationship (g vs. pave) is generally used to correlate the average reservoir pressure pave with the cumulative gas production Gp by the use of reserves G. The determination of average reservoir pressure, from a certain perspective, is equivalent to the determination of gas reserves. The ways of estimating gas in place or gas reserves [1–3] mainly include analogy, volumetric method, material balance analysis, pressure transient analysis (PTA), rate transient analysis (RTA), history matching, and numerical simulation. The volumetric method [4–6] is often used in geological modeling of gas reservoirs or gas fields, and the estimated reserves by static data have relatively low accuracy. The traditional material balance method [7–9] is widely used, but it requires a certain number of measurement points of formation pressure. Nowadays, PTA and RTA have found their extensive applications in determining the well-controlled reserves. Well testing must be carried out to obtain the accurate pressure data requisite for PTA [10–12] at the expense of production. The daily production data necessary for RTA [13–15], however, are convenient to glean and are analyzed by rigorous theoretical models for which these cost-effective methods have developed rapidly in recent years.The research on RTA can be traced back to the statistical law of rate data of oil and gas wells presented by Arps [16] who empirically summarized the rate-time relationship into three decline types, namely, exponential, hyperbolic, and harmonic. The conventional hyperbolic decline model (HDM) as a classical tool to analyze the production data of constant-BHP cases is easy to apply, but its theoretical foundation is not fully reckoned with.Rate decline analysis came into widespread use during the early 1980s when Fetkovich [17, 18] proved that dimensionless rate qdD against time tdD for decline curve analysis during late BDF is subject to the exponential decline law under the constant-BHP condition. Fetkovich’s type curves, composed of the dimensionless curves generated by the analytical solution in the transient flow period and Arps’ empirical curves in the BDF stage, have created a new direction of modern decline curve analysis where the analytical solution to unsteady flow model and the log-log plot as used in well test analysis are presented. Fetkovich’s approach removes the empiricism of Arps’ models, but its application is still confined to the slightly compressible fluid because it reckons without the change in gas properties.Following Fetkovich’s investigation, Carter [19, 20] explained the variations of gas viscosity μ and compressibility Cg with pressure by introducing the product ratio λλ=μiCgi/μCg and proposed a new set of type curves for gas reservoir analysis. The approximate drawdown parameter takes into account the effect of initial pressure drawdown on the product of gas viscosity and compressibility which enables Carter’s approach to be applied to constant-BHP gas flow systems.In order to address the influence of the changes in production scenarios on data analysis, Blasingame et al. [21] introduced an idea to analyze the production data for variable rate/variable pressure drop systems by using the time function of constant pressure analog which can be obtained by the constant rate analog (or material balance pseudotime tca) for variable-rate flow under BDF conditions. This approach, for gas well analysis, entails the iterative calculations of tca to determine the parameters mbdf and b as shown in Equation (4). After that, Palacio and Blasingame [22] demonstrated that the relation of dimensionless rate qdD vs. dimensionless time tdD for the decline curve analysis during BDF coincides with Arps’ harmonic decline model based on the material balance time (or pseudotime) function, and then, a new kind of type curves was developed. The innovation of Blasingame’s approach is that it not only considers the changes in fluid properties but also fully explains the variations of production schedules as long as the boundary effect can dominate the fluid flow. In addition, other advanced production decline analysis methods such as the normalized pressure integral method (NPI) [23], Agarwal and Gardner’s type curves [24, 25] are also applied to the production data analysis of gas wells. The differences between these approaches and Blasingame’s type curves (the most representative RTA method) mainly lie in the plot functions; hence, they, in essence, are all consistent with the BDF solutions [26, 27] whether to place emphasis on pressure analysis or rate analysis.The type curve analysis approaches need an iterative procedure for the determination of gas reserves or material balance pseudotime before curve fitting and parameter interpretation. In addition to type curve analysis, some researchers also proposed other RTA methods to conduct reserve estimation and performance prediction without curve fitting. Blasingame and Lee [27], for instance, presented an iterative method for the determination of gas in place, termed variable-rate reservoir limit testing (VRRLT) of gas wells, where SMBE and iterations of average reservoir pressure are combined to estimate reserves of the volumetric gas reservoir. Zhang et al. [28] developed the VRRLT for abnormally pressured gas reservoirs by considering the compressibilities of rock pore and bound water in the mathematical model of gas flow and SMBE and proposed the selection strategy of initial iteration value.Mattar and McNeil [29, 30] put forward the “flowing” material balance procedure only applicable to constant-rate cases. This approach, however, may not apply to gas reservoirs because the variation of gas property parameters (such as viscosity, deviation factor, and compressibility) with pressure is not considered and its theoretical basis is not rigorously demonstrated, too.Mattar et al. [31, 32] developed the dynamic material balance (DMB) procedure, an extension of the flowing material balance, to cope with the rate fluctuations of the gas well. This approach is derived from the “stabilized flow” solution of the gas flow model and the material balance equation of the volumetric gas reservoir. Mattar’s dynamic material balance procedure boasts a much stricter fundamental basis than the previous flowing material balance by the use of pseudovariables to tackle the changes in production rate and gas properties. The dynamic material balance equation (DMBE), afterwards, was strictly deduced by Zhang et al. [33] who extended the DMB to abnormally pressured gas reservoirs by incorporating the compressibility effects of irreducible water and rock into SMBE and the gas flow model. Similarly, SMBE, DMBE, and iterations of material balance pseudotime or average reservoir pressure can be combined to estimate gas reserves where the selection strategy of the initial value can refer to Zhang et al. [28].The VRRLT and DMB function well without type curve fitting, but they still entail the iterations of material balance pseudotime or average reservoir pressure for determination of reserves in an implicit fashion. In order to avoid the iterative procedure, some researchers began to attempt more direct and explicit approaches for reserve estimation with the unremitting application. Ye and Ayala [34] proposed a method to linearize the gas diffusion equation by using density function (ρ) and pseudotime factor (i.e., depletion-driven time rescaling factor β). The density function, in reality, ignores the viscosity term in pseudopressure, while the treatment of pseudotime factor is consistent with the traditional definition of pseudotime. They evade the iteration of gas reserves for the first time and use the coincidence of liquid solution and numerical simulation results (or gas well data) in the early stage when changes in gas properties are not obvious to select the dimensionless reservoir radius in liquid-type curves, and then, the differences during BDF are employed to determine the λ and β profiles. An approximate relation of λ to ρ is developed to construct the material balance plot, namely, λ1/B vs. Gp, for estimation of reserves where the slope of log-log plot of μ·Cg vs. ρ is approximately reckoned a constant −B. This approach is targeted at the rate-time analysis of the gas well in the volumetric gas reservoir under a constant-BHP specification.Stumpf and Ayala [35] also presented an explicit method for determining gas in place under constant-BHP conditions by the use of q1−n vs. Gp∗ straight line based on the hyperbolic window (hyperbolic decline period) identification. They stated that the hyperbolic window only exists in the early period of BDF when the decline exponent n can be considered as a constant determined by the initial pseudopressure drawdown ratio rmi and the integral mean value of α with respect to pseudopressure mp so the rate data conform to the hyperbolic decline law. The hyperbolic window can be determined by adjusting the rate and time multipliers via the hyperbolic type curve fitting. This approach, however, seems little easy in terms of the whole analysis process and inapplicable to geopressured gas reservoirs where the elastic influences of irreducible water and rock are not negligible.Alom et al. [36] proposed a type curve fitting method to analyze production data during BDF by using the modified pseudotime function defined by the cumulative gas production Gp. This approach, similar to Blasingame’s type curves in principle, uses polynomials to match the relation of μiCgi/μCg to Gp, but in fact, the relationship between μiCgi/μCg and time t is hard to determine because the profile of average reservoir pressure pave and the relation of Gp vs. pave are unknown.Zhang et al. [37] developed the material balance-quasipressure approximation condition method for determining gas reserves which combines the material balance condition derived from the material-balance principle with the quasipressure approximation condition obtained by solving the mathematical model of gas flow through porous media. This approach also circumvents the determination of material balance pseudotime and works well for both gas reservoirs with abnormal pressure and those with normal pressure under BDF conditions.Recently, Wang and Ayala [38] extended the hyperbolic decline model under constant-BHP conditions presented by Stumpf and Ayala [35] to the special variable-BHP condition and deduced the decline exponent n for the volumetric gas reservoir. Following the analysis steps of Stumpf and Ayala, they also employed the type-curve fitting and q1−n~Gp∗ straight-line relation to implement estimation of reserves. It is worth noting that the ratio of pseudopressure at BHP to that at average reservoir pressure is postulated constant. Under this condition, the hyperbolic decline model can be applied to all data during BDF rather than only the early BDF. These approaches, however, reckoned without the elastic effects of irreducible water and rock. Furthermore, the decline exponent in the hyperbolic decline model (Stumpf and Ayala; Wang and Ayala) is evaluated in terms of variable α on average, that is, n is deemed the integral mean of α with respect to pseudopressure mp from BHP to initial pressure, but in reality, α is a function of average reservoir pressure (larger than BHP).Jongkittinarukorn et al. [39] presented a correlation between production rate, decline rate, and decline exponent that relates q vs. Dn·dlnλ/dpD−1 curve to reserve estimation, based on the viscosity-compressibility ratio λ defined by Carter [20] and the stabilized gas flow equation proposed by Ansah et al. [40, 41], where pD, a dimensionless pressure, is defined as p/Z/pi/Zi; D denotes decline rate and n is decline exponent. They considered the changes in decline exponent for the volumetric gas reservoir at constant-BHP conditions. The synthetic data for validation, nevertheless, is generated by Ansah et al.’s model which is the rationale of Jongkittinarukorn’s approach in the first place.To sum up, the current rate transient analysis methods can be roughly divided into five categories: (1) traditional Arps’ decline curve analysis; (2) Fetkovich’s and Carter’s type curves for constant-BHP production data analysis; (3) modern decline curve analysis represented by Blasingame’s type curves capable of analyzing the production system under variable-rate/variable-BHP conditions; (4) the variable-rate reservoir limits testing [27, 28] and dynamic material balance method [31–33] based on the iterations of material balance pseudotime tca without the curve fitting; (5) some explicit methods [34, 35] to avoid pseudotime calculations. The first method is simple to use, but the theoretical basis is not rigorous. The second is limited to constant-BHP conditions, and curve fitting is complicated. The third and fourth methods are effective in dealing with the BDF data under various production conditions such as constant rate, constant pressure, and even variable rate/variable pressure, but the repeated iterations on average reservoir pressure or pseudotime are indispensable. The fifth category of methods can accomplish reserve estimation in an iteration-free way which has gradually become a research hotspot in the field of gas reservoir performance analysis; however, their applicability to complicated production schedules and abnormally pressured reservoirs remains to be investigated.Therefore, taking into consideration the compressibility effects of rock and irreducible water, this article is aimed at exploring the explicit methodologies for determining the average reservoir pressure and reserves of a gas reservoir, presenting the theoretical proof of the applicability of the traditional hyperbolic decline model to the gas well with a constant-BHP specification, and eliminating the empiricism of HDM and the defect that it fails to estimate gas in place. These methodologies can evade the iteration calculations oftca or pave and lay the foundation for the performance analysis of gas wells.The overall layout of this paper is as follows. Firstly, the variable-flowrate solutions (or correlations between rate and pressure during BDF) are introduced. Secondly, the static material balance equation is introduced to deduce the analytic expressions of decline rate and decline exponent; then, we explain the rationales behind the traditional hyperbolic model to characterize the production decline data of constant-BHP gas well and the fundamentals of the decline parameter method available for reserve estimation. Thirdly, the flow integral method is presented for estimating the average reservoir pressure by a reference point during BDF based on the flow integral equation derived from the DMBE, and the implementation steps of the decline parameter method are described. Then, several case investigations from both simulated and realistic data are presented to demonstrate and validate the proposed methodologies. Finally, the innovations of the article are outlined and the conclusions are drawn.


## 2. Model Development

### 2.1. Correlations of Rate vs. Pressure during BDF

To linearize the diffusivity equation of gas flow, it is usually necessary to introduce pseudovariables such as(1)pp=pi+μiZipi∫pipξμξZξdξ,(2)ta=μiCti∫0t1μCtdt,(3)Ct=eCϕp−piCϕ+1−SwcCg+SwcCw,wherepp (Pa) denotes pseudopressure, pi (Pa) is initial reservoir pressure, μ (Pa·s) represents natural gas viscosity, μi (Pa·s) denotes gas viscosity at pi, Z is gas deviation factor or compressibility factor, Zi denotes gas compressibility factor at pi, t (s) is time, ta (s) denotes pseudotime, Ct (Pa-1) represents system compressibility, Cti (Pa-1) represents system compressibility at pi, Swc represents irreducible water saturation, Cg (Pa-1) denotes gas compressibility, Cϕ (Pa-1) is rock compressibility, and Cw (Pa-1) is water compressibility. All the symbols in this paper are explained in international standard units (SI) for convenience.By solving the mathematical model of gas flow with variable rates for a gas well in the center of a closed reservoir, one obtains the solution of boundary-dominated flow (BDF) represented by(4)ppi−ppwfq=Δppi−wfq=mbdf⋅tca+b,(5)mbdf=2παtαp⋅BgiAhϕiCti=1G⋅2παtαp⋅1−SwciCti,(6)b=μiBgiαpKh4reD4lnreD−3reD4+4reD2−14reD2−12+∑i=1mqi−qi−1qm∑n=1∞2λn2⋅fta,λn,(7)fta,λn=e−λn2·αtK/ϕiμiCtirw2ta−ta,i−1J12reDλnJ12reDλn−J12λn,(8)tca=μiCtiqt∫0tqtμpaveCtpavedt,where ppi or pppi (Pa) denotes the pseudopressure at pi, ppwf (Pa) denotes the pseudopressure at bottomhole pressure (BHP) pwf, and q (m3/s) is the surface production rate of the gas well; tca (s) represents material balance pseudotime, mbdf (Pa/m3) denotes the slope of the Δppi−wf/q vs. tca straight line, and b (Pa·s/m3) is its intercept which is generally reckoned a constant during BDF; G (m3) is gas in place or gas reserves; αp=2π and αt=1 when the international system of units is used.The relation of pseudopressure at the average reservoir pressure to pseudopressure at BHP under boundary-dominated flow is subject to the dynamic material balance equation [33] given by
(9)ppave=ppwf+q⋅b,where pave (Pa) is the average reservoir pressure, and ppave (Pa) denotes the pseudopressure at pave.According to the definition of pseudopressure, Equation (1), we can rewrite Equation (9) as
(10)q=ppave−ppwfb=1b⋅μiZipi∫pipaveξμξZξdξ−∫pipwfξμξZξdξ.Equation (10) has the same variable “b” as Equation (4).


### 2.2. Conventional Decline Parameters

#### 2.2.1. Decline Rate

Taking the derivative of Equation (10) with respect to time, we have
(11)dqdt=1b⋅μiZipipaveμpaveZpave⋅dpavedt−pwfμpwfZpwf⋅dpwfdt.When a gas well operates at a constant BHP,dpwf/dt is equal to 0; then, Equation (11) is transformed into
(12)dqdt=1b⋅μiZipi⋅paveμpaveZpave⋅dpavedt.Taking into account the compressibilities of rock pore and bound water, the static material balance equation (SMBE) of gas reservoirs captures the intrinsic relationship between average formation pressure and cumulative gas production, which is given by(13)gpave=piZi1−GpG,(14)gp=pZp⋅eCϕp−pi−Swcie−Cwp−pi1−Swci,where gp denotes the unified pressure function in the SMBE, Pa; Gp is the cumulative gas production, m3; and Swci represents the irreducible water saturation at pi.Differentiating Equation (13) with respect to t gives
(15)dgpavedt=−piZiGdGpdt=−piZiG⋅q.Differentiating Equation (14) with respect to p yields
(16)dgpdp=Cp1−Swci⋅pZp,(17)Cp=eCϕp−piCgp+Cϕ+Swcie−Cwp−piCw−Cgp,where Cp denotes the compressibility function dependent on Cw, Cϕ, and Cg. Though Cp defined in Equation (17) is different from Ctp defined in Equation (3) in form, they are identical in essence and both can be transformed into 1−Swc·Cg when Cw and Cϕ approximate to 0.According to the chain rule, the derivative of the average reservoir pressurepave with respect to time t can be expressed as
(18)dpavedt=dpavedgpave⋅dgpavedt=dgpdpp=pave−1⋅dgpavedt.By substituting Equations (15) and (16) into Equation (18), one obtains
(19)dpavedt=1−SwciCpave⋅Zpavepave⋅−piZiG⋅q.The decline rateD is defined as
(20)D=−1qdqdt.Substituting Equations (12) and (19) into Equation (20) gives
(21)D=μi1−Swcib⋅G⋅1μpaveCpave.Equation (21) shows that the decline rate D dependent upon μC is a function of pave and decreases with time because 1/μC continues to fall with the falling formation pressure.


#### 2.2.2. Decline Exponent

The decline exponentn is defined as
(22)n=dD−1dt.Substituting Equation (21) into Equation (22) yields
(23)n=bGμi1−Swci⋅dμpaveCpavedpave⋅dpavedt.By putting Equation (19) into Equation (23), we obtain
(24)n=−piμiZi⋅qb⋅dμpaveCpavedpppavedpppavedpave⋅1Cpave⋅Zpavepave.From Equation (1), it follows
(25)dpppavedpave=μiZipi⋅paveμpaveZpave.Substituting Equations (9) and (25) into Equation (24) gives
(26)n=−dμpaveCpavedpppave⋅ppave−ppwfμpaveCpave=dlnμ⋅C−1dlnpp⋅pp−ppwfppp=pave.The decline exponentn in Equation (26) is similar to the formulations given by Chen and Teufel [42, 43] and the difference lies in the new compressibility function Cp in the paper. Equation (26) can be further reduced to
(27)n=spave⋅γpave,(28)sp=ddlnppln1μ⋅C,(29)γp=pp−ppwfpp,where sp denotes the derivative of the logarithm of 1/μ·C with respect to the logarithm of pp, and γp represents the instantaneous pseudopressure level.It is found that the decline exponent is the product ofspave and γpave, all of which are dependent on pave which decreases with time. With falling pave, γpave falls; however, spave tends to rise as a whole. Therefore, n changes slowly with time during the depletion of formation energy and ultimately approaches 0 when gas production is subject to the exponential decline. An equivalent constant is employed to replace the dynamic decline exponent given the slow change in n, and the hyperbolic decline model is used to match the production data of the gas well at a constant BHP in the paper.Note thatb in Equation (9) is considered a constant in the above derivation process and it indeed holds true for a quite long time during BDF. In fact, the value of b, nevertheless, also varies slightly with time.


#### 2.2.3. Estimation of Gas Reserves

The conventional hyperbolic decline model developed by Arps is delineated by three parameters (qi, Di, and n), and the production rate and decline rate are represented by
(30)q=qi⋅1+nDi⋅t−1/n,(31)Dt=Di1+n̄Di⋅t,where qi denotes the initial gas flowrate during BDF, m3/s; Di is the initial decline rate, s-1; and n̄ represents the equivalent decline exponent.From Equation (21), it follows
(32)Dj=μi1−Swcib⋅G⋅1μpjCpj.According to Equation (10), it follows
(33)qj=pppj−ppwfb,where pj is the average formation pressure at a certain time tj, Pa; qj denotes the production rate at tj, m3/s; and Dj is the decline rate at tj, 1/s.Based on Equations (32) and (33), qj/Dj can be written as
(34)Djqj=μi1−SwciG⋅1μpjCpj⋅1pppj−ppwf.From Equation (34), it follows
(35)Gj=qjDj⋅μi1−SwciμpjCpj⋅1pppj−ppwf,where Gj denotes the calculated reserves from the production data (rate, decline rate, and average reservoir pressure) at tj, m3.Equation (35) reveals that gas reserves can be estimated according to the intrinsic relationship between production rate, decline rate, and average reservoir pressure, but it needs to know the relation of formation pressure to time. Measurements of average reservoir pressure are expensive and time-consuming so it is unrealistic to glean or record the whole profile of pave. The lack of pave data prompts us to infer the unknown data points from a limited number of measurement points of average reservoir pressure. What follows in the passage will introduce a method to determine the pave profile merely based on a certain measurement point during BDF.


### 2.3. The Flow Integral Method for Determination of Average Reservoir Pressure

From Equation (10), it follows
(36)qj=1b⋅μiZipi∫pwf,jpjξμξZξdξ.If the formation pressure attk is known, then the production rate at that time qk can be expressed as
(37)qk=1b⋅μiZipi∫pwf,kpkξμξZξdξ.Combining Equation (36) with Equation (37), we obtain
(38)∫pwf,jpjξμξZξdξ=qjqk∫pwf,kpkξμξZξdξ,where pk is the average reservoir pressure at tj, Pa; pwf,k denotes the bottomhole pressure at tj, Pa.Equation (38), termed “flow integral equation,” indicates that the average reservoir pressure at tj (i.e., pj) can be solved by the integral operation of p/μZ with respect to p. To speed up the solution process, we utilize the 6th-degree polynomial to approximate the nonlinear function p/μZ, that is,
(39)fintp=pμpZp≈a6p6+a5p5+a4p4+a3p3+a2p2+a1p+a0.From Equation (39), it follows
(40)Fintp=∫ξμξZξdξ≈a67p7+a56p6+a45p5+a34p4+a23p3+a12p2+a0p+c,where Fintp is the antiderivative of fintp with respect to p, c is an arbitrary constant (for convenience, let c=0 in the paper), and ai (i=0,1,2⋯6) denotes the coefficient of the polynomial function used to approximate p/μZ.Substituting Equation (40) into Equation (38), the rearranged equation is represented by
(41)Fintpj=qjqkFintpk−Fintpwf,k+Fintpwf,j.Equation (41) is applicable to all cases under BDF, not limited to the constant-BHP condition. If the BHP remains unchanged, then this equation collapses to
(42)Fintpj=qjqkFintpk−Fintpwf+Fintpwf.The right-side terms of Equation (41) are known and its left side is a 7th-degree polynomial; thus, pj can be easily solved in the range of pwf,j<pj<pi. The explicit method for estimating the average reservoir pressure data, based on the flow integral equation and one measure point during BDF, is called the “flow integral method” in the paper.


### 2.4. The Decline Parameter Method for Determination of Gas Reserves

The profile of average reservoir pressure can be determined by the flow integral method so that the terms related topave in Equation (35) can also be determined. The profiles of q, D, and pave can, therefore, be used to conduct the explicit estimation of gas reserves only based on production data. The specific steps are as follows:
(1)
Sort out the gas reservoir and fluid property parameters, and determine the multinomial coefficientsa0~a6 in fintp(2)
Use the HDM to fit the rate-time data of gas well at a constant BHP and determine three model parameters (qi, Di, and n̄)(3)
Use Equation (31) to determine the decline rate D(4)
Glean the availablepave data and select a certain data point at tk during BDF to calculate the values of average reservoir pressure with time by the flow integral method(5)
Useqj, Dj, and pj to determine Gj based on Equation (35). The average value of Gj is the estimated reserves GHD given by(43)GHD=1N∑j=1NGj,where N denotes the number of data points.The above iteration-free method, termed the “decline parameter method,” can explicitly calculate gas reserves of a constant-BHP system based on the hyperbolic decline parameters and the flow integral equation. The static material balance relationship can be also employed to estimate gas in place after the profile ofpave is obtained by the flow integral method. The intercept of gpave vs. Gp straight line is pi/Zi, and its slope −m is related to gas reserves Gsmb defined as
(44)Gsmb=piZi⋅1m.The value ofGsmb determined by the SMBE is comparable to GHD estimated by the decline parameter method.


## 2.1. Correlations of Rate vs. Pressure during BDF

To linearize the diffusivity equation of gas flow, it is usually necessary to introduce pseudovariables such as(1)pp=pi+μiZipi∫pipξμξZξdξ,(2)ta=μiCti∫0t1μCtdt,(3)Ct=eCϕp−piCϕ+1−SwcCg+SwcCw,wherepp (Pa) denotes pseudopressure, pi (Pa) is initial reservoir pressure, μ (Pa·s) represents natural gas viscosity, μi (Pa·s) denotes gas viscosity at pi, Z is gas deviation factor or compressibility factor, Zi denotes gas compressibility factor at pi, t (s) is time, ta (s) denotes pseudotime, Ct (Pa-1) represents system compressibility, Cti (Pa-1) represents system compressibility at pi, Swc represents irreducible water saturation, Cg (Pa-1) denotes gas compressibility, Cϕ (Pa-1) is rock compressibility, and Cw (Pa-1) is water compressibility. All the symbols in this paper are explained in international standard units (SI) for convenience.By solving the mathematical model of gas flow with variable rates for a gas well in the center of a closed reservoir, one obtains the solution of boundary-dominated flow (BDF) represented by(4)ppi−ppwfq=Δppi−wfq=mbdf⋅tca+b,(5)mbdf=2παtαp⋅BgiAhϕiCti=1G⋅2παtαp⋅1−SwciCti,(6)b=μiBgiαpKh4reD4lnreD−3reD4+4reD2−14reD2−12+∑i=1mqi−qi−1qm∑n=1∞2λn2⋅fta,λn,(7)fta,λn=e−λn2·αtK/ϕiμiCtirw2ta−ta,i−1J12reDλnJ12reDλn−J12λn,(8)tca=μiCtiqt∫0tqtμpaveCtpavedt,where ppi or pppi (Pa) denotes the pseudopressure at pi, ppwf (Pa) denotes the pseudopressure at bottomhole pressure (BHP) pwf, and q (m3/s) is the surface production rate of the gas well; tca (s) represents material balance pseudotime, mbdf (Pa/m3) denotes the slope of the Δppi−wf/q vs. tca straight line, and b (Pa·s/m3) is its intercept which is generally reckoned a constant during BDF; G (m3) is gas in place or gas reserves; αp=2π and αt=1 when the international system of units is used.The relation of pseudopressure at the average reservoir pressure to pseudopressure at BHP under boundary-dominated flow is subject to the dynamic material balance equation [33] given by
(9)ppave=ppwf+q⋅b,where pave (Pa) is the average reservoir pressure, and ppave (Pa) denotes the pseudopressure at pave.According to the definition of pseudopressure, Equation (1), we can rewrite Equation (9) as
(10)q=ppave−ppwfb=1b⋅μiZipi∫pipaveξμξZξdξ−∫pipwfξμξZξdξ.Equation (10) has the same variable “b” as Equation (4).


## 2.2. Conventional Decline Parameters

### 2.2.1. Decline Rate

Taking the derivative of Equation (10) with respect to time, we have
(11)dqdt=1b⋅μiZipipaveμpaveZpave⋅dpavedt−pwfμpwfZpwf⋅dpwfdt.When a gas well operates at a constant BHP,dpwf/dt is equal to 0; then, Equation (11) is transformed into
(12)dqdt=1b⋅μiZipi⋅paveμpaveZpave⋅dpavedt.Taking into account the compressibilities of rock pore and bound water, the static material balance equation (SMBE) of gas reservoirs captures the intrinsic relationship between average formation pressure and cumulative gas production, which is given by(13)gpave=piZi1−GpG,(14)gp=pZp⋅eCϕp−pi−Swcie−Cwp−pi1−Swci,where gp denotes the unified pressure function in the SMBE, Pa; Gp is the cumulative gas production, m3; and Swci represents the irreducible water saturation at pi.Differentiating Equation (13) with respect to t gives
(15)dgpavedt=−piZiGdGpdt=−piZiG⋅q.Differentiating Equation (14) with respect to p yields
(16)dgpdp=Cp1−Swci⋅pZp,(17)Cp=eCϕp−piCgp+Cϕ+Swcie−Cwp−piCw−Cgp,where Cp denotes the compressibility function dependent on Cw, Cϕ, and Cg. Though Cp defined in Equation (17) is different from Ctp defined in Equation (3) in form, they are identical in essence and both can be transformed into 1−Swc·Cg when Cw and Cϕ approximate to 0.According to the chain rule, the derivative of the average reservoir pressurepave with respect to time t can be expressed as
(18)dpavedt=dpavedgpave⋅dgpavedt=dgpdpp=pave−1⋅dgpavedt.By substituting Equations (15) and (16) into Equation (18), one obtains
(19)dpavedt=1−SwciCpave⋅Zpavepave⋅−piZiG⋅q.The decline rateD is defined as
(20)D=−1qdqdt.Substituting Equations (12) and (19) into Equation (20) gives
(21)D=μi1−Swcib⋅G⋅1μpaveCpave.Equation (21) shows that the decline rate D dependent upon μC is a function of pave and decreases with time because 1/μC continues to fall with the falling formation pressure.


### 2.2.2. Decline Exponent

The decline exponentn is defined as
(22)n=dD−1dt.Substituting Equation (21) into Equation (22) yields
(23)n=bGμi1−Swci⋅dμpaveCpavedpave⋅dpavedt.By putting Equation (19) into Equation (23), we obtain
(24)n=−piμiZi⋅qb⋅dμpaveCpavedpppavedpppavedpave⋅1Cpave⋅Zpavepave.From Equation (1), it follows
(25)dpppavedpave=μiZipi⋅paveμpaveZpave.Substituting Equations (9) and (25) into Equation (24) gives
(26)n=−dμpaveCpavedpppave⋅ppave−ppwfμpaveCpave=dlnμ⋅C−1dlnpp⋅pp−ppwfppp=pave.The decline exponentn in Equation (26) is similar to the formulations given by Chen and Teufel [42, 43] and the difference lies in the new compressibility function Cp in the paper. Equation (26) can be further reduced to
(27)n=spave⋅γpave,(28)sp=ddlnppln1μ⋅C,(29)γp=pp−ppwfpp,where sp denotes the derivative of the logarithm of 1/μ·C with respect to the logarithm of pp, and γp represents the instantaneous pseudopressure level.It is found that the decline exponent is the product ofspave and γpave, all of which are dependent on pave which decreases with time. With falling pave, γpave falls; however, spave tends to rise as a whole. Therefore, n changes slowly with time during the depletion of formation energy and ultimately approaches 0 when gas production is subject to the exponential decline. An equivalent constant is employed to replace the dynamic decline exponent given the slow change in n, and the hyperbolic decline model is used to match the production data of the gas well at a constant BHP in the paper.Note thatb in Equation (9) is considered a constant in the above derivation process and it indeed holds true for a quite long time during BDF. In fact, the value of b, nevertheless, also varies slightly with time.


### 2.2.3. Estimation of Gas Reserves

The conventional hyperbolic decline model developed by Arps is delineated by three parameters (qi, Di, and n), and the production rate and decline rate are represented by
(30)q=qi⋅1+nDi⋅t−1/n,(31)Dt=Di1+n̄Di⋅t,where qi denotes the initial gas flowrate during BDF, m3/s; Di is the initial decline rate, s-1; and n̄ represents the equivalent decline exponent.From Equation (21), it follows
(32)Dj=μi1−Swcib⋅G⋅1μpjCpj.According to Equation (10), it follows
(33)qj=pppj−ppwfb,where pj is the average formation pressure at a certain time tj, Pa; qj denotes the production rate at tj, m3/s; and Dj is the decline rate at tj, 1/s.Based on Equations (32) and (33), qj/Dj can be written as
(34)Djqj=μi1−SwciG⋅1μpjCpj⋅1pppj−ppwf.From Equation (34), it follows
(35)Gj=qjDj⋅μi1−SwciμpjCpj⋅1pppj−ppwf,where Gj denotes the calculated reserves from the production data (rate, decline rate, and average reservoir pressure) at tj, m3.Equation (35) reveals that gas reserves can be estimated according to the intrinsic relationship between production rate, decline rate, and average reservoir pressure, but it needs to know the relation of formation pressure to time. Measurements of average reservoir pressure are expensive and time-consuming so it is unrealistic to glean or record the whole profile of pave. The lack of pave data prompts us to infer the unknown data points from a limited number of measurement points of average reservoir pressure. What follows in the passage will introduce a method to determine the pave profile merely based on a certain measurement point during BDF.


## 2.2.1. Decline Rate

Taking the derivative of Equation (10) with respect to time, we have
(11)dqdt=1b⋅μiZipipaveμpaveZpave⋅dpavedt−pwfμpwfZpwf⋅dpwfdt.When a gas well operates at a constant BHP,dpwf/dt is equal to 0; then, Equation (11) is transformed into
(12)dqdt=1b⋅μiZipi⋅paveμpaveZpave⋅dpavedt.Taking into account the compressibilities of rock pore and bound water, the static material balance equation (SMBE) of gas reservoirs captures the intrinsic relationship between average formation pressure and cumulative gas production, which is given by(13)gpave=piZi1−GpG,(14)gp=pZp⋅eCϕp−pi−Swcie−Cwp−pi1−Swci,where gp denotes the unified pressure function in the SMBE, Pa; Gp is the cumulative gas production, m3; and Swci represents the irreducible water saturation at pi.Differentiating Equation (13) with respect to t gives
(15)dgpavedt=−piZiGdGpdt=−piZiG⋅q.Differentiating Equation (14) with respect to p yields
(16)dgpdp=Cp1−Swci⋅pZp,(17)Cp=eCϕp−piCgp+Cϕ+Swcie−Cwp−piCw−Cgp,where Cp denotes the compressibility function dependent on Cw, Cϕ, and Cg. Though Cp defined in Equation (17) is different from Ctp defined in Equation (3) in form, they are identical in essence and both can be transformed into 1−Swc·Cg when Cw and Cϕ approximate to 0.According to the chain rule, the derivative of the average reservoir pressurepave with respect to time t can be expressed as
(18)dpavedt=dpavedgpave⋅dgpavedt=dgpdpp=pave−1⋅dgpavedt.By substituting Equations (15) and (16) into Equation (18), one obtains
(19)dpavedt=1−SwciCpave⋅Zpavepave⋅−piZiG⋅q.The decline rateD is defined as
(20)D=−1qdqdt.Substituting Equations (12) and (19) into Equation (20) gives
(21)D=μi1−Swcib⋅G⋅1μpaveCpave.Equation (21) shows that the decline rate D dependent upon μC is a function of pave and decreases with time because 1/μC continues to fall with the falling formation pressure.


## 2.2.2. Decline Exponent

The decline exponentn is defined as
(22)n=dD−1dt.Substituting Equation (21) into Equation (22) yields
(23)n=bGμi1−Swci⋅dμpaveCpavedpave⋅dpavedt.By putting Equation (19) into Equation (23), we obtain
(24)n=−piμiZi⋅qb⋅dμpaveCpavedpppavedpppavedpave⋅1Cpave⋅Zpavepave.From Equation (1), it follows
(25)dpppavedpave=μiZipi⋅paveμpaveZpave.Substituting Equations (9) and (25) into Equation (24) gives
(26)n=−dμpaveCpavedpppave⋅ppave−ppwfμpaveCpave=dlnμ⋅C−1dlnpp⋅pp−ppwfppp=pave.The decline exponentn in Equation (26) is similar to the formulations given by Chen and Teufel [42, 43] and the difference lies in the new compressibility function Cp in the paper. Equation (26) can be further reduced to
(27)n=spave⋅γpave,(28)sp=ddlnppln1μ⋅C,(29)γp=pp−ppwfpp,where sp denotes the derivative of the logarithm of 1/μ·C with respect to the logarithm of pp, and γp represents the instantaneous pseudopressure level.It is found that the decline exponent is the product ofspave and γpave, all of which are dependent on pave which decreases with time. With falling pave, γpave falls; however, spave tends to rise as a whole. Therefore, n changes slowly with time during the depletion of formation energy and ultimately approaches 0 when gas production is subject to the exponential decline. An equivalent constant is employed to replace the dynamic decline exponent given the slow change in n, and the hyperbolic decline model is used to match the production data of the gas well at a constant BHP in the paper.Note thatb in Equation (9) is considered a constant in the above derivation process and it indeed holds true for a quite long time during BDF. In fact, the value of b, nevertheless, also varies slightly with time.


## 2.2.3. Estimation of Gas Reserves

The conventional hyperbolic decline model developed by Arps is delineated by three parameters (qi, Di, and n), and the production rate and decline rate are represented by
(30)q=qi⋅1+nDi⋅t−1/n,(31)Dt=Di1+n̄Di⋅t,where qi denotes the initial gas flowrate during BDF, m3/s; Di is the initial decline rate, s-1; and n̄ represents the equivalent decline exponent.From Equation (21), it follows
(32)Dj=μi1−Swcib⋅G⋅1μpjCpj.According to Equation (10), it follows
(33)qj=pppj−ppwfb,where pj is the average formation pressure at a certain time tj, Pa; qj denotes the production rate at tj, m3/s; and Dj is the decline rate at tj, 1/s.Based on Equations (32) and (33), qj/Dj can be written as
(34)Djqj=μi1−SwciG⋅1μpjCpj⋅1pppj−ppwf.From Equation (34), it follows
(35)Gj=qjDj⋅μi1−SwciμpjCpj⋅1pppj−ppwf,where Gj denotes the calculated reserves from the production data (rate, decline rate, and average reservoir pressure) at tj, m3.Equation (35) reveals that gas reserves can be estimated according to the intrinsic relationship between production rate, decline rate, and average reservoir pressure, but it needs to know the relation of formation pressure to time. Measurements of average reservoir pressure are expensive and time-consuming so it is unrealistic to glean or record the whole profile of pave. The lack of pave data prompts us to infer the unknown data points from a limited number of measurement points of average reservoir pressure. What follows in the passage will introduce a method to determine the pave profile merely based on a certain measurement point during BDF.


## 2.3. The Flow Integral Method for Determination of Average Reservoir Pressure

From Equation (10), it follows
(36)qj=1b⋅μiZipi∫pwf,jpjξμξZξdξ.If the formation pressure attk is known, then the production rate at that time qk can be expressed as
(37)qk=1b⋅μiZipi∫pwf,kpkξμξZξdξ.Combining Equation (36) with Equation (37), we obtain
(38)∫pwf,jpjξμξZξdξ=qjqk∫pwf,kpkξμξZξdξ,where pk is the average reservoir pressure at tj, Pa; pwf,k denotes the bottomhole pressure at tj, Pa.Equation (38), termed “flow integral equation,” indicates that the average reservoir pressure at tj (i.e., pj) can be solved by the integral operation of p/μZ with respect to p. To speed up the solution process, we utilize the 6th-degree polynomial to approximate the nonlinear function p/μZ, that is,
(39)fintp=pμpZp≈a6p6+a5p5+a4p4+a3p3+a2p2+a1p+a0.From Equation (39), it follows
(40)Fintp=∫ξμξZξdξ≈a67p7+a56p6+a45p5+a34p4+a23p3+a12p2+a0p+c,where Fintp is the antiderivative of fintp with respect to p, c is an arbitrary constant (for convenience, let c=0 in the paper), and ai (i=0,1,2⋯6) denotes the coefficient of the polynomial function used to approximate p/μZ.Substituting Equation (40) into Equation (38), the rearranged equation is represented by
(41)Fintpj=qjqkFintpk−Fintpwf,k+Fintpwf,j.Equation (41) is applicable to all cases under BDF, not limited to the constant-BHP condition. If the BHP remains unchanged, then this equation collapses to
(42)Fintpj=qjqkFintpk−Fintpwf+Fintpwf.The right-side terms of Equation (41) are known and its left side is a 7th-degree polynomial; thus, pj can be easily solved in the range of pwf,j<pj<pi. The explicit method for estimating the average reservoir pressure data, based on the flow integral equation and one measure point during BDF, is called the “flow integral method” in the paper.


## 2.4. The Decline Parameter Method for Determination of Gas Reserves

The profile of average reservoir pressure can be determined by the flow integral method so that the terms related topave in Equation (35) can also be determined. The profiles of q, D, and pave can, therefore, be used to conduct the explicit estimation of gas reserves only based on production data. The specific steps are as follows:
(1)
Sort out the gas reservoir and fluid property parameters, and determine the multinomial coefficientsa0~a6 in fintp(2)
Use the HDM to fit the rate-time data of gas well at a constant BHP and determine three model parameters (qi, Di, and n̄)(3)
Use Equation (31) to determine the decline rate D(4)
Glean the availablepave data and select a certain data point at tk during BDF to calculate the values of average reservoir pressure with time by the flow integral method(5)
Useqj, Dj, and pj to determine Gj based on Equation (35). The average value of Gj is the estimated reserves GHD given by(43)GHD=1N∑j=1NGj,where N denotes the number of data points.The above iteration-free method, termed the “decline parameter method,” can explicitly calculate gas reserves of a constant-BHP system based on the hyperbolic decline parameters and the flow integral equation. The static material balance relationship can be also employed to estimate gas in place after the profile ofpave is obtained by the flow integral method. The intercept of gpave vs. Gp straight line is pi/Zi, and its slope −m is related to gas reserves Gsmb defined as
(44)Gsmb=piZi⋅1m.The value ofGsmb determined by the SMBE is comparable to GHD estimated by the decline parameter method.


## 3. Model Validation

In this section, several synthetic cases with different BHPs and a field case are used to illustrate the effectiveness and applicability of the proposed methodologies above. Furthermore, the modified Stumpf method is also introduced for comparison.


### 3.1. Synthetic Cases: Gas Reservoir with Normal Pressure

A commercial numerical simulator generates the production data of a gas well at constant BHPs centered in a circular reservoir. The simulated gas reservoir has an original pressure (pi) of 25 MPa, a temperature (Ti) of 70°C (or 343.15 K), and a reference depth of 2500 m. The molar mass of natural gas (Mg) is 18 g/mol, and the pseudocritical temperature and pressure determined by Sutton’s correlations [44, 45] are 201.34 K and 4.60 MPa, respectively. Hall and Yarborough’s approach [46] is used for calculations of Z-factor and gas compressibility, and Londono et al.’s correlations [47, 48] for gas viscosity. Other property parameters are shown in Table 1.Table 1
Reservoir and gas properties for numerical cases with normal pressure.
Property parametersProperty valuesProperty parametersProperty valuesGrid300×120×1h12 mϕi0.18rw0.1 mKr9 mDρsc0.750146 kg/m3Kθ9 mDMg18 g/molKz0.9 mDZsc0.997515Swci0.22Zi0.896867dr1.3 mCϕ5.462×10−4MPa−1dθ3°Cw4.178×10−4MPa−1dz12 mμw0.432089 cppi25 MPaμi2.193343823×10−2cpTi343.15 Kμsc1.075008625×10−2cpTsc293.15 KCgi3.152713092×10−2MPa−1psc0.101325 MPaCti2.522927812×10−2MPa−1Tpc201.34 KBgi4.265593097×10−3m3/m3ppc4.60 MPaVpi1 032 656 m3re390 mG188 829 866 m3Figure1 illustrates the relationships of p/μ·Z vs. p and μ·Z vs. p for the simulated gas reservoir with normal pressure, and their nonlinearities may result in large error when real pressure or pressure-squared approach is adopted. Equation (45), a 6th-degree polynomial as represented by the dotted line in Figure 1, is used to approximate the relation of p/μ·Z to p, and the polynomial boasts high precision within the range of 0.2 to 25 MPa.
(45)pμZ=−3.854104×10−6p6+2.969799×10−4p5−6.396372×10−3p4−8.938681×10−3p3−1.885397×10−1p2+8.271143×101p−1.192612×10−1.Figure 1
Relation curves ofp/μZ and μZ vs. pressure (p) for the normal pressure cases.The curve ofgp vs. p shown in Figure 2 reveals the static material balance relationship, and 1/μ·C vs. p in Figure 3 reflects the change in the decline rate to some extent. The maximum value of the instantaneous decline exponent sp dependent on p is displayed in Figure 4, and the relationship between the pseudopressure pp and p is shown in Figure 5.Figure 2
The relation curve ofgp vs. p for the normal pressure cases.Figure 3
The relation curve of1/μ·C vs. p for the normal pressure cases.Figure 4
The relation curve ofsp vs. p for the normal pressure cases.Figure 5
The relation curve ofppp vs. p for the normal pressure cases.There are 2064 time steps set in the numerical simulation, and the gas well is produced under various constant-BHP conditions for 7000 days (about 19.16 years). Four simulated production schedules correspond to four synthetic cases as shown in Table2.Table 2
Production schedules for the numerical cases with normal pressure at constant BHPs.
t (d)Time stepΔt (h)Δt (d)pwf (MPa)1~21~1210.041667Constant:Case 1-1: 20Case 1-2: 15Case 1-3: 10Case 1-4: 051330.1251490.37515~16120.53~30017~314241301-1100315~7144821101-2000715~10147232001~30001015~12649643001~70001265~20641205


#### 3.1.1. Case 1-1

The decline parameter method is employed to estimate the gas reserves for Case 1-1 with a BHP of 20 MPa. The solution steps are as follows:Step 1. Determination of property parameters.Table1, Figures 1–5, and Equation (45) display the corresponding results.Step 2. Determination of hyperbolic decline model.Equation (30) is used to match the rate-time data of the first 6000 days, and three model parameters are determined as follows: qi=1.165088×104m3/d, Di=3.898168×10−4d−1, and n̄=0.133561. The production decline law characterized by the hyperbolic decline model almost coincides with the actual production rate profile as shown in Figure 6 where an equivalent decline exponent constant could scarcely give rise to large error in the production rate. The instantaneous decline exponent n would approximate to 0 only when the rate q is very low and close to 0 though it varies slowly with time.Figure 6
The real production rate profile and the predicted profile by the HDM for Case 1-1.Step 3. Calculation of decline rate.Equation (31) is employed to estimate the corresponding decline rate Dt at each time point t, as shown in Figure 7.Figure 7
The relation curve ofD vs. t for Case 1-1.Step 4. Estimation of average reservoir pressure.The flow integral method is utilized to determine the average reservoir pressure profile. A reference point is necessary for that method. Taketk=1000d as an example here. Equation (42) is adopted to determine the average reservoir pressure at each time point; then, the predicted pave profile is shown in Figure 8.Figure 8
The estimated profile of average reservoir pressure (pave) for Case 1-1.Step 5. Estimation of gas reserves.The value ofGj in Equation (35) can be calculated by the decline rate in Step 3 and the average reservoir pressure in Step 4. Figure 9 shows the profile of estimated reserves, and the resulting GHD is 1.907956937×108m3 with an error of 1.041%.Figure 9
The estimated profile of gas reserves (Gj) for Case 1-1.Step 6. Verification by the static material balance equation.Thepave profile in Step 4 can also be used to draw gpave vs. Gp curve as shown in Figure 10 where pi/Zi=27.874801234MPa and −m=−1.480821926×10−3 MPa/(104 m3). The value of Gsmb calculated by Equation (44) is 1.882387122×108m3 with an error of -0.313%. The small difference between Gsmb and GHD validates the decline parameter method.Figure 10
The static material balance relation ofgpave to Gp for Case 1-1.


#### 3.1.2. Case 1-2

As previously mentioned, Equation (30) is used to match the 6000 d rate profile for Case 2-1 with a constant BHP of 15 MPa, and three model parameters are determined as follows: qi=2.326642×104m3/d, Di=4.040669×10−4d−1, and n̄=0.254215. The production rate predicted by the hyperbolic decline model is consistent with the real production profile as shown in Figure 11.Figure 11
The realq profile and the predicted profile by the hyperbolic decline model for Case 1-2.Similarly, take thepave and q at 1000 d as the reference data. Figure 12 displays the estimated reserve profile, and the resulting value of GHD is 1.906270758×108m3 with an error of 0.952%. In addition, the straight line of gpave vs. Gp has a negative slope −m of −1.490280824×10−3 MPa/(104 m3) and Gsmb determined by the SMBE is 1.870439503×108m3 with an error of -0.946%.Figure 12
The estimated profile ofGj for Case 1-2.


#### 3.1.3. Case 1-3

Figure13 shows the real and predicted profiles of production rate with a BHP of 10 MPa where the dashed line is obtained by the hyperbolic decline model with qi=3.486867×104m3/d, Di=4.217685×10−4d−1, and n̄=0.325809. The ultimate GHD of 1.897272945×108m3 determined by the decline parameter method, as shown in Figure 14, compares favourably with Gsmb of 1.859010973×108m3 determined by the SMBE.Figure 13
The realq profile and the predicted profile by the hyperbolic decline model for Case 1-3.Figure 14
The estimated profile ofGj for Case 1-3.


#### 3.1.4. Case 1-4

The real and predicted profiles of production rate for Case 1-4 with a BHP of 5 MPa are exhibited in Figure15 where the dotted line derives from the hyperbolic decline model with qi=4.645448×104m3/d, Di=4.402959×10−4d−1, and n̄=0.331420. Figure 16 shows the estimated profile of Gj. Similar to the previous cases, GHD of 1.915021414×108m3 coincides with Gsmb of 1.867302158973 × 108 m3.Figure 15
The realq profile and the predicted profile by the hyperbolic decline model for Case 1-4.Figure 16
The estimated profile ofGj for Case 1-4.Table3(a) presents the calculation results of the decline parameter method for the synthetic gas reservoir with normal pressure where Error1 denotes the calculation error of GHD obtained by the decline parameter method and Error2 represents the corresponding error of Gsmb determined by the static material balance equation.Table 3(a)
The analysis results of the decline parameter method for numerical cases with normal pressure
Normal pressure casesCase 1-1Case 1-2Case 1-3Case 1-4pwf (MPa)2015105qi (104 m3/d)1.1650882.3266423.4868674.645448Di (10-4 d-1)3.8981684.0406694.2176854.402959n̄0.1335610.2542150.3258090.331420GHD (108 m3)1.9079569371.9062707581.8972729451.915 021414Error1 (%)1.0410.9520.4751.415m (10-3 MPa/104 m3)1.4808219261.4902808241.4994425341.492784717Gsmb (108 m3)1.8823871221.8704395031.8590109731.867302158Error2 (%)-0.313-0.946-1.551-1.112(b)
The analysis results of the modified Stumpf method for numerical cases with normal pressure
Normal pressure casesCase 1-1Case 1-2Case 1-3Case 1-4pwf (MPa)2015105qi (104 m3/d)1.1629162.3199573.4892274.705304Di (10-4 d-1)3.8529203.9681734.2282814.715305na0.1156820.2274450.3252680.408158Gm (108 m3)1.9216768671.9560241991.9946865972.042141878Error3 (%)1.7683.5875.6348.147For the sake of contrast, the modified Stumpf method is developed in the appendix following the research of Stumpf and Ayala [35]. It is necessary for this method to evaluate the average decline exponent (na) on the basis of Equations (28) and (29); then, the straight line of q1−na vs. Gp can be drawn. Its slope and intercept can be used to estimate the decline parameters such as qi, Di, and G (see Appendix). The corresponding analysis results for numerical cases with normal pressure are illustrated in Table 3(b) and Figures 33 and 34 (see Appendix), where Gm denotes the value of gas reserves calculated by the modified Stumpf method and Error3 represents its error of estimation.It can be seen from Tables3(a) and 3(b) that the decline parameter method provides desirable estimations of gas reserves and decline parameters though the results of the modified Stumpf method are also acceptable except Case 1-4. The differences between two methods lie in the calculation approach of decline exponent and the processing of production data. The former seems more effective in the performance analysis of the relatively low-pressure gas reservoir than the latter which is mainly applicable to early BDF period.


### 3.2. Synthetic Cases: Gas Reservoir with High Pressure

The gas reservoir with a high initial pressure of 80 MPa is simulated at the temperature of 120°C. Table4 presents the reservoir and gas properties for the high-pressure cases. Four constant-BHP schedules are generated as shown in Table 5.Table 4
Reservoir and gas properties for numerical cases with high pressure.
Property parametersProperty valuesProperty parametersProperty valuesGrid300×120×1h12 mϕi0.15rw0.1 mKr5 mDρsc0.750146 kg/m3Kθ5 mDMg18 g/molKz0.5 mDZsc0.997515Swci0.22Zi1.518784dr1.3 mCϕ1.929×10−3MPa−1dθ3°Cw4.134×10−4MPa−1dz12 mμw0.235499 cppi80 MPaμi3.705514139×10−2 cpTi393.15 Kμsc1.075008625×10−2cpTsc293.15 KCgi5.088279959×10−3MPa−1psc0.101325 MPaCti5.988806368×10−3MPa−1Tpc201.34 KBgi2.586255933×10−3m3/m3ppc4.60 MPaVpi860 546 m3re390 mG259 535 855 m3Table 5
Production schedules for the numerical cases with high pressure at constant BHPs.
t (d)TSΔt (h)Δt (d)pwf (MPa)1~21~1210.041667Constant:Case 2-1: 64Case 2-2: 52Case 2-3: 40Case 2-4: 241330.1251490.37515~16120.53~30017~314241301-1100315~7144821101-2000715~10147232001~30001015~12649643001~70001265~20641205Figures17–21 exhibit the relation curves of p/μ·Z and μ·Z vs. p, gp vs. p, 1/μ·C vs. p, sp vs. p, and pp vs. p for the synthetic cases with high pressure, respectively. The nonlinear relationship of p/μ·Z vs. p within the range of 0.2≤p≤80 MPa in Figure 17 can approximately be represented by the following polynomial:
(46)pμZ=7.003288×10−8p6−1.815134×10−5p5+1.719386×10−3p4−6.336605×10−2p3−2.343059×10−1p2+7.586930×101p−4.957952.Figure 17
Relation curves ofp/μZ and μZ vs. pressure (p) for the high-pressure cases.Figure 18
The relation curve ofgp vs. p for the high-pressure cases.Figure 19
The relation curve of1/μ·C vs. p for the high-pressure cases.Figure 20
The relation curve ofsp vs. p for the high-pressure cases.Figure 21
The relation curve ofppp vs. p for the high-pressure cases.


#### 3.2.1. Case 2-1

Figure22 displays the simulated production profile and the predicted profile by the HDM with three parameters: qi=2.027258×104m3/d, Di=6.193329×10−4d−1, and n̄=0.030078 for Case 2-1 with a BHP of 64 MPa. The estimated profile of Gj is shown in Figure 23, and GHD is 2.631125906×108m3. The value of Gsmb estimated by the SMBE is 2.596071194×108m3 which compares well with the result of the decline parameter method.Figure 22
The realq profile and the predicted profile by the HDM for Case 2-1.Figure 23
The estimated profile ofGj for Case 2-1.


#### 3.2.2. Case 2-2

The real and predicted profiles of production rate for Case 2-2 with a BHP of 52 MPa are illustrated in Figure24 where the dashed line is obtained by the hyperbolic decline model with qi=3.548412×104m3/d, Di=6.102143×10−4d−1, and n̄=0.071877. Figure 25 shows the estimated profile of Gj determined by Equation (35). Similarly, GHD of 2.650991436×108m3 coincides with Gsmb of 2.596831903×108m3.Figure 24
The realq profile and the predicted profile by the hyperbolic decline model for Case 2-2.Figure 25
The estimated profile ofGj for Case 2-2.


#### 3.2.3. Case 2-3

Figure26 shows the real and predicted profiles of production rate for Case 2-3 with a BHP of 40 MPa in which the dashed line represents the matched hyperbolic decline model with qi=5.073846×104m3/d, Di=6.059150×10−4d−1, and n̄=0.137588. Then, the ultimate GHD of 2.667464699×108m3 determined by the decline parameter method, as shown in Figure 27, is comparable with Gsmb of 2.596426085×108m3 determined by the gpave vs. Gp curve.Figure 26
The realq profile and the predicted profile by the hyperbolic decline model for Case 2-3.Figure 27
The estimated profile ofGj for Case 2-3.


#### 3.2.4. Case 2-4

As mentioned earlier, Equation (30) is used to match the rate profile of the first 6000 days for Case 2-4 with a constant BHP of 24 MPa, and three model parameters are determined as follows: qi=7.120844×104m3/d, Di=6.082153×10−4d−1, and n̄=0.258816. The production rate predicted by the hyperbolic decline model is in line with the real production profile as shown in Figure 28. The decline parameter method is applied, and Figure 29 displays the estimated reserve profile with the resulting GHD of 2.664141330×108m3 in line with Gsmb of 2.596254338×108m3 determined by the SMBE.Figure 28
The realq profile and the predicted profile by the hyperbolic decline model for Case 2-4.Figure 29
The estimated profile ofGj for Case 2-4.The computation performance of the decline parameter method for the synthetic gas reservoir with high pressure is summarized in Table6(a). The corresponding analysis results of the modified Stumpf method are displayed in Table 6(b) and Figures 35 and 36 (see Appendix). Two methods both give similar estimations of gas reserves, but the latter shows an increasing error trend with a decrease in pwf.Table 6(a)
The analysis results of the decline parameter method for numerical cases with high pressure
High-pressure casesCase 2-1Case 2-2Case 2-3Case 2-4pwf (MPa)64524024qi (104 m3/d)2.0272583.5484125.0738467.120844Di (10-4 d-1)6.1933296.1021436.0591506.082153n̄0.0300780.0718770.1375880.258816GHD (108 m3)2.6311259062.6509914362.6674646992.664141330Error1 (%)1.3782.1442.7782.650m (10-3 MPa/104 m3)2.0289782092.0283838452.0287008792.028835081Gsmb (108 m3)2.5960711942.5968319032.5964260852.596254338Error2 (%)0.0270.0570.0410.035(b)
The analysis results of modified Stumpf method for numerical cases with high pressure
High-pressure casesCase 2-1Case 2-2Case 2-3Case 2-4pwf (MPa)64524024qi (104 m3/d)2.0331033.5639585.0982277.158882Di (10-4 d-1)6.2808796.2297616.1949826.235047na0.0428620.0906990.1572940.281530Gm (108 m3)2.6101385522.6230317362.6394887272.673866274Error3 (%)0.5691.0661.7003.025It is found, from the previous synthetic cases with whether normal pressure or high pressure, that an equivalent decline exponentn̄ is capable of matching the production decline data of the gas well with a constant BHP though the instantaneous decline exponent n as a function of pave changes slowly with time. The numerical cases fully demonstrate that
(1)
The decline parameters of HDM are affected by the BHP of the gas well, and the initial production rate during BDFqi and the equivalent decline exponent n̄ are more sensitive to BHP values than the initial decline rate Di. The lower the bottomhole flowing pressure the gas well is produced at, the greater qi and n̄ are, while the change in Di is little obvious(2)
The flow integral method, an explicit approach for calculations of average reservoir pressure based on the flow integral equation, can estimate thepave profile effectively and conveniently because only one data point during BDF is necessary, and in reality, these estimations are not confined to the constant-BHP situations(3)
The decline parameter method, an explicit technique for calculations of gas in place (or gas reserves) based on the conventional decline parameters in HDM can reliably evaluate reserves of the production system with a constant BHP. Furthermore, the value of gas reserves calculated by the static material balance equation is available for the mutual verification with the estimate of the decline parameter method(4)
The production data matching in the decline parameter method could be superior to the analytic calculations ofqi, Di, and na in the modified Stumpf method based on the straight line of q1−na vs. Gp which is only suitable for early BDF data. The former extends the range of application of the latter. These approaches, however, can be verified against each other(5)
The decline parameter method is valid for both high-pressure and normal pressure gas reservoirs. The calculation errors for the high-pressure system, however, are slightly larger than those of the low-pressure system which may be attributed to the more decrease in formation pressure and a bit severer fluctuation inb during BDF for high-pressure gas reservoirs


### 3.3. Field Case: Gas Well A in West Virginia

Case 3, a field example provided by Fetkovich et al. [18], reflects the production scenario of Gas Well A in West Virginia. Fraim and Wattenbarger [49] also investigated this case and digitized the raw rate-time data as shown in Figure 30. The hyperbolic decline model is used to fit the rate profile, and the resulting three parameters are as follows: qi=6.509603×104m3/d, Di=2.534238×10−3d−1, and n̄=0.907180.Figure 30
The realq profile and the predicted profile by the hyperbolic decline model for Case 3.Since the specific gravity of the natural gas is 0.57, we obtain the pseudocritical temperatureTpc and pressure ppc are 191.22 K and 4.61 MPa, respectively, according to Sutton [44, 45]. Other property parameters are shown in Table 7.Table 7
Reservoir and gas properties for Case 3.
Property parametersProperty valuesProperty parametersProperty valuespi28.785612 MPah21.336 mTi344.26 Krw0.1079 mpwf3.447379 MPaρsc0.687754 kg/m3ϕi0.06Mg16.509480 g/molSwci0.35Zsc0.997913Tsc293.15 KZi0.952170psc0.101325 MPaμi2.209414019×10−2cpTpc191.218772 Kμsc1.100728475×10−2cpppc4.605256 MPaCgi2.589484787×10−2MPa−1Cϕ8.816×10−4MPa−1Cti1.785752111×10−2MPa−1Cw4.122×10−4MPa−1Bgi3.944216913×10−3m3/m3As revealed by the dotted line in Figure31, the relation of p/μ·Z to p for Case 3 can be represented by the following polynomial in the pressure range of 0.2 to 30:
(47)pμZ=−1.612592×10−6p6+1.292460×10−4p5−2.355782×10−3p4−3.709335×10−2p3−9.545551×10−2p2+8.033991×101p+9.824220×10−2.Figure 31
Relation curves ofp/μZ and μZ vs. pressure (p) for Case 3.According to Fetkovich et al., the average reservoir pressure at 200 days is 3268 psi (or 22.532067 MPa), and the corresponding production rate is estimated to be4.289933×104m3/d by the matched hyperbolic model. Then, the formation pressure profile is determined by the flow integral method. Figure 32 displays the estimated reserve profile by implementing the above steps of the decline parameter method.Figure 32
The estimated profile ofGj for Case 3.The average value ofGj for Case 3 leads to a GHD of 0.857228577×108m3; however, Gsmb is not estimated because of unavailable cumulative production. The negligible difference between GHD in the paper and the estimate for gas in place of 0.303451×107 Mscf (or 0.859277542×108m3) predicted by Fraim and Wattenbarger [49] using normalized time matching demonstrates the applicability and feasibility of the decline parameter method.


## 3.1. Synthetic Cases: Gas Reservoir with Normal Pressure

A commercial numerical simulator generates the production data of a gas well at constant BHPs centered in a circular reservoir. The simulated gas reservoir has an original pressure (pi) of 25 MPa, a temperature (Ti) of 70°C (or 343.15 K), and a reference depth of 2500 m. The molar mass of natural gas (Mg) is 18 g/mol, and the pseudocritical temperature and pressure determined by Sutton’s correlations [44, 45] are 201.34 K and 4.60 MPa, respectively. Hall and Yarborough’s approach [46] is used for calculations of Z-factor and gas compressibility, and Londono et al.’s correlations [47, 48] for gas viscosity. Other property parameters are shown in Table 1.Table 1
Reservoir and gas properties for numerical cases with normal pressure.
Property parametersProperty valuesProperty parametersProperty valuesGrid300×120×1h12 mϕi0.18rw0.1 mKr9 mDρsc0.750146 kg/m3Kθ9 mDMg18 g/molKz0.9 mDZsc0.997515Swci0.22Zi0.896867dr1.3 mCϕ5.462×10−4MPa−1dθ3°Cw4.178×10−4MPa−1dz12 mμw0.432089 cppi25 MPaμi2.193343823×10−2cpTi343.15 Kμsc1.075008625×10−2cpTsc293.15 KCgi3.152713092×10−2MPa−1psc0.101325 MPaCti2.522927812×10−2MPa−1Tpc201.34 KBgi4.265593097×10−3m3/m3ppc4.60 MPaVpi1 032 656 m3re390 mG188 829 866 m3Figure1 illustrates the relationships of p/μ·Z vs. p and μ·Z vs. p for the simulated gas reservoir with normal pressure, and their nonlinearities may result in large error when real pressure or pressure-squared approach is adopted. Equation (45), a 6th-degree polynomial as represented by the dotted line in Figure 1, is used to approximate the relation of p/μ·Z to p, and the polynomial boasts high precision within the range of 0.2 to 25 MPa.
(45)pμZ=−3.854104×10−6p6+2.969799×10−4p5−6.396372×10−3p4−8.938681×10−3p3−1.885397×10−1p2+8.271143×101p−1.192612×10−1.Figure 1
Relation curves ofp/μZ and μZ vs. pressure (p) for the normal pressure cases.The curve ofgp vs. p shown in Figure 2 reveals the static material balance relationship, and 1/μ·C vs. p in Figure 3 reflects the change in the decline rate to some extent. The maximum value of the instantaneous decline exponent sp dependent on p is displayed in Figure 4, and the relationship between the pseudopressure pp and p is shown in Figure 5.Figure 2
The relation curve ofgp vs. p for the normal pressure cases.Figure 3
The relation curve of1/μ·C vs. p for the normal pressure cases.Figure 4
The relation curve ofsp vs. p for the normal pressure cases.Figure 5
The relation curve ofppp vs. p for the normal pressure cases.There are 2064 time steps set in the numerical simulation, and the gas well is produced under various constant-BHP conditions for 7000 days (about 19.16 years). Four simulated production schedules correspond to four synthetic cases as shown in Table2.Table 2
Production schedules for the numerical cases with normal pressure at constant BHPs.
t (d)Time stepΔt (h)Δt (d)pwf (MPa)1~21~1210.041667Constant:Case 1-1: 20Case 1-2: 15Case 1-3: 10Case 1-4: 051330.1251490.37515~16120.53~30017~314241301-1100315~7144821101-2000715~10147232001~30001015~12649643001~70001265~20641205


### 3.1.1. Case 1-1

The decline parameter method is employed to estimate the gas reserves for Case 1-1 with a BHP of 20 MPa. The solution steps are as follows:Step 1. Determination of property parameters.Table1, Figures 1–5, and Equation (45) display the corresponding results.Step 2. Determination of hyperbolic decline model.Equation (30) is used to match the rate-time data of the first 6000 days, and three model parameters are determined as follows: qi=1.165088×104m3/d, Di=3.898168×10−4d−1, and n̄=0.133561. The production decline law characterized by the hyperbolic decline model almost coincides with the actual production rate profile as shown in Figure 6 where an equivalent decline exponent constant could scarcely give rise to large error in the production rate. The instantaneous decline exponent n would approximate to 0 only when the rate q is very low and close to 0 though it varies slowly with time.Figure 6
The real production rate profile and the predicted profile by the HDM for Case 1-1.Step 3. Calculation of decline rate.Equation (31) is employed to estimate the corresponding decline rate Dt at each time point t, as shown in Figure 7.Figure 7
The relation curve ofD vs. t for Case 1-1.Step 4. Estimation of average reservoir pressure.The flow integral method is utilized to determine the average reservoir pressure profile. A reference point is necessary for that method. Taketk=1000d as an example here. Equation (42) is adopted to determine the average reservoir pressure at each time point; then, the predicted pave profile is shown in Figure 8.Figure 8
The estimated profile of average reservoir pressure (pave) for Case 1-1.Step 5. Estimation of gas reserves.The value ofGj in Equation (35) can be calculated by the decline rate in Step 3 and the average reservoir pressure in Step 4. Figure 9 shows the profile of estimated reserves, and the resulting GHD is 1.907956937×108m3 with an error of 1.041%.Figure 9
The estimated profile of gas reserves (Gj) for Case 1-1.Step 6. Verification by the static material balance equation.Thepave profile in Step 4 can also be used to draw gpave vs. Gp curve as shown in Figure 10 where pi/Zi=27.874801234MPa and −m=−1.480821926×10−3 MPa/(104 m3). The value of Gsmb calculated by Equation (44) is 1.882387122×108m3 with an error of -0.313%. The small difference between Gsmb and GHD validates the decline parameter method.Figure 10
The static material balance relation ofgpave to Gp for Case 1-1.


### 3.1.2. Case 1-2

As previously mentioned, Equation (30) is used to match the 6000 d rate profile for Case 2-1 with a constant BHP of 15 MPa, and three model parameters are determined as follows: qi=2.326642×104m3/d, Di=4.040669×10−4d−1, and n̄=0.254215. The production rate predicted by the hyperbolic decline model is consistent with the real production profile as shown in Figure 11.Figure 11
The realq profile and the predicted profile by the hyperbolic decline model for Case 1-2.Similarly, take thepave and q at 1000 d as the reference data. Figure 12 displays the estimated reserve profile, and the resulting value of GHD is 1.906270758×108m3 with an error of 0.952%. In addition, the straight line of gpave vs. Gp has a negative slope −m of −1.490280824×10−3 MPa/(104 m3) and Gsmb determined by the SMBE is 1.870439503×108m3 with an error of -0.946%.Figure 12
The estimated profile ofGj for Case 1-2.


### 3.1.3. Case 1-3

Figure13 shows the real and predicted profiles of production rate with a BHP of 10 MPa where the dashed line is obtained by the hyperbolic decline model with qi=3.486867×104m3/d, Di=4.217685×10−4d−1, and n̄=0.325809. The ultimate GHD of 1.897272945×108m3 determined by the decline parameter method, as shown in Figure 14, compares favourably with Gsmb of 1.859010973×108m3 determined by the SMBE.Figure 13
The realq profile and the predicted profile by the hyperbolic decline model for Case 1-3.Figure 14
The estimated profile ofGj for Case 1-3.


### 3.1.4. Case 1-4

The real and predicted profiles of production rate for Case 1-4 with a BHP of 5 MPa are exhibited in Figure15 where the dotted line derives from the hyperbolic decline model with qi=4.645448×104m3/d, Di=4.402959×10−4d−1, and n̄=0.331420. Figure 16 shows the estimated profile of Gj. Similar to the previous cases, GHD of 1.915021414×108m3 coincides with Gsmb of 1.867302158973 × 108 m3.Figure 15
The realq profile and the predicted profile by the hyperbolic decline model for Case 1-4.Figure 16
The estimated profile ofGj for Case 1-4.Table3(a) presents the calculation results of the decline parameter method for the synthetic gas reservoir with normal pressure where Error1 denotes the calculation error of GHD obtained by the decline parameter method and Error2 represents the corresponding error of Gsmb determined by the static material balance equation.Table 3(a)
The analysis results of the decline parameter method for numerical cases with normal pressure
Normal pressure casesCase 1-1Case 1-2Case 1-3Case 1-4pwf (MPa)2015105qi (104 m3/d)1.1650882.3266423.4868674.645448Di (10-4 d-1)3.8981684.0406694.2176854.402959n̄0.1335610.2542150.3258090.331420GHD (108 m3)1.9079569371.9062707581.8972729451.915 021414Error1 (%)1.0410.9520.4751.415m (10-3 MPa/104 m3)1.4808219261.4902808241.4994425341.492784717Gsmb (108 m3)1.8823871221.8704395031.8590109731.867302158Error2 (%)-0.313-0.946-1.551-1.112(b)
The analysis results of the modified Stumpf method for numerical cases with normal pressure
Normal pressure casesCase 1-1Case 1-2Case 1-3Case 1-4pwf (MPa)2015105qi (104 m3/d)1.1629162.3199573.4892274.705304Di (10-4 d-1)3.8529203.9681734.2282814.715305na0.1156820.2274450.3252680.408158Gm (108 m3)1.9216768671.9560241991.9946865972.042141878Error3 (%)1.7683.5875.6348.147For the sake of contrast, the modified Stumpf method is developed in the appendix following the research of Stumpf and Ayala [35]. It is necessary for this method to evaluate the average decline exponent (na) on the basis of Equations (28) and (29); then, the straight line of q1−na vs. Gp can be drawn. Its slope and intercept can be used to estimate the decline parameters such as qi, Di, and G (see Appendix). The corresponding analysis results for numerical cases with normal pressure are illustrated in Table 3(b) and Figures 33 and 34 (see Appendix), where Gm denotes the value of gas reserves calculated by the modified Stumpf method and Error3 represents its error of estimation.It can be seen from Tables3(a) and 3(b) that the decline parameter method provides desirable estimations of gas reserves and decline parameters though the results of the modified Stumpf method are also acceptable except Case 1-4. The differences between two methods lie in the calculation approach of decline exponent and the processing of production data. The former seems more effective in the performance analysis of the relatively low-pressure gas reservoir than the latter which is mainly applicable to early BDF period.


## 3.1.1. Case 1-1

The decline parameter method is employed to estimate the gas reserves for Case 1-1 with a BHP of 20 MPa. The solution steps are as follows:Step 1. Determination of property parameters.Table1, Figures 1–5, and Equation (45) display the corresponding results.Step 2. Determination of hyperbolic decline model.Equation (30) is used to match the rate-time data of the first 6000 days, and three model parameters are determined as follows: qi=1.165088×104m3/d, Di=3.898168×10−4d−1, and n̄=0.133561. The production decline law characterized by the hyperbolic decline model almost coincides with the actual production rate profile as shown in Figure 6 where an equivalent decline exponent constant could scarcely give rise to large error in the production rate. The instantaneous decline exponent n would approximate to 0 only when the rate q is very low and close to 0 though it varies slowly with time.Figure 6
The real production rate profile and the predicted profile by the HDM for Case 1-1.Step 3. Calculation of decline rate.Equation (31) is employed to estimate the corresponding decline rate Dt at each time point t, as shown in Figure 7.Figure 7
The relation curve ofD vs. t for Case 1-1.Step 4. Estimation of average reservoir pressure.The flow integral method is utilized to determine the average reservoir pressure profile. A reference point is necessary for that method. Taketk=1000d as an example here. Equation (42) is adopted to determine the average reservoir pressure at each time point; then, the predicted pave profile is shown in Figure 8.Figure 8
The estimated profile of average reservoir pressure (pave) for Case 1-1.Step 5. Estimation of gas reserves.The value ofGj in Equation (35) can be calculated by the decline rate in Step 3 and the average reservoir pressure in Step 4. Figure 9 shows the profile of estimated reserves, and the resulting GHD is 1.907956937×108m3 with an error of 1.041%.Figure 9
The estimated profile of gas reserves (Gj) for Case 1-1.Step 6. Verification by the static material balance equation.Thepave profile in Step 4 can also be used to draw gpave vs. Gp curve as shown in Figure 10 where pi/Zi=27.874801234MPa and −m=−1.480821926×10−3 MPa/(104 m3). The value of Gsmb calculated by Equation (44) is 1.882387122×108m3 with an error of -0.313%. The small difference between Gsmb and GHD validates the decline parameter method.Figure 10
The static material balance relation ofgpave to Gp for Case 1-1.


## 3.1.2. Case 1-2

As previously mentioned, Equation (30) is used to match the 6000 d rate profile for Case 2-1 with a constant BHP of 15 MPa, and three model parameters are determined as follows: qi=2.326642×104m3/d, Di=4.040669×10−4d−1, and n̄=0.254215. The production rate predicted by the hyperbolic decline model is consistent with the real production profile as shown in Figure 11.Figure 11
The realq profile and the predicted profile by the hyperbolic decline model for Case 1-2.Similarly, take thepave and q at 1000 d as the reference data. Figure 12 displays the estimated reserve profile, and the resulting value of GHD is 1.906270758×108m3 with an error of 0.952%. In addition, the straight line of gpave vs. Gp has a negative slope −m of −1.490280824×10−3 MPa/(104 m3) and Gsmb determined by the SMBE is 1.870439503×108m3 with an error of -0.946%.Figure 12
The estimated profile ofGj for Case 1-2.


## 3.1.3. Case 1-3

Figure13 shows the real and predicted profiles of production rate with a BHP of 10 MPa where the dashed line is obtained by the hyperbolic decline model with qi=3.486867×104m3/d, Di=4.217685×10−4d−1, and n̄=0.325809. The ultimate GHD of 1.897272945×108m3 determined by the decline parameter method, as shown in Figure 14, compares favourably with Gsmb of 1.859010973×108m3 determined by the SMBE.Figure 13
The realq profile and the predicted profile by the hyperbolic decline model for Case 1-3.Figure 14
The estimated profile ofGj for Case 1-3.


## 3.1.4. Case 1-4

The real and predicted profiles of production rate for Case 1-4 with a BHP of 5 MPa are exhibited in Figure15 where the dotted line derives from the hyperbolic decline model with qi=4.645448×104m3/d, Di=4.402959×10−4d−1, and n̄=0.331420. Figure 16 shows the estimated profile of Gj. Similar to the previous cases, GHD of 1.915021414×108m3 coincides with Gsmb of 1.867302158973 × 108 m3.Figure 15
The realq profile and the predicted profile by the hyperbolic decline model for Case 1-4.Figure 16
The estimated profile ofGj for Case 1-4.Table3(a) presents the calculation results of the decline parameter method for the synthetic gas reservoir with normal pressure where Error1 denotes the calculation error of GHD obtained by the decline parameter method and Error2 represents the corresponding error of Gsmb determined by the static material balance equation.Table 3(a)
The analysis results of the decline parameter method for numerical cases with normal pressure
Normal pressure casesCase 1-1Case 1-2Case 1-3Case 1-4pwf (MPa)2015105qi (104 m3/d)1.1650882.3266423.4868674.645448Di (10-4 d-1)3.8981684.0406694.2176854.402959n̄0.1335610.2542150.3258090.331420GHD (108 m3)1.9079569371.9062707581.8972729451.915 021414Error1 (%)1.0410.9520.4751.415m (10-3 MPa/104 m3)1.4808219261.4902808241.4994425341.492784717Gsmb (108 m3)1.8823871221.8704395031.8590109731.867302158Error2 (%)-0.313-0.946-1.551-1.112(b)
The analysis results of the modified Stumpf method for numerical cases with normal pressure
Normal pressure casesCase 1-1Case 1-2Case 1-3Case 1-4pwf (MPa)2015105qi (104 m3/d)1.1629162.3199573.4892274.705304Di (10-4 d-1)3.8529203.9681734.2282814.715305na0.1156820.2274450.3252680.408158Gm (108 m3)1.9216768671.9560241991.9946865972.042141878Error3 (%)1.7683.5875.6348.147For the sake of contrast, the modified Stumpf method is developed in the appendix following the research of Stumpf and Ayala [35]. It is necessary for this method to evaluate the average decline exponent (na) on the basis of Equations (28) and (29); then, the straight line of q1−na vs. Gp can be drawn. Its slope and intercept can be used to estimate the decline parameters such as qi, Di, and G (see Appendix). The corresponding analysis results for numerical cases with normal pressure are illustrated in Table 3(b) and Figures 33 and 34 (see Appendix), where Gm denotes the value of gas reserves calculated by the modified Stumpf method and Error3 represents its error of estimation.It can be seen from Tables3(a) and 3(b) that the decline parameter method provides desirable estimations of gas reserves and decline parameters though the results of the modified Stumpf method are also acceptable except Case 1-4. The differences between two methods lie in the calculation approach of decline exponent and the processing of production data. The former seems more effective in the performance analysis of the relatively low-pressure gas reservoir than the latter which is mainly applicable to early BDF period.


## 3.2. Synthetic Cases: Gas Reservoir with High Pressure

The gas reservoir with a high initial pressure of 80 MPa is simulated at the temperature of 120°C. Table4 presents the reservoir and gas properties for the high-pressure cases. Four constant-BHP schedules are generated as shown in Table 5.Table 4
Reservoir and gas properties for numerical cases with high pressure.
Property parametersProperty valuesProperty parametersProperty valuesGrid300×120×1h12 mϕi0.15rw0.1 mKr5 mDρsc0.750146 kg/m3Kθ5 mDMg18 g/molKz0.5 mDZsc0.997515Swci0.22Zi1.518784dr1.3 mCϕ1.929×10−3MPa−1dθ3°Cw4.134×10−4MPa−1dz12 mμw0.235499 cppi80 MPaμi3.705514139×10−2 cpTi393.15 Kμsc1.075008625×10−2cpTsc293.15 KCgi5.088279959×10−3MPa−1psc0.101325 MPaCti5.988806368×10−3MPa−1Tpc201.34 KBgi2.586255933×10−3m3/m3ppc4.60 MPaVpi860 546 m3re390 mG259 535 855 m3Table 5
Production schedules for the numerical cases with high pressure at constant BHPs.
t (d)TSΔt (h)Δt (d)pwf (MPa)1~21~1210.041667Constant:Case 2-1: 64Case 2-2: 52Case 2-3: 40Case 2-4: 241330.1251490.37515~16120.53~30017~314241301-1100315~7144821101-2000715~10147232001~30001015~12649643001~70001265~20641205Figures17–21 exhibit the relation curves of p/μ·Z and μ·Z vs. p, gp vs. p, 1/μ·C vs. p, sp vs. p, and pp vs. p for the synthetic cases with high pressure, respectively. The nonlinear relationship of p/μ·Z vs. p within the range of 0.2≤p≤80 MPa in Figure 17 can approximately be represented by the following polynomial:
(46)pμZ=7.003288×10−8p6−1.815134×10−5p5+1.719386×10−3p4−6.336605×10−2p3−2.343059×10−1p2+7.586930×101p−4.957952.Figure 17
Relation curves ofp/μZ and μZ vs. pressure (p) for the high-pressure cases.Figure 18
The relation curve ofgp vs. p for the high-pressure cases.Figure 19
The relation curve of1/μ·C vs. p for the high-pressure cases.Figure 20
The relation curve ofsp vs. p for the high-pressure cases.Figure 21
The relation curve ofppp vs. p for the high-pressure cases.


### 3.2.1. Case 2-1

Figure22 displays the simulated production profile and the predicted profile by the HDM with three parameters: qi=2.027258×104m3/d, Di=6.193329×10−4d−1, and n̄=0.030078 for Case 2-1 with a BHP of 64 MPa. The estimated profile of Gj is shown in Figure 23, and GHD is 2.631125906×108m3. The value of Gsmb estimated by the SMBE is 2.596071194×108m3 which compares well with the result of the decline parameter method.Figure 22
The realq profile and the predicted profile by the HDM for Case 2-1.Figure 23
The estimated profile ofGj for Case 2-1.


### 3.2.2. Case 2-2

The real and predicted profiles of production rate for Case 2-2 with a BHP of 52 MPa are illustrated in Figure24 where the dashed line is obtained by the hyperbolic decline model with qi=3.548412×104m3/d, Di=6.102143×10−4d−1, and n̄=0.071877. Figure 25 shows the estimated profile of Gj determined by Equation (35). Similarly, GHD of 2.650991436×108m3 coincides with Gsmb of 2.596831903×108m3.Figure 24
The realq profile and the predicted profile by the hyperbolic decline model for Case 2-2.Figure 25
The estimated profile ofGj for Case 2-2.


### 3.2.3. Case 2-3

Figure26 shows the real and predicted profiles of production rate for Case 2-3 with a BHP of 40 MPa in which the dashed line represents the matched hyperbolic decline model with qi=5.073846×104m3/d, Di=6.059150×10−4d−1, and n̄=0.137588. Then, the ultimate GHD of 2.667464699×108m3 determined by the decline parameter method, as shown in Figure 27, is comparable with Gsmb of 2.596426085×108m3 determined by the gpave vs. Gp curve.Figure 26
The realq profile and the predicted profile by the hyperbolic decline model for Case 2-3.Figure 27
The estimated profile ofGj for Case 2-3.


### 3.2.4. Case 2-4

As mentioned earlier, Equation (30) is used to match the rate profile of the first 6000 days for Case 2-4 with a constant BHP of 24 MPa, and three model parameters are determined as follows: qi=7.120844×104m3/d, Di=6.082153×10−4d−1, and n̄=0.258816. The production rate predicted by the hyperbolic decline model is in line with the real production profile as shown in Figure 28. The decline parameter method is applied, and Figure 29 displays the estimated reserve profile with the resulting GHD of 2.664141330×108m3 in line with Gsmb of 2.596254338×108m3 determined by the SMBE.Figure 28
The realq profile and the predicted profile by the hyperbolic decline model for Case 2-4.Figure 29
The estimated profile ofGj for Case 2-4.The computation performance of the decline parameter method for the synthetic gas reservoir with high pressure is summarized in Table6(a). The corresponding analysis results of the modified Stumpf method are displayed in Table 6(b) and Figures 35 and 36 (see Appendix). Two methods both give similar estimations of gas reserves, but the latter shows an increasing error trend with a decrease in pwf.Table 6(a)
The analysis results of the decline parameter method for numerical cases with high pressure
High-pressure casesCase 2-1Case 2-2Case 2-3Case 2-4pwf (MPa)64524024qi (104 m3/d)2.0272583.5484125.0738467.120844Di (10-4 d-1)6.1933296.1021436.0591506.082153n̄0.0300780.0718770.1375880.258816GHD (108 m3)2.6311259062.6509914362.6674646992.664141330Error1 (%)1.3782.1442.7782.650m (10-3 MPa/104 m3)2.0289782092.0283838452.0287008792.028835081Gsmb (108 m3)2.5960711942.5968319032.5964260852.596254338Error2 (%)0.0270.0570.0410.035(b)
The analysis results of modified Stumpf method for numerical cases with high pressure
High-pressure casesCase 2-1Case 2-2Case 2-3Case 2-4pwf (MPa)64524024qi (104 m3/d)2.0331033.5639585.0982277.158882Di (10-4 d-1)6.2808796.2297616.1949826.235047na0.0428620.0906990.1572940.281530Gm (108 m3)2.6101385522.6230317362.6394887272.673866274Error3 (%)0.5691.0661.7003.025It is found, from the previous synthetic cases with whether normal pressure or high pressure, that an equivalent decline exponentn̄ is capable of matching the production decline data of the gas well with a constant BHP though the instantaneous decline exponent n as a function of pave changes slowly with time. The numerical cases fully demonstrate that
(1)
The decline parameters of HDM are affected by the BHP of the gas well, and the initial production rate during BDFqi and the equivalent decline exponent n̄ are more sensitive to BHP values than the initial decline rate Di. The lower the bottomhole flowing pressure the gas well is produced at, the greater qi and n̄ are, while the change in Di is little obvious(2)
The flow integral method, an explicit approach for calculations of average reservoir pressure based on the flow integral equation, can estimate thepave profile effectively and conveniently because only one data point during BDF is necessary, and in reality, these estimations are not confined to the constant-BHP situations(3)
The decline parameter method, an explicit technique for calculations of gas in place (or gas reserves) based on the conventional decline parameters in HDM can reliably evaluate reserves of the production system with a constant BHP. Furthermore, the value of gas reserves calculated by the static material balance equation is available for the mutual verification with the estimate of the decline parameter method(4)
The production data matching in the decline parameter method could be superior to the analytic calculations ofqi, Di, and na in the modified Stumpf method based on the straight line of q1−na vs. Gp which is only suitable for early BDF data. The former extends the range of application of the latter. These approaches, however, can be verified against each other(5)
The decline parameter method is valid for both high-pressure and normal pressure gas reservoirs. The calculation errors for the high-pressure system, however, are slightly larger than those of the low-pressure system which may be attributed to the more decrease in formation pressure and a bit severer fluctuation inb during BDF for high-pressure gas reservoirs


## 3.2.1. Case 2-1

Figure22 displays the simulated production profile and the predicted profile by the HDM with three parameters: qi=2.027258×104m3/d, Di=6.193329×10−4d−1, and n̄=0.030078 for Case 2-1 with a BHP of 64 MPa. The estimated profile of Gj is shown in Figure 23, and GHD is 2.631125906×108m3. The value of Gsmb estimated by the SMBE is 2.596071194×108m3 which compares well with the result of the decline parameter method.Figure 22
The realq profile and the predicted profile by the HDM for Case 2-1.Figure 23
The estimated profile ofGj for Case 2-1.


## 3.2.2. Case 2-2

The real and predicted profiles of production rate for Case 2-2 with a BHP of 52 MPa are illustrated in Figure24 where the dashed line is obtained by the hyperbolic decline model with qi=3.548412×104m3/d, Di=6.102143×10−4d−1, and n̄=0.071877. Figure 25 shows the estimated profile of Gj determined by Equation (35). Similarly, GHD of 2.650991436×108m3 coincides with Gsmb of 2.596831903×108m3.Figure 24
The realq profile and the predicted profile by the hyperbolic decline model for Case 2-2.Figure 25
The estimated profile ofGj for Case 2-2.


## 3.2.3. Case 2-3

Figure26 shows the real and predicted profiles of production rate for Case 2-3 with a BHP of 40 MPa in which the dashed line represents the matched hyperbolic decline model with qi=5.073846×104m3/d, Di=6.059150×10−4d−1, and n̄=0.137588. Then, the ultimate GHD of 2.667464699×108m3 determined by the decline parameter method, as shown in Figure 27, is comparable with Gsmb of 2.596426085×108m3 determined by the gpave vs. Gp curve.Figure 26
The realq profile and the predicted profile by the hyperbolic decline model for Case 2-3.Figure 27
The estimated profile ofGj for Case 2-3.


## 3.2.4. Case 2-4

As mentioned earlier, Equation (30) is used to match the rate profile of the first 6000 days for Case 2-4 with a constant BHP of 24 MPa, and three model parameters are determined as follows: qi=7.120844×104m3/d, Di=6.082153×10−4d−1, and n̄=0.258816. The production rate predicted by the hyperbolic decline model is in line with the real production profile as shown in Figure 28. The decline parameter method is applied, and Figure 29 displays the estimated reserve profile with the resulting GHD of 2.664141330×108m3 in line with Gsmb of 2.596254338×108m3 determined by the SMBE.Figure 28
The realq profile and the predicted profile by the hyperbolic decline model for Case 2-4.Figure 29
The estimated profile ofGj for Case 2-4.The computation performance of the decline parameter method for the synthetic gas reservoir with high pressure is summarized in Table6(a). The corresponding analysis results of the modified Stumpf method are displayed in Table 6(b) and Figures 35 and 36 (see Appendix). Two methods both give similar estimations of gas reserves, but the latter shows an increasing error trend with a decrease in pwf.Table 6(a)
The analysis results of the decline parameter method for numerical cases with high pressure
High-pressure casesCase 2-1Case 2-2Case 2-3Case 2-4pwf (MPa)64524024qi (104 m3/d)2.0272583.5484125.0738467.120844Di (10-4 d-1)6.1933296.1021436.0591506.082153n̄0.0300780.0718770.1375880.258816GHD (108 m3)2.6311259062.6509914362.6674646992.664141330Error1 (%)1.3782.1442.7782.650m (10-3 MPa/104 m3)2.0289782092.0283838452.0287008792.028835081Gsmb (108 m3)2.5960711942.5968319032.5964260852.596254338Error2 (%)0.0270.0570.0410.035(b)
The analysis results of modified Stumpf method for numerical cases with high pressure
High-pressure casesCase 2-1Case 2-2Case 2-3Case 2-4pwf (MPa)64524024qi (104 m3/d)2.0331033.5639585.0982277.158882Di (10-4 d-1)6.2808796.2297616.1949826.235047na0.0428620.0906990.1572940.281530Gm (108 m3)2.6101385522.6230317362.6394887272.673866274Error3 (%)0.5691.0661.7003.025It is found, from the previous synthetic cases with whether normal pressure or high pressure, that an equivalent decline exponentn̄ is capable of matching the production decline data of the gas well with a constant BHP though the instantaneous decline exponent n as a function of pave changes slowly with time. The numerical cases fully demonstrate that
(1)
The decline parameters of HDM are affected by the BHP of the gas well, and the initial production rate during BDFqi and the equivalent decline exponent n̄ are more sensitive to BHP values than the initial decline rate Di. The lower the bottomhole flowing pressure the gas well is produced at, the greater qi and n̄ are, while the change in Di is little obvious(2)
The flow integral method, an explicit approach for calculations of average reservoir pressure based on the flow integral equation, can estimate thepave profile effectively and conveniently because only one data point during BDF is necessary, and in reality, these estimations are not confined to the constant-BHP situations(3)
The decline parameter method, an explicit technique for calculations of gas in place (or gas reserves) based on the conventional decline parameters in HDM can reliably evaluate reserves of the production system with a constant BHP. Furthermore, the value of gas reserves calculated by the static material balance equation is available for the mutual verification with the estimate of the decline parameter method(4)
The production data matching in the decline parameter method could be superior to the analytic calculations ofqi, Di, and na in the modified Stumpf method based on the straight line of q1−na vs. Gp which is only suitable for early BDF data. The former extends the range of application of the latter. These approaches, however, can be verified against each other(5)
The decline parameter method is valid for both high-pressure and normal pressure gas reservoirs. The calculation errors for the high-pressure system, however, are slightly larger than those of the low-pressure system which may be attributed to the more decrease in formation pressure and a bit severer fluctuation inb during BDF for high-pressure gas reservoirs


## 3.3. Field Case: Gas Well A in West Virginia

Case 3, a field example provided by Fetkovich et al. [18], reflects the production scenario of Gas Well A in West Virginia. Fraim and Wattenbarger [49] also investigated this case and digitized the raw rate-time data as shown in Figure 30. The hyperbolic decline model is used to fit the rate profile, and the resulting three parameters are as follows: qi=6.509603×104m3/d, Di=2.534238×10−3d−1, and n̄=0.907180.Figure 30
The realq profile and the predicted profile by the hyperbolic decline model for Case 3.Since the specific gravity of the natural gas is 0.57, we obtain the pseudocritical temperatureTpc and pressure ppc are 191.22 K and 4.61 MPa, respectively, according to Sutton [44, 45]. Other property parameters are shown in Table 7.Table 7
Reservoir and gas properties for Case 3.
Property parametersProperty valuesProperty parametersProperty valuespi28.785612 MPah21.336 mTi344.26 Krw0.1079 mpwf3.447379 MPaρsc0.687754 kg/m3ϕi0.06Mg16.509480 g/molSwci0.35Zsc0.997913Tsc293.15 KZi0.952170psc0.101325 MPaμi2.209414019×10−2cpTpc191.218772 Kμsc1.100728475×10−2cpppc4.605256 MPaCgi2.589484787×10−2MPa−1Cϕ8.816×10−4MPa−1Cti1.785752111×10−2MPa−1Cw4.122×10−4MPa−1Bgi3.944216913×10−3m3/m3As revealed by the dotted line in Figure31, the relation of p/μ·Z to p for Case 3 can be represented by the following polynomial in the pressure range of 0.2 to 30:
(47)pμZ=−1.612592×10−6p6+1.292460×10−4p5−2.355782×10−3p4−3.709335×10−2p3−9.545551×10−2p2+8.033991×101p+9.824220×10−2.Figure 31
Relation curves ofp/μZ and μZ vs. pressure (p) for Case 3.According to Fetkovich et al., the average reservoir pressure at 200 days is 3268 psi (or 22.532067 MPa), and the corresponding production rate is estimated to be4.289933×104m3/d by the matched hyperbolic model. Then, the formation pressure profile is determined by the flow integral method. Figure 32 displays the estimated reserve profile by implementing the above steps of the decline parameter method.Figure 32
The estimated profile ofGj for Case 3.The average value ofGj for Case 3 leads to a GHD of 0.857228577×108m3; however, Gsmb is not estimated because of unavailable cumulative production. The negligible difference between GHD in the paper and the estimate for gas in place of 0.303451×107 Mscf (or 0.859277542×108m3) predicted by Fraim and Wattenbarger [49] using normalized time matching demonstrates the applicability and feasibility of the decline parameter method.


## 4. Conclusions

The flow integral equation during BDF is derived in the paper from the dynamic material balance equation of the gas reservoir; then, the flow integral method based on it is proposed for explicit determination of the average reservoir pressure profile. The decline parameter method, a novel approach for explicit estimation of gas reserves, is developed to evaluate the reserves of the gas well at a constant BHP on the basis of the conventional hyperbolic decline model and the predicted formation pressure data by the flow integral method; and the static material balance equation, on the other hand, can be used to verify the results. These methods are suitable for both abnormally pressured gas reservoirs and normal pressure systems because the compressibilities of rock and bound water are considered in the rationales. The following conclusions are drawn:(i)
The hyperbolic decline model with a constant decline exponentn̄ can accurately delineate the production rate data of the gas well at a constant BHP, although the actual value n varies slowly with time(ii)
The values of BHP influence the decline parameters (qi, Di, and n̄) in HDM, and the initial decline rate Di is relatively insensitive to BHPs compared to the other two parameters. The smaller the BHP is, the larger qi and n̄ are(iii)
The flow integral method, not limited to the constant-BHP conditions, is effective in estimating the profile of average reservoir pressure for various gas production systems during BDF in an explicit fashion(iv)
The decline parameter method, capturing the intrinsic correlation between production rate, decline rate, formation pressure, and gas reserves, is mainly applicable to the production well with a constant BHP and can also provide results of reference for those gas wells with slight fluctuations in BHP(v)
The production data matching in the decline parameter method could be superior to the analytic calculations ofqi, Di, and na in the modified Stumpf method based on the straight line of q1−na vs. Gp which is limited to relatively early BDF period. These approaches, nonetheless, can be verified against each other(vi)
The static material balance equation, combined with the estimated formation pressure data by the flow integral method, is also capable of determining the reserves accurately, which can also be employed to verify the decline parameter method



---

*Source: 5144920-2021-12-16.xml*