# An Explicit Method for Estimating the Average Reservoir Pressure and Gas in Place Based on the Gas Property Polynomial Function

Authors: Lixia Zhang; Hao Jiang; Yong Li; Xinmin Song; Yang Li
Journal: Geofluids
 (2023)

---

## Abstract

The fast and accurate determination of average reservoir pressure (ARP) and gas reserves is vital for the analysis and forecasting of gas well performance. An explicit and relatively rigorous method for estimating the ARP and gas in place, based on the dynamic material balance equation and the gas property polynomial function (GPPF), is presented to circumvent such drawbacks as high cost, empiricism, poor adaptability to variable production schedules, and the necessity for iteration processes, by which current approaches usually have been limited. The gas flow model is solved by introducing pseudofunctions and considering the compressibility effects of rock and irreducible water, and the pressure-rate correlations during boundary-dominated flow (BDF) which constitute the theoretical proofs of this method are derived from the superposition principle coupled with the constant rate solution. Production data under different production scenarios prove it effective. The error of estimation for formation pressure and gas reserves hardly ever goes beyond 4% given that the BDF condition is generally satisfied by the rate and bottom hole flowing pressure (BHP) data. We employ the GPPF to capture the nonlinear variations of gas viscosity andZ-factor, helpful to implement a quick conversion between pressure and pseudopressure and to solve the integral equation concerning pseudopressure and thus valid in overcoming the limitations of previous methods concerning pressure or pressure squared. The proposed methodology boasts its simpleness and practicability, which not only dispenses with iterations on pseudotime but also applies to various production systems of gas wells under constant BHP, constant rate, or variable BHP/variable rate conditions.


---

## Body

## 1. Introduction

The determination of ARP is an essential part of the analysis and prediction of well performance. Gas reservoir pressure data can be related to cumulative gas production through the static material balance equation (SMBE). Gas reserves can be estimated directly by that relationship if pressures are obtained by multiple measurements which are unfortunately time-consuming and have a high cost. The ARP profile in turn can also be determined by the SMBE when reserves are available.Generally, the ARP is extrapolated by the pressure transient analysis (PTA) [1–4], including the methods proposed by Horner [5], Miller et al. [6], Matthews et al. [7], Dietz [8], Ramey and Cobb [9], and Odeh and Al-Hussainy [10]. Formation pressure can also be estimated by empirical exponential [11] or binomial [12, 13] productivity equations using systematic well test (flow-after-flow test) data and can be measured directly after shutting in the well. These methods, however, are not suitable to be used frequently due to their high cost, which brings difficulties in continuous pressure estimation. Researchers, hence, began to explore cheaper methods for production data analysis (PDA).PDA involves the estimation of reservoir information through theoretical or empirical models by using daily measurement data such as rate, pressure, and cumulative production. The rate transient analysis (RTA), receiving wide attention due to its economical and practical characteristics, is an important type of PDA. The research on the RTA can be traced back to the production decline analysis by Arps [14] using statistics. The reported method is simple and applies to constant BHP production systems, but it was not proved theoretically at that time.Fetkovich [15, 16] first introduced the idea of PTA into RTA. He incorporated the constant BHP analytical solution of the unsteady flow and Arps’ curves into the log-log plots and performed the production decline analysis by a type curve fitting similar to PTA. Fetkovich proved theoretically that the relationship between the dimensionless rate qdD and the dimensionless time tdD for constant pressure systems follows the exponential decline law in late BDF, which, as a representative of RTA, overcomes the drawback (empiricism) of Arps’ curves. However, Fetkovich curves can only be applied to the constant pressure production systems of slightly compressible fluids for ignoring the changes in fluid properties.To consider the variations of gas properties during depletion, Carter [17, 18] introduced the variable, λ=μiCgi/μCg, to explain changes in viscosity and compressibility, and presented Carter curves to analyze the gas well data. The pressure drop parameter, λ, nevertheless, is an average from BHP pwf to the original reservoir pressure pi and does not change with the ARP. Therefore, it is only an approximate method for analyzing gas well data under the constant BHP condition.To tackle fluctuations in production schedules and changes in gas properties, pseudopressure and pseudotime functions were introduced by Blasingame and Lee [19], and he proposed, based on the pressure-rate relation in BDF and SMBE for constant volume gas reservoirs, a method for iterative calculation of reserves and pseudotime, where the ARP can be estimated by the SMBE.The material balance pseudotime functiontca can greatly eliminate the influence of production fluctuations on the analysis, but its determination requires iterations since the ARP is implicit. Blasingame et al. [20] proposed the concept of constant BHP equivalent time function. One can use the constant BHP solution to analyze data under variable BHP/variable rate conditions. Its innovation is the conversion from the constant rate analog time (tca) into the constant pressure analog time function; nonetheless, for gas wells, determining the parameters, mbdfandb, still requires iterations on tca.Palacio and Blasingame [21], based on the previous BDF pressure-rate correlation [19], proved that the relationship between qdDandtdD defined by the pseudopressure and tca is always consistent with the harmonic decline, thus developing a type curve analysis method that can explain variable rate/variable pressure data. As another representative of RTA, Blasingame’s type curves boast the ability to explain the changes in production schedules and fluid properties during BDF.In addition, other type curve analysis methods such as Agarwal-Gardner curves [22, 23] and normalized pressure integral (NPI) curves [24] actually change the drawing functions, while their analysis principles are basically the same.Aiming at constant rate gas wells in volumetric gas reservoirs, Mattar and McNeil [25, 26] proposed the flowing material balance (FMB) method for estimating reserves by using the ratio p/Z at the BHP to replace the value at the ARP. In fact, the slope of the straight line plotted by p/Z and the cumulative gas production when the reference pressure is the BHP is not equal to the result at ARP. Later, some researchers [27, 28] improved it, but their approaches still could not meet the requirements for analyzing variable rate/variable BHP data.To deal with fluctuations in flow rate, Mattar and Anderson [29, 30] improved the previous FMB, based on the “stabilized flow” solution and the SMBE, and developed the dynamic material balance (DMB) method (or variable rate FMB), which considers changes in property parameters and production scenarios. In terms of theory, it is more rigorous but requires iterations on the ARP and tca.Those methods for determining the ARP or reserves are usually implicit sincetca needs to be evaluated iteratively at the ARP when one analyzes gas production data using the classical PDA/RTA methods.Later, some researchers began to explore explicit methods. The research of Ye and Ayala [31] is quite representative. They ignored the viscosity term in the pseudopressure and introduced gas density (ρ) and pseudotime factor (β) to describe the gas flow. The dimensionless boundary, reD, is determined by the solution characterizing the liquid flow due to its small deviation away from the solution of gas flow in the early stage when gas property changes are small, while the values of λ, i.e., μiCgi/μCg, and β are estimated by the large deviation between the two solutions during BDF. Gas reserves are determined by the approximate relation between λ and ρ. This work assumes a constant straight-line slope in the μ·Cg vs. ρ logarithmic curve, and it is mainly targeted at the constant pressure production systems.Stumpf and Ayala [32] explained the physical connotations of the Arps’ hyperbolic-decline model for constant BHP gas wells and proposed a method for determining reserves using the rate vs. cumulative production data in early BDF (the hyperbolic window). It is necessary to identify the start/end moments of the hyperbolic window through curve fitting and determine the integral average of the decline exponent, which is cumbersome and only applicable to constant BHP conditions.Zhang et al. [33] presented an approach to determine reserves without pseudotime iterations by combining the material balance condition derived from the SMBE with the pseudopressure approximation condition derived from the DMB equation.Following Stumpf and Ayala’s research [32], Wang and Ayala [34] extended the hyperbolic decline model by postulating a constant ratio of the pseudo-BHP to the pseudo-ARP and determined reserves by using curve fitting and linear analysis of the rate vs. cumulative production curve. Both Stumpf and Wang used the integral mean from the pseudo-BHP to the pseudo-ARP when calculating the decline exponent, but it is actually a function of ARP. It seems more reasonable that the integration range is evaluated by the variation range of ARP during the target period.Jongkittinarukorn et al. [35] also discussed the changes in decline rate and decline exponent for constant BHP gas wells in volumetric reservoirs. They explained the relationship between production rate, decline rate, decline exponent, and λ vs. p/Z/pi/Zi based on the stabilized flow equation developed by Ansah et al. [36, 37] and estimated gas reserves by that explanation. The workflow, which infers the decline rate profile by the rate-time profile and then determines the decline exponent profile, however, may cause undesirable results.In conclusion, currently available PDA methods for determining gas reservoir pressure and reserves mainly include type curve analysis [38, 39], FMB [25–28], pseudotime iteration methods (such as variable rate reservoir limits testing developed by Blasingame and Lee [19] and DMB developed by Mattar and Anderson [29, 30]), and various explicit methods. Among them, Carter’s type curves are limited to constant pressure conditions. Furthermore, variable rate reservoir limits testing, DMB, and Blasingame’s type curves can dispose fluctuations in production conditions, but the pseudotime needs to be calculated iteratively. The modified FMB approach is only suitable for constant rate gas wells, while explicit methods that rely on traditional decline parameters can merely handle constant BHP conditions or specific pressure conditions [32, 34, 35]. In addition, most of these methods do not consider the compressibility of pores and irreducible water during pressure depletion.Therefore, based on the DMB equation and the gas property polynomial function (GPPF), this research is aimed at developing an explicit method for estimating the ARP and gas in place under variable BHP/variable rate conditions, which also captures the similar superiority of implicit techniques and considers the compressibility effects.


## 2. Model Development

The developed methodology here is based on the pressure-rate relationships during BDF (especially the dynamic material balance equation, DMBE) for a gas well in a bounded reservoir and the gas property polynomial function for pseudopressure approximation. This paper takes the radial reservoir system as an example to derive the DMBE also applicable to noncircular boundary reservoirs. The theoretical bases of the proposed method for the determination of ARP and gas reserves will be demonstrated below.


### 2.1. Pressure-Rate Correlations in BDF

The gas flow through porous media conforms to Darcy’s law, namely,(1)v⟶g=−Kμ∇p.The state equations of rock [40, 41] are represented by
(2)ϕ=ϕi⋅eCϕp−pi,(3)Vp=Vpi⋅eCϕp−pi.The state equations of irreducible water can be written as(4)Vw=Vwi⋅e−Cwp−pi,(5)Swc=VwVp=Swci⋅e−Cw+Cϕp−pi.The state of gas is mainly characterized by its density function and gas compressibility:(6)ρg=pMZRT,(7)Cg=ZTp⋅ddppZT≈Zp⋅ddppZ.The mass conservation equation of gas can be obtained by the combination of the motion equation and the state equations, i.e.,(8)−∇⋅ρgv⟶g+q̇g=1−Swc∂ρgϕ∂t+ρgϕ∂1−Swc∂t.Substituting Equations (1)–(7) into Equation (8), the governing equation of gas flow in the porous media can be expressed as
(9)∇⋅pμZ∇p=ϕiμiCtiK⋅μCtμiCti⋅pμZ∂p∂t,(10)Ct=eCϕp−pi1−SwcCg+Cϕ+SwcCw.It is difficult to directly solve the nonlinear differential equation, Equation (9), because the gas viscosity, the deviation factor, and the gas compressibility are all functions of pressure and temperature. Even if under the isothermal condition, they all change with pressure. Hence, the pseudofunctions of real time and real pressure are introduced:
(11)ta=μiCti∫0t1μCtdt,(12)pp=pi+μiZipi∫pipξμξZξdξ.By substituting Equations (11) and (12) into Equation (9), one obtains
(13)∇2pp=ϕiμiCtiK⋅∂pp∂ta.To simplify the calculation of the integral term in the pseudopressure in the field practice, the integrandfp in the two cases of high pressure and low pressure is often approximated as follows, respectively:
(14)fp=pμpZp≈pC0p<13.8MPa,1C1p>20.7MPa,Fp=∫pμpZpdp=12C0p2+C0′p<13.8MPa,1C1p+C1′p>20.7MPa.It should be noted that this approach of using pressure or pressure squared to approximate pseudopressure is difficult to capture the pressure drop process that spans two pressure ranges and may cause large errors. Therefore, the mathematical model of gas flow in a circular bounded reservoir under the constant production rate condition is given in the rigorous form of pseudopressure (in the international system of units,αt=1,αp=2π):
(15)∂2pp∂r2+1r⋅∂pp∂r=ϕiμiCtiαtK⋅∂pp∂ta,r∂pp∂rr=rw=qμiBgiαpKh,∂pp∂rr=re=0,ppta=0=pppi.By solving nondimensionalized Equation (15), the dimensionless pressure-time relationship is given by
(16)pD=2tDreD2−1−reD2reD2−1lnrD+4reD4lnreD−3reD4+2rD2reD2−1+2reD2+14reD2−12−π∑n=1∞e−λn2⋅tD⋅J12reDλnY1λnJ0rDλn−J1λnY0rDλnλn⋅J12reDλn−J12λn,(17)rD=rrw,(18)tD=αtKϕiμiCti⋅rw2ta,(19)pD=αpKhppi−ppqμiBgi.According to the superposition principle, pseudopressure drop caused by production fluctuations can be written as(20)ppi−ppr,ta=∑i=1mμiBgiαpKh⋅qi−qi−1⋅pDta−ta,i−1.Substituting Equation (16) into Equation (20), we obtain
(21)ppi−ppwfq=mbdf⋅tca+b,(22)mbdf=2παtαp⋅BgiAhϕiCti,(23)tca=μiCtiqt∫0tqtμpaveCtpavedt,(24)b=μiBgiαpKh4reD4lnreD−3reD4+4reD2−14reD2−12+∑i=1mqi−qi−1qm∑n=1∞2λn2⋅lta,λn,(25)lta,λn=e−λn2⋅αtK/ϕiμiCtirw2ta−ta,i−1J12reDλnJ12reDλn−J12λn.Equation (21) reveals the relationship between BHP, production rate, and material-balance pseudotime where tca is the reference value under the ARP. The average pseudopressure of the gas reservoir is defined as
(26)ppave=1πre2−rw2∫rwrepp⋅2πr⋅dr.Differentiating Equation (20) with respect to pseudotime ta and combining it with Equation (15), we can get
(27)pp−ppwf=qmμiBgiαpKhre2−rw2re2lnrrw−r2−rw22−μiBgiαpKh×∑i=1mqi−qi−1∑n=1∞lta,λn−πλnY1λnJ0λnrwr+πλnJ1λnY0λnrwr−2λn2.Substituting Equation (27) into Equation (26) gives
(28)ppave=ppwf+qmμiBgiαpKh⋅4reD2+4reD4lnreD−3reD4−14reD2−12+μiBgiαpKh∑i=1mqi−qi−1∑n=1∞lta,λn⋅2λn2.By applyingpppave=ppave≈ppave, Equation (28) can be rewritten as
(29)ppave=ppwf+q⋅b.Equation (29), termed “dynamic material balance equation” (DMBE), reveals the internal relationship between the ARP, BHP, and production rate. The value of b in Equations (21) and (29) changes very slowly with time in the BDF stage and can be approximated to a certain constant during the evaluation period; that is to say, the “boundary dominated flow condition” holds.


### 2.2. Determination of ARP

According to the definition of pseudopressure Equation (12), one can rewrite Equation (29) as
(30)q=ppave−ppwfb=1b⋅μiZipi∫pipavefξdξ−∫pipwffξdξ.Knowing the average formation pressure (pk) at some time (tk), we have
(31)qk=1b⋅μiZipi∫pwf,kpkfξdξ.From Equations (30) and (31), it follows
(32)∫pwfpavefξdξ=qqk∫pwf,kpkfξdξ.Equation (32) is an integral equation derived from the DMBE. According to the reservoir pressure at a certain moment in the BDF stage and the recorded production data of the gas well (such as production rate and BHP), the reservoir pressure at each moment can be inferred by solving Equation (32). Directly using the pressure or pressure square method described above to solve it may cause large errors, while using numerical integration (such as the Romberg algorithm) may be time-consuming. Hence, this paper uses a polynomial to approximate the integrand f, namely,
(33)fp=pμpZp≈fpnp=c6p6+c5p5+c4p4+c3p3+c2p2+c1p+c0.So, the original function offp can be expressed as
(34)Fp=∫pμpZpdp≈c67p7+c56p6+c45p5+c34p4+c23p3+c12p2+c0p+C.For convenience, the arbitrary constantC is taken to be zero in the paper. Equation (32) can be expressed analytically as
(35)Fpave=qqkFpk−Fpwf,k+Fpwf.The left term of Equation (35), Fpave, is the polynomial function dependent on gas properties. Knowledge of the right term of Equation (35) can solve the values of pave in the interval of pwf,pi by the polynomial rooting, bisection, or Newton iteration method.


### 2.3. Estimation of Gas Reserves

The material balance for a bounded gas reservoir involves the elastic expansion of natural gas, the shrinkage of pores, and the expansion of irreducible water. The three elastic energy supplies the output of natural gas, namely,(36)GBg−GBgi+ΔVp+ΔVw=GpBg.The gas volume factor,BgiandBgpave, can be expressed as
(37)Bgi=ZiTipi⋅pscZscTsc,Bgpave=ZTpave⋅pscZscTsc.The volume of irreducible water and rock pores under original formation conditions can be written as(38)Vwi=GBgi1−SwciSwci,Vpi=GBgi1−Swci.It is reasonably assumed that the gas reservoir temperature remains unchanged during gas production. Substituting Equations (3), (4), (37), and (38) into Equation (36) yields
(39)paveZpave⋅eCϕpave−pi−Swcie−Cwpave−pi1−Swci=piZi1−GpG.For expression convenience, an elastic effect functiong related to pressure and compressibilities (of the rock and irreducible water) is introduced as follows:
(40)gp=pZp⋅eCϕp−pi−Swcie−Cwp−pi1−Swci.So, Equation (39) can be reduced to
(41)gpave=−piZi⋅1G⋅Gp+piZi,(42)msmb=piZi⋅1G.Equation (41), termed “static material balance equation” (SMBE), reflects the relationship between cumulative gas production and average reservoir pressure. The ARP profile can be inferred from the GPPF fpnp and production data, and the gas reserves can be estimated by the negative slope (msmb) of the gpavevs.Gp curve.


### 2.4. Main Steps for ARP Estimation

The pressure-rate correlations during BDF (including the dynamic material balance equation) are derived rigorously from the gas flow theory for the closed gas reservoirs, and the theoretical bases for the developed methodology for the determination of ARP and gas reserves are also explained in the earlier article. Note that the shape factor [19, 33] can be introduced into Equations (21) and (29) to delineate the gas well performance in noncircular reservoirs, and it will not be discussed in detail here.The main steps of the explicit approach presented in the paper are described below. Firstly, the property parameters of the reservoir and natural gas are needed such as reservoir temperature, initial pressure, compressibility (Cw, Cϕ), and pseudocritical temperature and pressure of gas. So, the gas properties (such as viscosity, Z-factor, and compressibility Cg) can be determined by corresponding correlations, and the behavior of the elastic effect function gp with pressure can also be described by Equation (40). Secondly, it is easy to establish the relation of fp or p/μ·Z to p, and one can use the polynomial function fpnp to replace fp; then, Fp, as the primitive function of fp, can be determined. Thirdly, gas production data are necessary such as rate, BHP, and cumulative gas production. Then, we need to choose a reference point in BDF to estimate the profile of ARP by Equation (35). Lastly, gas reserves can be inferred from the gpave vs. Gp straight line according to Equation (41) by which in turn the obtained reserves are available to determine the formation pressure profile. The key procedures are summarized in Figure 1.Figure 1
The workflow for estimation of ARP and reserves.


## 2.1. Pressure-Rate Correlations in BDF

The gas flow through porous media conforms to Darcy’s law, namely,(1)v⟶g=−Kμ∇p.The state equations of rock [40, 41] are represented by
(2)ϕ=ϕi⋅eCϕp−pi,(3)Vp=Vpi⋅eCϕp−pi.The state equations of irreducible water can be written as(4)Vw=Vwi⋅e−Cwp−pi,(5)Swc=VwVp=Swci⋅e−Cw+Cϕp−pi.The state of gas is mainly characterized by its density function and gas compressibility:(6)ρg=pMZRT,(7)Cg=ZTp⋅ddppZT≈Zp⋅ddppZ.The mass conservation equation of gas can be obtained by the combination of the motion equation and the state equations, i.e.,(8)−∇⋅ρgv⟶g+q̇g=1−Swc∂ρgϕ∂t+ρgϕ∂1−Swc∂t.Substituting Equations (1)–(7) into Equation (8), the governing equation of gas flow in the porous media can be expressed as
(9)∇⋅pμZ∇p=ϕiμiCtiK⋅μCtμiCti⋅pμZ∂p∂t,(10)Ct=eCϕp−pi1−SwcCg+Cϕ+SwcCw.It is difficult to directly solve the nonlinear differential equation, Equation (9), because the gas viscosity, the deviation factor, and the gas compressibility are all functions of pressure and temperature. Even if under the isothermal condition, they all change with pressure. Hence, the pseudofunctions of real time and real pressure are introduced:
(11)ta=μiCti∫0t1μCtdt,(12)pp=pi+μiZipi∫pipξμξZξdξ.By substituting Equations (11) and (12) into Equation (9), one obtains
(13)∇2pp=ϕiμiCtiK⋅∂pp∂ta.To simplify the calculation of the integral term in the pseudopressure in the field practice, the integrandfp in the two cases of high pressure and low pressure is often approximated as follows, respectively:
(14)fp=pμpZp≈pC0p<13.8MPa,1C1p>20.7MPa,Fp=∫pμpZpdp=12C0p2+C0′p<13.8MPa,1C1p+C1′p>20.7MPa.It should be noted that this approach of using pressure or pressure squared to approximate pseudopressure is difficult to capture the pressure drop process that spans two pressure ranges and may cause large errors. Therefore, the mathematical model of gas flow in a circular bounded reservoir under the constant production rate condition is given in the rigorous form of pseudopressure (in the international system of units,αt=1,αp=2π):
(15)∂2pp∂r2+1r⋅∂pp∂r=ϕiμiCtiαtK⋅∂pp∂ta,r∂pp∂rr=rw=qμiBgiαpKh,∂pp∂rr=re=0,ppta=0=pppi.By solving nondimensionalized Equation (15), the dimensionless pressure-time relationship is given by
(16)pD=2tDreD2−1−reD2reD2−1lnrD+4reD4lnreD−3reD4+2rD2reD2−1+2reD2+14reD2−12−π∑n=1∞e−λn2⋅tD⋅J12reDλnY1λnJ0rDλn−J1λnY0rDλnλn⋅J12reDλn−J12λn,(17)rD=rrw,(18)tD=αtKϕiμiCti⋅rw2ta,(19)pD=αpKhppi−ppqμiBgi.According to the superposition principle, pseudopressure drop caused by production fluctuations can be written as(20)ppi−ppr,ta=∑i=1mμiBgiαpKh⋅qi−qi−1⋅pDta−ta,i−1.Substituting Equation (16) into Equation (20), we obtain
(21)ppi−ppwfq=mbdf⋅tca+b,(22)mbdf=2παtαp⋅BgiAhϕiCti,(23)tca=μiCtiqt∫0tqtμpaveCtpavedt,(24)b=μiBgiαpKh4reD4lnreD−3reD4+4reD2−14reD2−12+∑i=1mqi−qi−1qm∑n=1∞2λn2⋅lta,λn,(25)lta,λn=e−λn2⋅αtK/ϕiμiCtirw2ta−ta,i−1J12reDλnJ12reDλn−J12λn.Equation (21) reveals the relationship between BHP, production rate, and material-balance pseudotime where tca is the reference value under the ARP. The average pseudopressure of the gas reservoir is defined as
(26)ppave=1πre2−rw2∫rwrepp⋅2πr⋅dr.Differentiating Equation (20) with respect to pseudotime ta and combining it with Equation (15), we can get
(27)pp−ppwf=qmμiBgiαpKhre2−rw2re2lnrrw−r2−rw22−μiBgiαpKh×∑i=1mqi−qi−1∑n=1∞lta,λn−πλnY1λnJ0λnrwr+πλnJ1λnY0λnrwr−2λn2.Substituting Equation (27) into Equation (26) gives
(28)ppave=ppwf+qmμiBgiαpKh⋅4reD2+4reD4lnreD−3reD4−14reD2−12+μiBgiαpKh∑i=1mqi−qi−1∑n=1∞lta,λn⋅2λn2.By applyingpppave=ppave≈ppave, Equation (28) can be rewritten as
(29)ppave=ppwf+q⋅b.Equation (29), termed “dynamic material balance equation” (DMBE), reveals the internal relationship between the ARP, BHP, and production rate. The value of b in Equations (21) and (29) changes very slowly with time in the BDF stage and can be approximated to a certain constant during the evaluation period; that is to say, the “boundary dominated flow condition” holds.


## 2.2. Determination of ARP

According to the definition of pseudopressure Equation (12), one can rewrite Equation (29) as
(30)q=ppave−ppwfb=1b⋅μiZipi∫pipavefξdξ−∫pipwffξdξ.Knowing the average formation pressure (pk) at some time (tk), we have
(31)qk=1b⋅μiZipi∫pwf,kpkfξdξ.From Equations (30) and (31), it follows
(32)∫pwfpavefξdξ=qqk∫pwf,kpkfξdξ.Equation (32) is an integral equation derived from the DMBE. According to the reservoir pressure at a certain moment in the BDF stage and the recorded production data of the gas well (such as production rate and BHP), the reservoir pressure at each moment can be inferred by solving Equation (32). Directly using the pressure or pressure square method described above to solve it may cause large errors, while using numerical integration (such as the Romberg algorithm) may be time-consuming. Hence, this paper uses a polynomial to approximate the integrand f, namely,
(33)fp=pμpZp≈fpnp=c6p6+c5p5+c4p4+c3p3+c2p2+c1p+c0.So, the original function offp can be expressed as
(34)Fp=∫pμpZpdp≈c67p7+c56p6+c45p5+c34p4+c23p3+c12p2+c0p+C.For convenience, the arbitrary constantC is taken to be zero in the paper. Equation (32) can be expressed analytically as
(35)Fpave=qqkFpk−Fpwf,k+Fpwf.The left term of Equation (35), Fpave, is the polynomial function dependent on gas properties. Knowledge of the right term of Equation (35) can solve the values of pave in the interval of pwf,pi by the polynomial rooting, bisection, or Newton iteration method.


## 2.3. Estimation of Gas Reserves

The material balance for a bounded gas reservoir involves the elastic expansion of natural gas, the shrinkage of pores, and the expansion of irreducible water. The three elastic energy supplies the output of natural gas, namely,(36)GBg−GBgi+ΔVp+ΔVw=GpBg.The gas volume factor,BgiandBgpave, can be expressed as
(37)Bgi=ZiTipi⋅pscZscTsc,Bgpave=ZTpave⋅pscZscTsc.The volume of irreducible water and rock pores under original formation conditions can be written as(38)Vwi=GBgi1−SwciSwci,Vpi=GBgi1−Swci.It is reasonably assumed that the gas reservoir temperature remains unchanged during gas production. Substituting Equations (3), (4), (37), and (38) into Equation (36) yields
(39)paveZpave⋅eCϕpave−pi−Swcie−Cwpave−pi1−Swci=piZi1−GpG.For expression convenience, an elastic effect functiong related to pressure and compressibilities (of the rock and irreducible water) is introduced as follows:
(40)gp=pZp⋅eCϕp−pi−Swcie−Cwp−pi1−Swci.So, Equation (39) can be reduced to
(41)gpave=−piZi⋅1G⋅Gp+piZi,(42)msmb=piZi⋅1G.Equation (41), termed “static material balance equation” (SMBE), reflects the relationship between cumulative gas production and average reservoir pressure. The ARP profile can be inferred from the GPPF fpnp and production data, and the gas reserves can be estimated by the negative slope (msmb) of the gpavevs.Gp curve.


## 2.4. Main Steps for ARP Estimation

The pressure-rate correlations during BDF (including the dynamic material balance equation) are derived rigorously from the gas flow theory for the closed gas reservoirs, and the theoretical bases for the developed methodology for the determination of ARP and gas reserves are also explained in the earlier article. Note that the shape factor [19, 33] can be introduced into Equations (21) and (29) to delineate the gas well performance in noncircular reservoirs, and it will not be discussed in detail here.The main steps of the explicit approach presented in the paper are described below. Firstly, the property parameters of the reservoir and natural gas are needed such as reservoir temperature, initial pressure, compressibility (Cw, Cϕ), and pseudocritical temperature and pressure of gas. So, the gas properties (such as viscosity, Z-factor, and compressibility Cg) can be determined by corresponding correlations, and the behavior of the elastic effect function gp with pressure can also be described by Equation (40). Secondly, it is easy to establish the relation of fp or p/μ·Z to p, and one can use the polynomial function fpnp to replace fp; then, Fp, as the primitive function of fp, can be determined. Thirdly, gas production data are necessary such as rate, BHP, and cumulative gas production. Then, we need to choose a reference point in BDF to estimate the profile of ARP by Equation (35). Lastly, gas reserves can be inferred from the gpave vs. Gp straight line according to Equation (41) by which in turn the obtained reserves are available to determine the formation pressure profile. The key procedures are summarized in Figure 1.Figure 1
The workflow for estimation of ARP and reserves.


## 3. Model Validation

The reservoir simulator is used to generate production data including production rate, BHP, and cumulative gas production when production constraints of gas well are constant BHP, constant rate, and variable BHP/variable rate, respectively. The radial grid number is300×120×1. The initial formation pressure pi is 60 MPa, the original formation temperature Ti is 100°C, and the initial porosity ϕi is 0.19. The radial permeability Kr and Kθ are both 8 mD, while the vertical permeability Kz equals 0.8 mD. The initial irreducible water saturation Swci is 0.21, and the gas reservoir top depth is set to 2996 m with a reservoir thickness of 8 m. The equilibrium initialization indicates that Vpi=967412m3, G=271300553m3, and pi/Zi=45.908 MPa. The compressibilities of rock and irreducible water are set by empirical formulas [42, 43], and other gas reservoir parameters are shown in Table 1.Table 1
Property parameters for numerical simulations.
Property parametersProperty valuesϕi0.190Kr8 mDKθ8 mDKz0.800 mDSwci0.210h8 mdr1.500 mdθ3°dz8 mpi60 MPaTi373.150 KTsc293.150 Kpsc1.013×105 PaTpc205.760 Kppc4.590 MPare450 mrw0.100 mρsc0.778 kg/m3M18.660 g/molZsc0.997Zi1.307Cϕ5.335×10−4 MPa-1Cw4.002×10−4 MPa-1μw0.289 cpμi3.434×10−2 cpμsc1.064×10−2 cpCgi7.230×10−3 MPa-1Cti6.329×10−3 MPa-1Bgi2.817×10−3 m3/m3G271 300 553 m3The viscosity correlations proposed by Londono et al. [44, 45] are used to obtain the viscosity data of natural gas, and the Hall-Yarborough method [46] is used to calculate the deviation factor Z and the gas compressibility Cg. The Z-factor and gas compressibility curves generated by the Hall-Yarborough correlations are displayed in Figures 2 and 3, respectively, where Cpr in Figure 3 represents the pseudoreduced compressibility [47–50] given by
(43)Cg=Cprppc−1.Figure 2
(a) Compressibility factors (Z-factors) of natural gases obtained by the Hall-Yarborough correlations at pseudoreduced pressures of 0.2 to 15. (b) Compressibility factors (Z-factors) of natural gases obtained by the Hall-Yarborough correlations at pseudoreduced pressures of 15 to 30.
(a)(b)Figure 3
Pseudoreduced compressibility (Cpr) of natural gases obtained by the Hall-Yarborough correlations at pseudoreduced pressures of 0.2 to 30.The curves of gas formation volume factorBg vs. pressure p and viscosity μ vs. p are shown in Figure 4. Figure 5 reveals the changes in p/μZ and μZ with pressure p which also indicates that using pressure squared approximation may cause large errors, and μZ/p is only approximately a constant when the pressure is high.Figure 4
Gas formation volume factor (Bg) and gas viscosity (μ) vs. pressure (p) for case studies.Figure 5
p/μZ and μZ vs. pressure (p) for case studies.Figure6 displays the integrand in pseudopressure obtained by Equation (45) and its error range, where ε is defined as
(44)ε=fpn−ff×100%,(45)fpnp≈2.863218×10−7p6−6.009013×10−5p5+4.746580×10−3p4−1.595755×10−1p3+8.873368×10−1p2+7.437974×101p+2.473877.Figure 6
The relation curves ofp/μZ and ε vs. p calculated by the polynomial function fpn for case studies.It can be seen from Figure6 that the polynomial function, fpn, can accurately represent the values of fp or p/μZ in the range of 0.7~60 MPa, and the error is no more than 1%. Figure 7 reveals the relation of the pseudopressure ppp to pressure p and that of the elastic effect function gp to pressure p. All of them have the same dimension as the pressure.Figure 7
The relation curves of pseudopressure (pp) and gp vs. pressure (p) for case studies.


### 3.1. Constant BHP Production System

The gas well’s BHP is set to 28.884 MPa to simulate the 2000-day production process of the gas reservoir, in which the number of time steps adopted is 564. The production rate profile is shown in Figure8.Figure 8
Production history of the gas well producing at a constant BHP.Taking the production time of the 1000th day as the reference pointtk, we can get pk=47.477 MPa and qk=2.385×104 m3/d. Based on the rate and flowing pressure data, the ARP profile can be inferred by Equation (35), as shown in Figure 9.Figure 9
The real and estimated profiles of average reservoir pressure for the constant BHP gas well.After the gas reservoir pressures are inferred, the correlation between the elastic effect functiongpave and the cumulative gas production Gp can be determined by the gp vs. p curve. As shown in Figure 10, the SMBE shows msmb of 1.692×10−3 MPa/(104 m3), and the gas reserves are determined to be 2.713×108 m3, only with an error of -0.003%. If another known measured pressure in the BDF stage is selected as the reference point, the result will change slightly but will not cause a large error. It is recommended to use midsection data in the evaluation period as reference points.Figure 10
The relation ofgpave to Gp for the gas well producing at a constant BHP.


### 3.2. Constant Rate Production System

The production constraint of the gas well is set at a constant production rate, i.e.,q=2.500×104 m3/d, with other conditions unchanged. The simulated BHP profile is shown in Figure 11.Figure 11
Production history of the gas well producing at a constant production rate.Taking the production time of the 1000th day as the reference pointtk, then pk is 49.570 MPa with pwf,k of 30.053 MPa. Based on the production data, the ARP profile estimated by Equation (35) is illustrated in Figure 12.Figure 12
The real and estimated profiles of average reservoir pressure for the gas well producing at a constant production rate.Figure12 shows that the difference between the estimated pave and the actual ARP in the later simulation stage gradually increases due to the slow variations of b in the dynamic material balance equation. The change becomes relatively obvious with a relatively long production time. The static material balance correlation obtained from the estimated ARP is exhibited in Figure 13, indicating msmb=1.662×10−3MPa/104m3 and that subsequently the gas reserves are 2.762×108 m3 with an error of merely 1.814%. The results demonstrate that the proposed method for estimation of average reservoir pressure and gas in place is also applicable to production systems at a constant rate.Figure 13
The relation ofgpave to Gp for the gas well producing at a constant production rate.


### 3.3. Variable BHP and Variable Rate Production Systems

We keep the gas reservoir parameters and natural gas properties unchanged and periodically change the BHP and production rate of the gas well to simulate the complex variable BHP/variable rate production schedules. The production history of the gas well is shown in Figure14.Figure 14
Production history of the gas well under variable BHP and variable rate conditions.Similarly, taking the production time of the 1000th day as the reference point, we obtainpk=46.490MPa,qk=2.977×104m3/d,andpwf,k=23.290MPa. Then, the estimated profile of ARP and the static material balance curve are displayed in Figures 15 and 16, respectively.Figure 15
The real and estimated profiles of average reservoir pressure for the gas well under variable BHP and variable production rate conditions.Figure 16
The relation ofgpave to Gp for the gas well under variable BHP and variable production rate conditions.The slope of thegpave vs. Gp curve, msmb, is 1.663×10−3 MPa/(104 m3); hence, the gas in place is estimated to be 2.760×108 m3 with an error of 1.744%. It can be observed from Figures 15 and 16 that the value gap of b between the initial stage and the reference point is slightly larger than that between the reference point and medium late stage due to fluctuations of gas well production, changes in BHP, and long production period. Selecting the midsection data in the evaluation period as reference points may alleviate the fluctuations of b.The above three numerical examples have effectively proved the validity of the explicit method for estimating the average reservoir pressure and gas reserves based on the GPPF in the paper. Only one measured formation pressure is needed to give an understanding of the pressure drop history and gas in place based on the proposed method, which is widely applicable to various production systems as long as the BDF condition is satisfied. It is worth mentioning that the estimation error of the ARP and reserves may increase slightly if the measured pressure data are available merely at the beginning or end of production or if the formation pressure at the reference point is inferred by the initial pressure, production rate, and BHP data.


### 3.4. Field Case

The field case (well H-58) [51] is derived from Ibrahim et al. [52] who employed the superposition normalized pseudotime correlation to smooth field data and more accurately determine gas reserves. This gas well is from a tight gas reservoir with high pressure. The gas in place calculated from their normalized pseudotime plotting function gives 9.8 Bcf. Table 2 shows the reservoir properties and fluid parameters. The formation temperature Ti is 93.333°C, the original formation pressure pi is equal to 96.527 MPa, and pi/Zi is 55.240 MPa.Table 2
Reservoir and gas properties for the field case.
Property parametersProperty valuesϕi1.354×10−1Swci2.618×10−1h114.148 mpi96.527 MPaTi366.483 KTsc293.150 Kpsc1.013×105 PaM19.268 g/molTpc200 Kppc4.600 MPaρsc0.803 kg/m3Zsc0.998Zi1.747Cϕ8.702×10−4 MPa-1Cw5.802×10−4 MPa-1μi4.752×10−2 cpμsc1.055×10−2 cpCgi3.478×10−3 MPa-1Cti3.590×10−3 MPa-1Bgi2.299×10−3 m3/m3Figure17 illustrates the variations of fp and μZ with pressure p which indicates a possible large error caused by the pressure squared approximation and that μZ/p is only approximately constant at relatively high pressure. Figure 18 shows the integrand in pseudopressure calculated by the GPPF (fpn) and its error where the expression of fpn fitted by the least square method is given by
(46)fpnp≈2.364940×10−8p6−6.219404×10−6p5+4.933695×10−4p4+1.569904×10−3p3−2.069461p2+9.564624×101p−2.970716×101.Figure 17
The relation curves ofp/μZ and μZ vs. pressure (p) for the field case.Figure 18
The relation curves ofp/μZ and ε vs. p calculated by the polynomial function fpn for the field case.It can be seen from Figure18 that the polynomial function, fpn, can accurately represent the solution of fp=p/μZ within the pressure range of 2.3~96.5 MPa, and errors of most of the results do not exceed 1%. Figure 19 shows the changes in ppp and gp with p. The production rate and bottom hole flowing pressure data of well H-58 are exhibited in Figure 20.Figure 19
The relation curves of pseudopressure (pp) and gp vs. pressure (p) for the field case.Figure 20
Production history of the gas well H-58.Taking the production time of the 250th day as the reference point, it is inferred thatpk=31.089MPa,qk=22.724×104m3/d,andpwf,k=10.217 MPa. The estimated formation pressure profile is displayed in Figure 21. It can be observed that the initial production is not stable, and the estimated formation pressure does not drop, which may be caused by the mismatches of b between the initial points and the reference point. Part of the estimated pressure data during that period should be removed when one determines the gas reserves. Figure 22 shows the static material balance relationship of the gas well H-58 where the slope of the gpave vs. Gp relationship curve msmb equals 2.017×10−3 MPa/(104 m3), and subsequently, the gas reserves are estimated to be 2.739×108 m3 which is close to the estimate of Ibrahim et al. [52] (9.8Bcf=2.775×108 m3).Figure 21
The estimated profile of average reservoir pressure for the gas well H-58.Figure 22
The relation ofgpave to Gp for the gas well H-58.


## 3.1. Constant BHP Production System

The gas well’s BHP is set to 28.884 MPa to simulate the 2000-day production process of the gas reservoir, in which the number of time steps adopted is 564. The production rate profile is shown in Figure8.Figure 8
Production history of the gas well producing at a constant BHP.Taking the production time of the 1000th day as the reference pointtk, we can get pk=47.477 MPa and qk=2.385×104 m3/d. Based on the rate and flowing pressure data, the ARP profile can be inferred by Equation (35), as shown in Figure 9.Figure 9
The real and estimated profiles of average reservoir pressure for the constant BHP gas well.After the gas reservoir pressures are inferred, the correlation between the elastic effect functiongpave and the cumulative gas production Gp can be determined by the gp vs. p curve. As shown in Figure 10, the SMBE shows msmb of 1.692×10−3 MPa/(104 m3), and the gas reserves are determined to be 2.713×108 m3, only with an error of -0.003%. If another known measured pressure in the BDF stage is selected as the reference point, the result will change slightly but will not cause a large error. It is recommended to use midsection data in the evaluation period as reference points.Figure 10
The relation ofgpave to Gp for the gas well producing at a constant BHP.


## 3.2. Constant Rate Production System

The production constraint of the gas well is set at a constant production rate, i.e.,q=2.500×104 m3/d, with other conditions unchanged. The simulated BHP profile is shown in Figure 11.Figure 11
Production history of the gas well producing at a constant production rate.Taking the production time of the 1000th day as the reference pointtk, then pk is 49.570 MPa with pwf,k of 30.053 MPa. Based on the production data, the ARP profile estimated by Equation (35) is illustrated in Figure 12.Figure 12
The real and estimated profiles of average reservoir pressure for the gas well producing at a constant production rate.Figure12 shows that the difference between the estimated pave and the actual ARP in the later simulation stage gradually increases due to the slow variations of b in the dynamic material balance equation. The change becomes relatively obvious with a relatively long production time. The static material balance correlation obtained from the estimated ARP is exhibited in Figure 13, indicating msmb=1.662×10−3MPa/104m3 and that subsequently the gas reserves are 2.762×108 m3 with an error of merely 1.814%. The results demonstrate that the proposed method for estimation of average reservoir pressure and gas in place is also applicable to production systems at a constant rate.Figure 13
The relation ofgpave to Gp for the gas well producing at a constant production rate.


## 3.3. Variable BHP and Variable Rate Production Systems

We keep the gas reservoir parameters and natural gas properties unchanged and periodically change the BHP and production rate of the gas well to simulate the complex variable BHP/variable rate production schedules. The production history of the gas well is shown in Figure14.Figure 14
Production history of the gas well under variable BHP and variable rate conditions.Similarly, taking the production time of the 1000th day as the reference point, we obtainpk=46.490MPa,qk=2.977×104m3/d,andpwf,k=23.290MPa. Then, the estimated profile of ARP and the static material balance curve are displayed in Figures 15 and 16, respectively.Figure 15
The real and estimated profiles of average reservoir pressure for the gas well under variable BHP and variable production rate conditions.Figure 16
The relation ofgpave to Gp for the gas well under variable BHP and variable production rate conditions.The slope of thegpave vs. Gp curve, msmb, is 1.663×10−3 MPa/(104 m3); hence, the gas in place is estimated to be 2.760×108 m3 with an error of 1.744%. It can be observed from Figures 15 and 16 that the value gap of b between the initial stage and the reference point is slightly larger than that between the reference point and medium late stage due to fluctuations of gas well production, changes in BHP, and long production period. Selecting the midsection data in the evaluation period as reference points may alleviate the fluctuations of b.The above three numerical examples have effectively proved the validity of the explicit method for estimating the average reservoir pressure and gas reserves based on the GPPF in the paper. Only one measured formation pressure is needed to give an understanding of the pressure drop history and gas in place based on the proposed method, which is widely applicable to various production systems as long as the BDF condition is satisfied. It is worth mentioning that the estimation error of the ARP and reserves may increase slightly if the measured pressure data are available merely at the beginning or end of production or if the formation pressure at the reference point is inferred by the initial pressure, production rate, and BHP data.


## 3.4. Field Case

The field case (well H-58) [51] is derived from Ibrahim et al. [52] who employed the superposition normalized pseudotime correlation to smooth field data and more accurately determine gas reserves. This gas well is from a tight gas reservoir with high pressure. The gas in place calculated from their normalized pseudotime plotting function gives 9.8 Bcf. Table 2 shows the reservoir properties and fluid parameters. The formation temperature Ti is 93.333°C, the original formation pressure pi is equal to 96.527 MPa, and pi/Zi is 55.240 MPa.Table 2
Reservoir and gas properties for the field case.
Property parametersProperty valuesϕi1.354×10−1Swci2.618×10−1h114.148 mpi96.527 MPaTi366.483 KTsc293.150 Kpsc1.013×105 PaM19.268 g/molTpc200 Kppc4.600 MPaρsc0.803 kg/m3Zsc0.998Zi1.747Cϕ8.702×10−4 MPa-1Cw5.802×10−4 MPa-1μi4.752×10−2 cpμsc1.055×10−2 cpCgi3.478×10−3 MPa-1Cti3.590×10−3 MPa-1Bgi2.299×10−3 m3/m3Figure17 illustrates the variations of fp and μZ with pressure p which indicates a possible large error caused by the pressure squared approximation and that μZ/p is only approximately constant at relatively high pressure. Figure 18 shows the integrand in pseudopressure calculated by the GPPF (fpn) and its error where the expression of fpn fitted by the least square method is given by
(46)fpnp≈2.364940×10−8p6−6.219404×10−6p5+4.933695×10−4p4+1.569904×10−3p3−2.069461p2+9.564624×101p−2.970716×101.Figure 17
The relation curves ofp/μZ and μZ vs. pressure (p) for the field case.Figure 18
The relation curves ofp/μZ and ε vs. p calculated by the polynomial function fpn for the field case.It can be seen from Figure18 that the polynomial function, fpn, can accurately represent the solution of fp=p/μZ within the pressure range of 2.3~96.5 MPa, and errors of most of the results do not exceed 1%. Figure 19 shows the changes in ppp and gp with p. The production rate and bottom hole flowing pressure data of well H-58 are exhibited in Figure 20.Figure 19
The relation curves of pseudopressure (pp) and gp vs. pressure (p) for the field case.Figure 20
Production history of the gas well H-58.Taking the production time of the 250th day as the reference point, it is inferred thatpk=31.089MPa,qk=22.724×104m3/d,andpwf,k=10.217 MPa. The estimated formation pressure profile is displayed in Figure 21. It can be observed that the initial production is not stable, and the estimated formation pressure does not drop, which may be caused by the mismatches of b between the initial points and the reference point. Part of the estimated pressure data during that period should be removed when one determines the gas reserves. Figure 22 shows the static material balance relationship of the gas well H-58 where the slope of the gpave vs. Gp relationship curve msmb equals 2.017×10−3 MPa/(104 m3), and subsequently, the gas reserves are estimated to be 2.739×108 m3 which is close to the estimate of Ibrahim et al. [52] (9.8Bcf=2.775×108 m3).Figure 21
The estimated profile of average reservoir pressure for the gas well H-58.Figure 22
The relation ofgpave to Gp for the gas well H-58.


## 4. Conclusions

The correlation between pressure and flow rate during BDF is derived in this paper using the theory of gas flow through porous media and the superposition principle. Based on the dynamic material balance equation (DMBE), a novel explicit method for estimating the ARP is proposed by using the gas property polynomial to approximate the nonlinear function,f=p/μZ. The conclusions are drawn as follows:
(i)
The differences among many definitions of pseudopressure for the linearization of the governing equation for gas flow are the integral calculation with respect tof or p/μZ. In view of the limitations of pressure and pressure squared, high-order polynomials can be employed to describe the nonlinear characteristics of the integrand f(ii)
Knowledge of the formation pressure at a certain point in the BDF stage can infer the ARP profile based on the DMBE and the GPPF. This method avoids the iteration of pseudotime and is widely applicable to various production constraints such as constant BHP, constant rate, and variable BHP/variable rate as long as the BDF condition holds true (that is,b can be approximated as a constant). It is recommended to use the data in the middle evaluation period as the reference point(iii)
Once the ARP is estimated, the gas reserves can also be determined by the SMBE. This method is simple and efficient and can be in conjunction with other methods to achieve mutual verification(iv)
The presented method takes into account the compressibilities of natural gas, rock pores, and irreducible water and is suitable for both abnormally pressured gas reservoirs and normal pressure gas reservoirs(v)
The presented method assumes the BDF condition after the transient flow for the gas well in a closed reservoir. The explicit PDA (including RTA and PTA) methods dispensing with pseudotime iteration and repeated curve fitting, however, remain investigated for determination of ARP for unconventional gas reservoirs subject to long periods of transient flow conditions



---

*Source: 4003656-2023-01-30.xml*