Abstract

Determination of gas in place (GIP) is among the hotspot issues in the field of oil/gas reservoir engineering. The conventional material balance method and other relevant approaches have found widespread application in estimating GIP of a gas reservoir or well-controlled gas reserves, but they are normally not cost-effective. To calculate GIP of abnormally pressured gas reservoirs economically and accurately, this paper deduces an iteration method for GIP estimation from production data, taking into consideration the pore shrinkage of reservoir rock and the volume expansion of irreducible water, and presents a strategy for selecting an initial iteration value of GIP. The approach, termed DMBM-APGR (dynamic material balance method for abnormally pressured gas reservoirs) here, is based on two equations: dynamic material balance equation and static material balance equation for overpressured gas reservoirs. The former delineates the relationship between the quasipressure at bottomhole pressure and the one at average reservoir pressure, and the latter reflects the relationship between average reservoir pressure and cumulative gas production, both of which are rigidly demonstrated in the paper using the basic theory of gas flow through porous media and material balance principle. The method proves effective with several numerical cases under various production schedules and a field case under a variable rate/variable pressure schedule, and the calculation error of GIP does not go beyond 5% provided that the production data are credible. DMBM-APGR goes for gas reservoirs with abnormally high pressure as well as those with normal pressure in virtue of its strict theoretical foundation, which not only considers the compressibilities of rock and bound water, but also reckons with the changes in production rate and variations of gas properties as functions of pressure. The method may serve as a valuable and reliable tool in determining gas reserves.

1. Introduction

The determination of gas in place, as a basic problem in the field of oil and gas reservoir engineering, is related to the development planning and production design of gas wells. Common methods to estimate reserves include volumetric method [1], material balance method [2], transient well test analysis [3, 4], and production decline analysis [57], where the first two approaches are usually applied to gas field or gas reservoir, while pressure transient analysis and rate transient analysis are mainly used to calculate well controlled reserves.

A variety of existing GIP determination methods based on the material balance equation, unfortunately, are inconvenient and uneconomical to implement owing to the demand for accurate average formation pressure data often obtained by shutting-in the wells or well testing. Therefore, many researchers turn to the investigation of production data-driven methods for calculating GIP or gas reserves.

For gas reservoirs without water drive, Mattar and McNeil (1995, 1998) [8, 9] proposed the flowing material balance procedure when the gas well is flowing at a constant rate. It is, however, less desirable to obtain convincing results because the variations of gas property parameters (such as gas viscosity and compressibility factor or Z-factor) with pressure or time are not considered.

Blasingame and Lee (1988) [10] presented the variable rate reservoir limits testing of gas wells, a new method of estimating GIP from production data (bottomhole pressures, flow rates, and cumulative gas production), based on the theory of gas flow derived from slightly compressible liquid solution and the introduction of material balance quasitime function tca. The method combines the gas flow equation (that is an approximation for variable rate and post-transient flow) with the material balance equation (MBE) of volumetric gas reservoirs, and requires a computer program to iterate on average reservoir pressure and material balance quasitime. Nevertheless, it does not consider the effects from characteristics of geopressured gas reservoirs nor explain how to select the initial value.

Mattar and Anderson (2005, 2006) [11, 12] developed the concept of variable rate flowing material balance or “dynamic material balance” extending the flowing material balance method to cases where the production rate is not constant and correspondingly flowing pressure is also variable. Like the treating process of Blasingame and Lee [10], the dynamic material balance procedure, derived from MBE of volumetric gas reservoirs and gas flow equation for boundary-dominated flow (or stabilized flow), also needs a calculation program to iterate on gas in place. On account of the use of pseudovariables (i.e., pseudopressure and pseudotime) accounting for the pressure-dependent changes in gas properties, Matter’s dynamic material balance procedure seemingly has a more rigorous theoretical basis than the previous flowing material balance, but it does not incorporate compressibilities of rock and irreducible water yet.

Material balance analysis of production performance for an overpressured gas reservoir is complicated because of the influence of rock and water compressibility in addition to significant gas compressibility effects. This makes it vital to include both formation and water compressibility in the production performance analysis when trying to estimate GIP of abnormally pressured reservoirs. In general, the traditional material balance equation or related methods can be directly used to calculate GIP, whereas well shutdown or well testing to determine the average reservoir pressure results in loss of production and expensive cost. Alternatively, those methods using production data [1317] almost do not reckon in the characteristic of abnormally pressured gas reservoirs, a fact that the elastic effects of rock pore and irreducible water are often not negligible as the formation pressure decreases [18, 19]. Therefore, this paper, taking the pore reduction in reservoir rock and the volume increase in bound water into account, develops an iterative calculation technique for determining gas in place of abnormally pressured gas reservoirs, an improvement of Matter’s dynamic material balance, based on a new total compressibility definition and the new material balance equation. The dynamic material balance equation (DMBE) is rigorously derived from the gas flow theory with variable rates; then, it is applied to GIP determination coupled with the material balance principle of overpressured gas reservoirs.

The layout of this article is as follows. First, the gas flow model for a closed gas reservoir with abnormally high pressure is introduced and solved in order to obtain the dynamic material balance equation. Second, the corresponding static material balance equation is presented which is then coupled with DMBE to compute average reservoir pressure. Third, the iterative method (DMBM-APGR) for GIP estimation is developed and the detailed iteration steps are described. Fourth, several case studies illustrate and validate the method, and a brief discussion reveals the novelties of the work. Finally, conclusions are drawn.

2. Mathematical Model of Gas Flow through Porous Media

When it comes to the problem of liquid flow through porous media, the comprehensive differential equation (governing equation) is a typical diffusivity equation where the diffusivity coefficient can be approximately seen as a constant; thus, it is easy to solve; but for gas flow, the situation is quite different in that the linearization of the corresponding diffusivity equation is essential to consider the pressure-dependent changes in gas properties such as viscosity, Z-factor, and compressibility.

The key to solve the gas flow problem is the linearization of the diffusivity equation. Al-Hussainy et al. (1966) [20, 21] and Russell et al. (1966) [22] took the lead in scoring representative research results in this field. The concept of “pseudopressure” proposed by them can partially linearize the gas flow diffusivity equation. The pseudopressure presented by Russell et al. (1966) is, in fact, a normalized pseudopressure with the dimension of pressure.

Agarwal (1979) [23] proposed the concept of pseudotime, representing a new stage for solution of gas flow. Lee and Holditch (1982) [24] pointed out that the application of pseudopressure and pseudotime can effectively linearize the gas diffusivity equation, and it is feasible to apply the slightly compressible liquid solution to gas flow.

Fraim and Wattenbarger (1987) [25], introducing the real gas pseudopressure and a normalized time function, investigated the gas flow problem for a closed gas reservoir producing against constant wellbore pressure during boundary-dominated flow, where the normalized time should adopt the viscosity and compressibility at the average reservoir pressure rather than the bottomhole pressure.

The concepts of quasipressure and quasitime put forward by Meunier et al. (1984, 1987) [26, 27] retain the dimensions of real pressure and real time, and are also known as normalized pseudopressure and normalized pseudotime, respectively, more convenient to use than before.

In addition, some researchers [28, 29] employ the time-dependent dimensionless viscosity-compressibility ratio λ, the time-averaged λ (i.e., λ¯), and gas density to approximately and successfully linearize the density-diffusivity equation for unsteady state analysis of natural gas reservoirs. This density-based approach, in fact, can be expressed as the counterpart mathematically equivalent to the pseudofunction-based approach for the gas flow diffusivity equation.

In this section, firstly, the gas flow model including water, rock, and gas compressibility effects is established on the basis of the mass conservation principle and state equations. Secondly, the constant rate solution for the gas flow model is obtained using Laplace transform and liquid-based analytical solutions rewritten in terms of quasivariables. Then, the variable rate solution is also acquired according to the principle of quasipressure drop superposition. Finally, the dynamic material balance equation for abnormally pressured gas reservoirs is proved based on the definition of average reservoir quasipressure and the preceding variable rate solution.

2.1. Gas Flow Governing Equation

The principle of continuity for gas flow through porous media in its differential form can be expressed as follows:
(1)ρgvxx+ρgvyy+ρgvzz+qg=ρg1Swcϕt,
where ρg (kg/m3) denotes the gas density, vx,vy,vz represent the component of velocity in the x, y, z directions, respectively, ϕ is the porosity, qg [kg/(m3·s)] denotes the change in fluid mass caused by the source or sink in the unit volume per unit time, Swc denotes the irreducible water saturation, and t (s) is the time.
The continuity equation of single-phase gas flow without source/sink term is as follows:
(2)ρgvg=1Swcρgϕt+ρgϕ1Swct.
Without regard to the non-Darcy effect by high flow velocity, or the additional threshold pressure gradient by the presence of adsorption or hydration films on the rock surface, the gas flow is subject to Darcy’s law:
(3)vg=Kμp,
where vg (m/s) is the gas velocity vector, μ (Pa·s) denotes the gas viscosity, K (m2) denotes the effective permeability, and p (Pa) represents the pressure.
According to the definition of rock compressibility, we obtain the following:
(4)Cϕ=1VpdVpdp=1ϕdϕdp,
where Cϕ (Pa-1) denotes the rock compressibility and Vp (m3) represents the pore volume.
Integrating above Equation (4) with respect to p yields porosity and pore volume given by the following:
(5)ϕ=ϕieCϕppi,Vp=VpieCϕppi,
where pi (Pa) denotes the initial reservoir pressure; ϕi and Vpi (m3) represents the porosity and pore volume under initial reservoir condition, respectively.
Isothermal compressibility of natural gas and gas density can be written as follows:
(6)Cg=ZTpddppZT=ZpddppZ,(7)ρg=pMZRT,
where Cg (Pa-1) denotes the gas isothermal compressibility, Z represents the deviation factor (Z-factor), T (K) is the temperature, M (kg/mol) denotes the gas molar mass, and R [8.3144598 J/(mol·k)] denotes the molar gas constant [30, 31].
According to the definition of water compressibility, there is
(8)Cw=1VwdVwdp.
By separating variables and integrating Equation (8) with respect to p, one obtains the following:
(9)Vw=VwieCwppi.
Then, the bound water saturation is given by the following:
(10)Swc=VwVp=VwieCwppiVpieCϕppi=SwcieCw+Cϕppi,
where Vw (m3) denotes the volume of bound water, Vwi represents the initial volume of bound water, Cw (Pa-1) is the water compressibility, and Swci is the initial bound water saturation.
Substituting Equations (3), (5), (6), (7), and (10) into Equation (2), the rearranged equation is represented by the following:
(11)pμZp=ϕiμCtKpμZpt,(12)Ct=eCϕppiCϕ+1SwcCg+SwcCw.
Equation (11) is the governing equation of unsteady single-phase gas flow in which gas property parameters (μ, Z, and Ct) are not constant and change with pressure. Thus, it is difficult to solve it directly. Referring to the concepts of pseudopressure proposed by Russell et al. (1966) [22] and quasitime developed by Meunier et al. (1984, 1987) [26, 27], quasivariables in this paper are defined by the following:
(13)pp=pi+μiρgipipρgξμξdξ=pi+μiZipipipξμξZξdξ,(14)ta=μiCti0t1μCtdt,
where pp (Pa) denotes the quasipressure function, ta (s) denotes the quasitime function, and Ct (Pa-1) denotes the total compressibility function defined by Equation (12); μi, Zi, ρgi, and Cti represent the gas viscosity, deviation factor, density, and total compressibility under initial reservoir condition, respectively.
Based on Equations (13) and (14), then Equation (11) can be transformed into the following:
(15)2pp=ϕiμiCtiKppta.

Equation (15) is characterized by the similar form of the diffusivity equation for liquid flow. It could be seen that the form of expression for gas flow behaviors is consistent with that of liquid based on the definitions of quasivariables in the paper, so the slightly compressible liquid solution could be applied to gas flow governing equation.

2.2. Constant Rate Solution

The mathematical model of gas flow for a single well producing in a bounded reservoir with a constant production rate can be written as follows:
(16)2ppr2+1rppr=ϕiμiCtiαtKppta,rpprr=rw=qμiBgiαpKh,pprr=re=0,ppta=0=pppi,
wherer (m) is the distance from a location to the center of the gas well, rw (m) is the wellbore radius, re (m) denotes the reservoir limit, Bgi (m3/m3) represents the gas formation volume factor under original formation condition, and h (m) denotes the net pay thickness; αp denotes the conversion factor of dimensionless pressure, and αt is the conversion factor of dimensionless time. Under the international system of units, αp=2π and αt=1; if other unit systems are used, their values are different. For the sake of simplicity, all the symbols in this article are explained in international standard units.
Dimensionless variables for quasipressure, time, and radius are defined by the following:
(17)pD=αpKhppippqμiBgi,(18)tD=αtKϕiμiCtirw2ta,(19)rD=rrw,reD=rerw,
where pD, tD, and rD are the dimensionless quasipressure, dimensionless quasitime, and dimensionless radius, respectively; reD denotes the dimensionless radial length at the reservoir limit.
Substituting Equations (17)–(19) into Equation (16) yields the dimensionless mathematical model of gas flow given by the following:
(20)2pDrD2+1rDpDrD=pDtD,rDpDrDrD=1=1,pDrDrD=reD=0,pDtD=0=0.

From the definitions of quasivariables, it can be easily noted that quasitime ta is a function of time t, while the gas viscosity μ and total compressibility Ct are the functions of pressure p dependent on location and time; so ta is a function of location rand time t. Similarly, quasipressure pp is only related to pressure; therefore, it is also a function of r and t. Strictly speaking, pp and ta correspond to the same pressure p, that is, the quasipressure and the quasitime correspond to the same position and the same time, for which we have a fat chance to solve Equations (16) and (20).

The general approach to address the above difficulties, nowadays, is to calculate the value of ta under a certain reference pressure irrespective of location r, i.e.,
(21)ta=μiCti0t1μCtdtμiCti0t1μprefCtprefdt,
where pref represents a reference pressure independent of location.
In this way, ta is approximately reckoned to be a function merely dependent on time t, so tD is only related to time, too. The Laplace transform on tD in Equation (20) enables us to gain the dimensionless pressure drop pD which can be expressed as a solution in Laplace domain:
(22)p~D=K1reDsI0rDs+I1reDsK0rDsssK1sI1reDsK1reDsI1s,
where p~D denotes the variable in Laplace domain that corresponds to pD and s is the variable in Laplace domain that corresponds to tD; I0 and I1 are the zero-order and first-order modified Bessel functions of the first kind, respectively; K0 and K1 are the zero-order and first-order modified Bessel functions of the second kind, respectively.
The analytical inversion expression of Equation (22) in the real time domain is as follows:
(23)pD=2tDreD21reD2reD21lnrD+4reD4lnreD3reD4+2rD2reD21+2reD2+14reD212πn=1eλn2tDJ12reDλnY1λnJ0rDλnJ1λnY0rDλnλnJ12reDλnJ12λn,
where J0 and J1 are the zero-order and first-order Bessel functions of the first kind, respectively; Y0 and Y1 are the zero-order and first-order Bessel functions of the second kind (Neumann functions), respectively. The variable λ=λn must satisfy the following condition:
(24)J1reDλY1λY1reDλJ1λ=0.

2.3. Variable Rate Solution

As shown in Figure 1, there are m times of production fluctuations in the production process of the gas well. According to the principle of pressure drop superposition, the quasipressure drop caused by the changes in production rate can be written as follows:
(25)ppippr,ta=i=1mΔppi=i=1mμiBgiαpKhqiqi1pDtata,i1,q0=0,ta,0=0,
where Δppi represents the quasipressure drop caused by ith production fluctuation, qi is the production rate during ith production period, qi1 denotes the production rate during (i-1)th production period, and ta,i1 is the quasitime for the start of ith production period.
Figure 1

Schematic diagram of variable production rates of a gas well.

Figure 1

Schematic diagram of variable production rates of a gas well.

Combining the definition of dimensionless time Equation (18) with Equations (23) and (25) gives the following:
(26)ppipp=μiBgiαpKh2re2rw2αtKϕiμiCtiGpa+re2re2rw2lnrrw+4re4lnre/rw3re4+2r2re2rw2+2re2rw2+rw44re2rw22qmi=1mqiqi1n=1πλneλn2αtK/ϕiμiCtirw2tata,i1J12rerwλnY1λnJ0r/rwλnJ1λnY0r/rwλnJ12re/rwλnJ12λn,(27)Gpa=i=1mqiqi1tata,i1=i=1mqita,ita,i1=0taqdta.
By dividing both ends of Equation (26) by qm and setting rD=1, i.e., r=rw, the downhole quasipressure drop correlation is obtained:
(28)ppippwfqm=μiBgiαpKh2re2rw2αtKϕiμiCtitca+4re4lnre/rw3re4+4re2rw2rw44re2rw22i=1mqiqi1qmn=1πλneλn2αtK/ϕiμiCtirw2tata,i1J12rerwλnY1λnJ0λnJ1λnY0λnJ12re/rwλnJ12λn,(29)tca=Gpaqm=1qm0taqdta=1qm0tqμiCtiμprefCtprefdt,
where tca denotes the material-balance quasitime, Gpa represents a variable similar to cumulative gas production defined by Equation (27), ppwf is the quasipressure corresponding to the bottomhole pressure pwf, and ppi is the quasipressure corresponding to initial reservoir pressure (that is equal to pi).
For the convenience of writing, the following expression is introduced here:
(30)fta,λn=eλn2αtK/ϕiμiCtirw2tata,i1J12re/rwλnJ12re/rwλnJ12λn.
Due to the identical relation [32] given by
(31)Y0xJ1xY1xJ0x=2π1x,
Equation (28) can be transformed into the following:
(32)ppippwfqm=μiBgiαpKh2re2rw2αtKϕiμiCtitca+μiBgiαpKh4re4lnre/rw3re4+4re2rw2rw44re2rw22μiBgiαpKhi=1mqiqi1qmn=1πλnfta,λn2π1λn.
Then, Equation (32) can be rewritten as follows:
(33)ppippwfqm=2παtαpBgiπre2rw2hϕiCtitca+μiBgiαpKh4re4lnre/rw3re4+4re2rw2rw44re2rw22+μiBgiαpKhi=1mqiqi1qmn=12λn2fta,λn.
The symbols mbdf and b are used to denote the slope and intercept of the straight line in Δppiwf/q vs. tca curve, respectively, i.e.,
(34)ppippwfq=Δppiwfq=mbdftca+b,(35)mbdf=2παtαpBgiAhϕiCti,(36)b=μiBgiαpKh4re4lnre/rw3re4+4re2rw2rw44re2rw22+i=1mqiqi1qmn=12λn2fta,λn,(37)A=πre2rw2,
where Δppiwf represents the difference between the quasipressure corresponding to original reservoir pressure pi and that corresponding to bottomhole pressure pwf, and A is the gas reservoir area.
Since transient and stabilized flow [33] may alternate for each production rate change, it can be suggested that stabilized flow will dominate for a small rate change. According to Equation (35), the relationship between gas in place (or well controlled reserves) G and the slope mbdf can be expressed as follows:
(38)G=Ahϕi1SwciBgi=1mbdf2παtαp1SwciCti.

It is worth noting that once the Δppiwf/q vs. tca graph is drawn and the slope mbdf of the straight line segment is determined, one can calculate GIP based on this slope, which is also the basis for selecting the initial iteration value in the paper. Since Δppiwf/q vs. tca curve cannot be obtained in advance, the authors innovatively use the slope of pipwf/q vs. t curve to replace the initial value of mbdf, and GIP estimated by this slope is taken as the initial value of iteration G0.

It should also be kept in mind that the quasitime ta, in practice, implies a reference pressure independent of location; thus, the material-balance quasitime function tca is also supposed to be determined at the same reference pressure. At first, Agarwal (1979) [23] and Lee and Holditch (1982) [24] took the bottomhole flowing pressure as the reference pressure, but currently, the reasonable procedure is to adopt the average reservoir pressure [10, 25, 28, 3436] instead, that is,
(39)tca=μiCtiqt0tqtμpaveCtpavedt,
where pave denotes the average reservoir pressure.

2.4. Dynamic Material Balance Equation for Abnormally Pressured Gas Reservoirs

By differentiating pp in Equation (26) with respect to ta, one obtains the following:
(40)ppta=μiBgiαpKh2re2rw2αtKϕiμiCtiqmμiBgiαpKh×i=1mqiqi1n=1πλnfta,λnY1λnJ0rrwλnJ1λnY0rrwλnλn2αtKϕiμiCtirw2.
Then, the governing equation of gas flow can be written as follows:
(41)1rrrppr=ϕiμiCtiαtKppta=qmμiBgiαpKh2re2rw2μiBgiαpKh×i=1mqiqi1n=1πλnrw2fta,λnY1λnJ0rrwλnJ1λnY0rrwλn.
Considering
(42)rJ0ardr=rJ1ara,rY0ardr=rY1ara,
integrating Equation (41) with respect to r gives the following:
(43)rppr=qmμiBgiαpKhr2re2rw2+C1μiBgiαpKh×i=1mqiqi1n=1πrwfta,λnY1λnrJ1λnrwrJ1λnrY1λnrwr.
Equation (43) could be rewritten as follows:
(44)ppr=qmμiBgiαpKhrre2rw2+C1rμiBgiαpKh×i=1mqiqi1n=1πrwfta,λnY1λnJ1λnrwrJ1λnY1λnrwr.
The external boundary condition coincides with the following formula:
(45)pprr=re=qmμiBgiαpKhrere2rw2+C1reμiBgiαpKh×i=1mqiqi1n=1πrwfta,λnY1λnJ1λnrwreJ1λnY1λnrwre=0.
From Equation (45), it follows the following equation:
(46)C1=re2re2rw2.
Owing to
(47)J1ardr=J0ara,Y1ardr=Y0ara,
by integrating Equation (44) with respect to r, one obtains the following:
(48)pp=qmμiBgiαpKh1re2rw2r22+C1lnr+C2μiBgiαpKh×i=1mqiqi1n=1πλnfta,λnY1λnJ0λnrwr+J1λnY0λnrwr.
Setting r=rw, the equation above can be transformed into the following:
(49)ppwf=qmμiBgiαpKh1re2rw2rw22+C1lnrw+C2μiBgiαpKh×i=1mqiqi1n=1πλnfta,λnY1λnJ0λn+J1λnY0λn.
Subtracting Equation (49) from Equation (48) and combining it with Equations (31) and (46), we obtain the following:
(50)ppppwf=qmμiBgiαpKhre2rw2re2lnrrwr2rw22μiBgiαpKh×i=1mqiqi1n=1fta,λnπλnY1λnJ0λnrwr+πλnJ1λnY0λnrwr2λn2.
Based on Equation (50), the quasipressure pp can be rewritten as follows:
(51)pp=Ar+Br,(52)Ar=ppwf+qmμiBgiαpKhre2rw2re2lnrrwr2rw22,(53)Br=μiBgiαpKhi=1mqiqi1n=1fta,λnπλnY1λnJ0λnrwr+πλnJ1λnY0λnrwr2λn2.
The definition of average reservoir quasipressure is given by the following:
(54)ppave=1πre2rw2rwrepp2πrdr.
By substituting Equation (51) into Equation (54), one obtains the following:
(55)ppave=I1+I2=2re2rw2rwreArrdr+2re2rw2rwreBrrdr,
where I1 and I2 are the integral formulas related to Ar and Br, respectively.
According to Equation (52), I1 can be expressed as follows:
(56)I1=2re2rw2rwreppwf+qmμiBgiαpKhre2rw2re2lnrrwr2rw22rdr=ppwf+2re2rw2qmμiBgiαpKhre2rw2rwrere2lnrrwr2rw22rdr.
Due to
(57)rwrere2lnrrwr2rw22rdr=re2rw22+re42lnrerw38re4rw48,
Equation (56) can be transformed into the following:
(58)I1=ppwf+qmμiBgiαpKh4re4lnre/rw3re4+4re2rw2rw44re2rw22.
According to Equation (53), analogously, I2 can be expressed as follows:
(59)I2=2re2rw2rwreBrrdr=μiBgiαpKh2re2rw2i=1mqiqi1n=1fta,λnrwreπλnY1λnJ0λnrwr+πλnJ1λnY0λnrwr2λn2rdr.
Then, the equation above can be transformed into the following:
(60)I2=μiBgiαpKh2re2rw2i=1mqiqi1n=1fta,λnre2rw2λn2=μiBgiαpKhi=1mqiqi1n=1fta,λn2λn2.
Substituting Equations (58) and (60) into Equation (55) yields the following:
(61)ppave=ppwf+qmμiBgiαpKh4re4lnre/rw3re4+4re2rw2rw44re2rw22+μiBgiαpKhi=1mqiqi1n=1fta,λn2λn2.
Assume ppaveppave=pppave, then divide both ends of Equation (61) by qm, and eventually acquire the following expression:
(62)ppaveppwfqm=μiBgiαpKh4re4lnre/rw3re4+4re2rw2rw44re2rw22+μiBgiαpKhi=1mqiqi1qmn=1fta,λn2λn2,
where ppave or pppave denotes the quasipressure corresponding to average reservoir pressure.
According to Equation (36), Equation (62) can be rewritten as follows:
(63)ppaveppwfqm=b,
or
(64)ppave=ppwf+qb.

The above derivation process has cogently proved the dynamic material balance equation for abnormally pressured gas reservoirs, i.e., Equation (64), which indicates that one can infer what the average reservoir quasipressure should be from the downhole quasipressure, indirectly capturing the intrinsic connection between the average formation pressure and the flowing bottomhole pressure. The dynamic material balance equation here is almost consistent with Matter’s variable rate flowing material balance equation derived from gas flow equation for boundary-dominated flow, namely, stabilized flow, and MBE of volumetric gas reservoirs; however, our derivation procedure embodies the pore shrinkage of reservoir rock and the volume expansion of irreducible water based on a new definition of total compressibility and variable rate gas flow theory.

According to Equation (39), it is necessary to comprehend the average reservoir pressure data over time in order to calculate tca, but pave values are often difficult to measure. Therefore, the iterative method, where an initial GIP value is postulated to calculate the pave and tca data and then GIP is updated by these data time and again, can be employed to determine the gas reserves. The average reservoir pressure pave could be determined by the (static) material balance equation of abnormally pressured gas reservoirs which will be introduced below.

3. Material Balance Equation for Abnormally Pressured Gas Reservoirs

For an overpressured gas reservoir, as a closed system, it is reasonable to assume that the invasion of edge or bottom water is negligible. Its material balance principle, therefore, can be phrased that natural gas production is mainly caused by the expansion of gas volume, bound water expansion, and the reduction in rock pore volume:
(65)GpBg=ΔVw+ΔVp+GBgBgi,
where Gp (m3) denotes the cumulative gas production, G (m3) represents the gas in place or gas reserves, Bg (m3/m3) represents the gas formation volume factor, Bgi (m3/m3) is the gas formation volume factor under original reservoir condition, ΔVp (m3) is the reduced pore volume, and ΔVw (m3) denotes the expanded bound water volume.
From Equations (5) and (9), one obtains the following:
(66)ΔVp=VpiVp=Vpi1eCϕpavepi,(67)ΔVw=VwVwi=VwieCwpavepi1.
Pore volume and bound water volume under initial reservoir condition can be expressed, respectively, as follows:
(68)Vpi=GBgi1Swci,Vwi=GBgi1SwciSwci.
Gas formation volume factor is given by the following:
(69)Bgpave=ZTpavepscZscTsc,Bgi=ZiTipipscZscTsc,
where psc (1.01325×105Pa), Tsc (293.15 K), and Zsc denote the pressure, temperature, and deviation factor (Z-factor) under the standard condition, respectively; Ti is the initial reservoir temperature.
It is deemed justifiable that the temperature of the gas reservoir during its depletion nearly remains unchanged. Substituting Equations (66)–(69) into Equation (65) yields the following:
(70)gpave=piZi1GpG,(71)gpave=paveZpaveeCϕpavepiSwcieCwpavepi1Swci,
where gpave denotes a function related to average reservoir pressure, defined by Equation (71).

The material balance equation for abnormally pressured gas reservoirs, i.e., Equation (70), expresses a relationship between the pave-dependent pressure function gpave in the overpressured gas reservoir and the cumulative gas production Gp, from which it can be observed that a straight-line section should appear on the gpave vs. Gp relationship curve with its intercept pi/Zi and the slope mmbe or pi/ZiG, helpful to determine the value of G. If the elastic effect of rock and bound water is insignificant, that is Cw=Cϕ0, then gpave can convert into pave/Zpave, which is applicable to volumetric gas reservoirs. For overpressured gas reservoirs, nevertheless, the function gpave is supposed to be used; that is to say, it seems unreasonable to rashly apply pave/Zpave~Gp relation to G estimation.

4. Dynamic Material Balance Method for Abnormally Pressured Gas Reservoirs

The relationships between Equations (34), (64), and (70) prompt us to calculate gas in place G, average reservoir pressure pave, and material balance quasitime tca by the iterative method shown in Figure 2. Its detailed steps are as follows:

  • (1)

    Collect and sort out the known data available, such as gas production rate q, cumulative gas production Gp, and bottomhole pressure pwf

  • (2)

    Use the slope of pipwf/q vs. t curve as mbdf and estimate the initial GIP iteration value G0 via Equation (38) (see Section 5.1.1 in the paper for details)

  • (3)

    Calculate the function gpave by Equation (70) and then determine pave according to Equation (71)

  • (4)

    Calculate tca by Equation (39), and total compressibility should be computed by Equation (12)

  • (5)

    Determine the intercept b of Δppiwf/q vs. tca curve according to Equation (34)

  • (6)

    Estimate pppave by Equation (64) and then convert pppave into pave according to Equation (13)

  • (7)

    Calculate gpave by Equation (71), draw gpave vs. Gp relationship curve, and set the intercept of the straight line as pi/Zi, then use the slope mmbe to determine the gas in place Ge

  • (8)

    Compare Ge calculated in Step 7 with the initial value G0 in Step 2. If their difference meets accuracy requirement, then the estimated gas in place is Ge; otherwise, assign Ge to G0 and repeat Steps 2 to 8

Figure 2

The iteration steps of dynamic material balance method for abnormally pressured gas reservoirs.

Figure 2

The iteration steps of dynamic material balance method for abnormally pressured gas reservoirs.

Owing to the application of dynamic material balance equation, i.e., Equation (64) and (static) material balance equation, i.e., Equation (70) to overpressured gas reservoirs, the above method for estimating gas in place is called DMBM-APGR (dynamic material balance method for abnormally pressured gas reservoirs) in the paper. It is inessential to shut the gas well or conduct a well testing destined for static pressure since the method only needs production data (q, pwf, Gp, etc.) available to determine gas in place G, which boasts low cost and considerable convenience in that the computing process can be easily completed by programs and data sheets.

5. Results and Discussion

5.1. Numerical Simulation Verification

Several numerical examples are listed below to illustrate the applicability of DMBM-APGR to overpressured gas reservoirs. We employ the black oil model to simulate the radial single-phase gas flow with a production well located in the center of a circular closed formation, considering the shrinkage of rock pore volume and the expansion of irreducible water volume. The radial grid type is 200×72×1; the original reservoir pressure pi is 55 MPa, the temperature Ti is 108°C, the initial rock porosity ϕi is 0.2, and the initial irreducible water saturation Swci is 0.2. The depths of the top face of each grid block are set to 2995 m. The pore volume at initial reservoir condition Vpi and actual gas in place G obtained by initialization are 565864 m3 and 1.51906310×108m3, respectively. Other parameters about gas reservoir and natural gas are shown in Table 1.

Table 1

Simulated reservoir and gas properties.

ParametersProperty valuesParametersProperty values
dr
1.5 m
ϕ
0.2
dθ
Swci
0.2
dz
10 m
Mg
19.576754 g/mol
Kr
3 mD
Zsc
0.997514
Kθ
3 mD
Cϕ
2.364×103MPa1
Kz
0.03 mD
Cw
4.230×104MPa1
pi
55 MPa
μw
0.264565 cp
Ti
381.15 K
μi
3.401434212×102cp
Tpc
201.325729 K
Cgi
8.924326991×103MPa1
ppc
4.597109 MPa
Cti
9.588061593×103MPa1
re
300 m
Bgi
2.980066777×103m3/m3
rw
0.1 m
Vpi
565864 m3
h
10 m
G
1.51906310×108m3
ParametersProperty valuesParametersProperty values
dr
1.5 m
ϕ
0.2
dθ
Swci
0.2
dz
10 m
Mg
19.576754 g/mol
Kr
3 mD
Zsc
0.997514
Kθ
3 mD
Cϕ
2.364×103MPa1
Kz
0.03 mD
Cw
4.230×104MPa1
pi
55 MPa
μw
0.264565 cp
Ti
381.15 K
μi
3.401434212×102cp
Tpc
201.325729 K
Cgi
8.924326991×103MPa1
ppc
4.597109 MPa
Cti
9.588061593×103MPa1
re
300 m
Bgi
2.980066777×103m3/m3
rw
0.1 m
Vpi
565864 m3
h
10 m
G
1.51906310×108m3

In Table 1, dr denotes the radial length of the grid, dθ denotes the angular size of each grid block, and dz denotes the grid height; Kr, Kθ, and Kz are the radial permeability, azimuthal permeability, and vertical permeability, respectively; Tpc represents the pseudocritical temperature of gas mixture and ppc represents the gas pseudocritical pressure.

The viscosity correlations proposed by Londono et al. (2002, 2005) [37, 38] are used to calculate the gas viscosity, and the Hall-Yarborough (1973) method [39] is used to determine the Z-factor and gas compressibility Cg. The resulting relation between gas formation volume factor Bg and pressure p and that between gas viscosity μ and p are shown in Figure 3. In line with Equation (13), Figure 4 reveals the relationship between quasipressure function pp and pressure p.

Figure 3

Gas formation volume factor (Bg) and gas viscosity (μ) vs. pressure (p) for numerical cases.

Figure 3

Gas formation volume factor (Bg) and gas viscosity (μ) vs. pressure (p) for numerical cases.

Figure 4

Gas quasipressure function (pp) vs. pressure (p) for numerical cases.

Figure 4

Gas quasipressure function (pp) vs. pressure (p) for numerical cases.

5.1.1. Constant Production Rate Simulation

We set 1111-time steps, and the simulation period lasts about 3 years where the time intervals represented by time steps 1~12 are one hour, time steps 13 and 14 are 3 h and 9 h, respectively, and the time intervals of time steps 15 and 16 are 0.5 d. Other time steps’ time intervals are all one day. The production rate of the gas well, a constant, was set to 1.6×104m3/d. The production schedules are shown in Table 2 where TS denotes time step and Δt represents time interval.

Table 2

Production schedules for the constant rate case.

t (d)TSΔt (h)Δt(d)q (104 m3/d)
1~10971~1210.0416671.6
1330.1251.6
1490.3751.6
15~16120.51.6
17~11112411.6
t (d)TSΔt (h)Δt(d)q (104 m3/d)
1~10971~1210.0416671.6
1330.1251.6
1490.3751.6
15~16120.51.6
17~11112411.6

The time-dependent data such as production rate q, bottomhole flowing pressurepwf, and cumulative gas production Gp, obtained by numerical simulation, can be readily recorded. Since the relationships between gas viscosity, deviation factor, and pressure are known, the quasipressure data corresponding to bottomhole pressure can be calculated according to Equation (13). Therefore, the ordinate plotting function Δppiwf/q in Equation (34) can be also determined. Now, we only need to iteratively calculate the values of the abscissa plotting function tca.

To determine the material-balance quasitime tca requires a given GIP value G0. This initial value of iteration can be determined by the slope of pipwf/q vs. tcurve, that is, this slope can be used to replace mbdf in Equation (38). As shown in Figure 5, when the quasipressure and quasitime are not used, the slope of pipwf/q vs. t curve is 5.092941367×103MPa·104m31, so the GIP estimated by the slope is 163828893 m3 according to Equation (39), which can be used as the initial iteration value. The relevant results by DMBM-APGR are shown in Table 3.

Figure 5

pipwf·q1 vs. t from real pressure (p) and real time (t) data for the constant flow rate case.

Figure 5

pipwf·q1 vs. t from real pressure (p) and real time (t) data for the constant flow rate case.

Table 3

Estimation of GIP by DMBM-APGR with G0=1.63828893×108m3 for the constant flow rate case.

No.G0 (m3)b [×10 MPa/(104 m3/d)]mmbe [×10-3 MPa/(104 m3)]Ge (m3)Ea (%)
11638288931.0848810822.8297604121566130293.098
21566130291.0846652172.8314328511565205223.038
31565205221.0846623072.8314555181565192693.037
41565192691.0846622682.8314557331565192583.037
51565192581.0846622682.8314557331565192583.037
No.G0 (m3)b [×10 MPa/(104 m3/d)]mmbe [×10-3 MPa/(104 m3)]Ge (m3)Ea (%)
11638288931.0848810822.8297604121566130293.098
21566130291.0846652172.8314328511565205223.038
31565205221.0846623072.8314555181565192693.037
41565192691.0846622682.8314557331565192583.037
51565192581.0846622682.8314557331565192583.037
In Table 3, Ge is the GIP determined by G0 during the same iteration step, and Ea is the estimation error of GIP defined by the following equation:
(72)Ea=GeGG×100%.

If the real pressure and real time are used directly, as mentioned above, the GIP value calculated by the slope of the straight line in pipwf/q vs. t graph is 163828893 m3 which is 7.849% larger than the real GIP. The GIP values estimated by DMBM-APGR, however, gradually stabilize to 156519258 m3, and the error is merely 3.037% satisfactory for reservoir engineering. The difference between two adjacent estimated values after two iterations, in fact, is small enough to terminate calculation. Though the initial iteration values of GIP are far away from the real value G, the estimation by DMBM-APGR will tend to be the same result (see Tables 4 and 5).

Table 4

Estimation of gas in place by DMBM-APGR with G0=0.3×108m3 for the constant flow rate case.

No.G0 (m3)b [×10 MPa/(104 m3/d)]mmbe [×10-3 MPa/(104 m3)]Ge (m3)Ea (%)
1300000001.0896010432.7932074471586625264.448
21586625261.0847287302.8309408571565477243.055
31565477241.0846631632.8314489871565196313.037
41565196311.0846622792.8314556661565192613.037
51565192611.0846622682.8314557331565192583.037
No.G0 (m3)b [×10 MPa/(104 m3/d)]mmbe [×10-3 MPa/(104 m3)]Ge (m3)Ea (%)
1300000001.0896010432.7932074471586625264.448
21586625261.0847287302.8309408571565477243.055
31565477241.0846631632.8314489871565196313.037
41565196311.0846622792.8314556661565192613.037
51565192611.0846622682.8314557331565192583.037
Table 5

Estimation of gas in place by DMBM-APGR with G0=7.5×108m3 for the constant flow rate case.

No.G0 (m3)b [×10 MPa/(104 m3/d)]mmbe [×10-3 MPa/(104 m3)]Ge (m3)Ea (%)
17500000001.0882358832.8037758351580644734.054
21580644731.0847103842.8310829291565398683.050
31565398681.0846629162.8314506811565195373.037
41565195371.0846622762.8314556781565192613.037
51565192611.0846622682.8314557331565192583.037
No.G0 (m3)b [×10 MPa/(104 m3/d)]mmbe [×10-3 MPa/(104 m3)]Ge (m3)Ea (%)
17500000001.0882358832.8037758351580644734.054
21580644731.0847103842.8310829291565398683.050
31565398681.0846629162.8314506811565195373.037
41565195371.0846622762.8314556781565192613.037
51565192611.0846622682.8314557331565192583.037

5.1.2. Constant Bottomhole Pressure Simulation

Setting the production schedule of gas well as pwf=30MPa without changes in other parameters generates the constant pressure production case as shown in Table 6. According to the simulation results, the slope of pipwf/q vs. t curve is 6.582477939×103MPa·(104 m3)-1 as shown in Figure 6, so the GIP estimation from this slope is 126756360m3 with an unacceptable error of -16.556%. The estimation, nonetheless, can serve as the initial value of iteration, and the resulting calculation results by DMBM-APGR are shown in Table 7.

Table 6

Production schedules for the constant pressure case.

t (d)TSΔt(h)Δt (d)pwf (MPa)
1~10971~1210.04166730
1330.12530
1490.37530
15~16120.530
17~111124130
t (d)TSΔt(h)Δt (d)pwf (MPa)
1~10971~1210.04166730
1330.12530
1490.37530
15~16120.530
17~111124130
Figure 6

pipwf·q1 vs. t for the constant pressure case.

Figure 6

pipwf·q1 vs. t for the constant pressure case.

Table 7

Estimation of GIP by DMBM-APGR with G0=126756360m3 for the constant pressure case.

No.G0 (m3)b×10MPa/104m3/dmmbe×103MPa/104m3Ge (m3)Ea (%)
11267563601.0771189882.8765287411540667201.422
21540667201.0784964862.8669987821545788411.759
31545788411.0785171202.8668559901545865401.764
41545865401.0785174292.8668539361545866511.764
51545866511.0785174342.8668539101545866531.764
61545866531.0785174342.8668539101545866531.764
No.G0 (m3)b×10MPa/104m3/dmmbe×103MPa/104m3Ge (m3)Ea (%)
11267563601.0771189882.8765287411540667201.422
21540667201.0784964862.8669987821545788411.759
31545788411.0785171202.8668559901545865401.764
41545865401.0785174292.8668539361545866511.764
51545866511.0785174342.8668539101545866531.764
61545866531.0785174342.8668539101545866531.764

Table 7 indicates that it only takes two iterations for DMBM-APGR to achieve the expected accuracy in that the difference between the GIP estimation in 2nd iteration and that in 1st iteration is merely 0.332%, and then the calculation could be terminated. If the iterative process continues, the ultimate results will become 154586653 m3 with an error of 1.764%. The calculation results, when selecting other values as the initial value of iteration, are shown in Tables 8 and 9.

Table 8

Estimation of GIP by DMBM-APGR with G0=0.3×108m3 for the constant pressure case.

No.G0 (m3)b×10MPa/104m3/dmmbe×103MPa/104m3Ge (m3)Ea (%)
1300000001.0875417202.8044938921580240024.027
21580240021.0786520402.8659228141546368751.798
31546368751.0785194492.8668398761545874091.765
41545874091.0785174642.8668537021545866641.764
51545866641.0785174342.8668539101545866531.764
No.G0 (m3)b×10MPa/104m3/dmmbe×103MPa/104m3Ge (m3)Ea (%)
1300000001.0875417202.8044938921580240024.027
21580240021.0786520402.8659228141546368751.798
31546368751.0785194492.8668398761545874091.765
41545874091.0785174642.8668537021545866641.764
51545866641.0785174342.8668539101545866531.764
Table 9

Estimation of GIP by DMBM-APGR with G0=7.5×108m3 for the constant pressure case.

No.G0 (m3)b×10MPa/104m3/dmmbe×103MPa/104m3Ge (m3)Ea (%)
17500000001.0829555792.8361687311562591622.865
21562591621.0785837582.8663950971546113971.781
31546113971.0785184272.8668470621545870221.765
41545870221.0785174492.8668537791545866601.764
51545866601.0785174342.8668539101545866531.764
No.G0 (m3)b×10MPa/104m3/dmmbe×103MPa/104m3Ge (m3)Ea (%)
17500000001.0829555792.8361687311562591622.865
21562591621.0785837582.8663950971546113971.781
31546113971.0785184272.8668470621545870221.765
41545870221.0785174492.8668537791545866601.764
51545866601.0785174342.8668539101545866531.764

It can be seen from Tables 8 and 9 that no matter how the initial iteration value is set, the final GIP estimation will stabilize to 154586653 m3. In general, the calculation process can just cease after three iterations because the difference between two adjacent iterations is small enough. The calculated results in Tables 79 prove that DMBM-APGR is also applicable to the production systems at constant bottomhole pressure.

5.1.3. Complex Production Schedules

It is not hard to discover from the previous two calculation examples that DMBM-APGR manifests a desirable application effect for both constant rate and constant pressure production systems. For further verification, the gas well production under complex schedules is also simulated by setting several combinations of constant rate production and constant pressure production as shown in Table 10 on condition that the reservoir and gas properties remain unchanged. Like the former example, the estimated value of GIP by the slope of pipwf/q vs. t relationship (see Figure 7) is 142187178 m3, which can be used as the initial value of iteration. The corresponding calculated results by DMBM-APGR are shown in Table 11.

Table 10

Complex production schedules.

t(d)TSΔt (h)Δt(d)Production scenarios
1~2511~1210.041667
q=2.0×104m3/d
1330.125
1490.375
15~16120.5
17~265241
252~592266~606241
pwf=30.9MPa
593~842607~856241
q=1.7×104m3/d
843~1097857-1111241
pwf=28.8MPa
t(d)TSΔt (h)Δt(d)Production scenarios
1~2511~1210.041667
q=2.0×104m3/d
1330.125
1490.375
15~16120.5
17~265241
252~592266~606241
pwf=30.9MPa
593~842607~856241
q=1.7×104m3/d
843~1097857-1111241
pwf=28.8MPa
Figure 7

pipwf·q1 vs. t for the complex production schedules.

Figure 7

pipwf·q1 vs. t for the complex production schedules.

Table 11

Estimation of GIP by DMBM-APGR with G0=142187178m3 for the complex production schedules.

No.G0 (m3)b×10MPa/104m3/dmmbe×103MPa/104m3Ge (m3)Ea (%)
11421871781.0811311472.8420971141559332182.651
21559332181.0817001282.8379552271561607972.801
31561607971.0817086082.8378934861561641942.803
41561641941.0817087352.8378925731561642442.803
51561642441.0817087372.8378925431561642462.803
61561642461.0817087372.8378925431561642462.803
No.G0 (m3)b×10MPa/104m3/dmmbe×103MPa/104m3Ge (m3)Ea (%)
11421871781.0811311472.8420971141559332182.651
21559332181.0817001282.8379552271561607972.801
31561607971.0817086082.8378934861561641942.803
41561641941.0817087352.8378925731561642442.803
51561642441.0817087372.8378925431561642462.803
61561642461.0817087372.8378925431561642462.803

It can be seen from Table 11 that the change rate in estimated GIP after two iterations is clearly less than 1% with an acceptable calculation error of 2.801%, when setting initial iteration value G0=142187178m3. If the iteration process continues, the estimated values will eventually stay at 156164246 m3. Tables 12 and 13 exhibit the calculation results when other initial values of iteration are given.

Table 12

Estimation of GIP by DMBM-APGR with G0=0.3×108m3 for the complex production schedules.

No.G0 (m3)b×10MPa/104m3/dmmbe×103MPa/104m3Ge (m3)Ea (%)
1300000001.0897225372.7796121571594385564.958
21594385561.0818277612.8370262541562119312.834
31562119311.0817105102.8378795931561649592.803
41561649591.0817087632.8378923761561642552.803
51561642551.0817087372.8378925431561642462.803
No.G0 (m3)b×10MPa/104m3/dmmbe×103MPa/104m3Ge (m3)Ea (%)
1300000001.0897225372.7796121571594385564.958
21594385561.0818277612.8370262541562119312.834
31562119311.0817105102.8378795931561649592.803
41561649591.0817087632.8378923761561642552.803
51561642551.0817087372.8378925431561642462.803
Table 13

Estimation of GIP by DMBM-APGR with G0=7.5×108m3 for the complex production schedules.

No.G0 (m3)b×10MPa/104m3/dmmbe×103MPa/104m3Ge (m3)Ea (%)
17500000001.0858748182.8075812861578502293.913
21578502291.0817707122.8374413721561890772.819
31561890771.0817096602.8378858211561646162.803
41561646161.0817087502.8378925091561642482.803
51561642481.0817087372.8378925431561642462.803
No.G0 (m3)b×10MPa/104m3/dmmbe×103MPa/104m3Ge (m3)Ea (%)
17500000001.0858748182.8075812861578502293.913
21578502291.0817707122.8374413721561890772.819
31561890771.0817096602.8378858211561646162.803
41561646161.0817087502.8378925091561642482.803
51561642481.0817087372.8378925431561642462.803

Similarly, the modifications to initial iteration values do not heavily affect the calculation accuracy of DMBM-APGR as can be seen from Tables 12 and 13, while the number of iterations to reach the termination condition may increase. The calculation process, in practice, can be terminated after 3 iterations, because the difference between two adjacent GIP estimations at that moment is less than 1%, and the calculation error is 2.803%. The calculation results for the complex production schedules, as shown in Tables 1113, provide convincing evidence for the feasibility and validity of DMBM-APGR in those situations where production rate and bottomhole flowing pressure fluctuate slightly.

5.2. Field Case Verification

Figure 8 shows the production history for gas well H-58 presented by Ibrahim et al. (2003) [40] and Wang and Ayala (2020) [17]. The reservoir and gas properties for the well are listed in Table 14. Similar to the previous numerical cases, the HY (1973) method [39] is employed to estimate the Z-factor and gas compressibility Cg, and viscosity correlations developed by Londono et al. (2002, 2005) [37, 38] are used for gas viscosity calculations. Figure 9 shows the relationship between Bg and p and that between μ and p. Figure 10 exhibits the relationship between pp and p, depicted by Equation (13).

Figure 8

q and pwf vs. t for gas well H-58.

Figure 8

q and pwf vs. t for gas well H-58.

Table 14

Reservoir and gas properties for the field case.

ParametersProperty values
Tpc
200 K
ppc
4.6 MPa
Swci
0.2618
ϕ
0.1354
pi
96.526602 MPa
Ti
3.664833×102K
Mg
19.268256 g/mol
Cϕ
8.702264×104MPa1
Cw
5.801510×104MPa1
μi
4.751752304×102cp
Zsc
0.997569589
Zi
1.747395656
Bgi
2.298698301×103m3/m3
Cti
3.589784792×103MPa1
ParametersProperty values
Tpc
200 K
ppc
4.6 MPa
Swci
0.2618
ϕ
0.1354
pi
96.526602 MPa
Ti
3.664833×102K
Mg
19.268256 g/mol
Cϕ
8.702264×104MPa1
Cw
5.801510×104MPa1
μi
4.751752304×102cp
Zsc
0.997569589
Zi
1.747395656
Bgi
2.298698301×103m3/m3
Cti
3.589784792×103MPa1
Figure 9

Bg and μ vs. p for the field case.

Figure 9

Bg and μ vs. p for the field case.

Figure 10

pp vs. p for the field case.

Figure 10

pp vs. p for the field case.

The estimated value of GIP by the slope of pipwf/q vs. t relationship shown in Figure 11 is 1.572152725×108m3, which can be used as the initial value of iteration. The corresponding calculated results by DMBM-APGR are shown in Table 15. The difference between two adjacent GIP estimations is denoted by Ed:
(73)Ed=GeG0G0×100%.
Figure 11

pipwf/q vs. t relationship for the field case.

Figure 11

pipwf/q vs. t relationship for the field case.

Table 15

Estimation of GIP by DMBM-APGR with G0 =1.572152725 × 108 m3 for the field case.

No.G0×108m3b×101MPa/104m3/dmmbe×103MPa/104m3Ge (×108 m3)Ed (%)
11.57215272513.338510261.5194388723.635570192131.248
23.6355701929.1478921471.8893446232.923779286-19.579
32.9237792868.3181608461.9739229262.798501704-4.285
42.7985017048.1319957771.9935139052.770999820-0.983
52.7709998208.0891882461.9980507582.764707878-0.227
No.G0×108m3b×101MPa/104m3/dmmbe×103MPa/104m3Ge (×108 m3)Ed (%)
11.57215272513.338510261.5194388723.635570192131.248
23.6355701929.1478921471.8893446232.923779286-19.579
32.9237792868.3181608461.9739229262.798501704-4.285
42.7985017048.1319957771.9935139052.770999820-0.983
52.7709998208.0891882461.9980507582.764707878-0.227

It can be seen from Table 15 that for a given initial iteration value of 1.572152725×108m3, the variation in estimated GIP after four iterations is less than 1% with an estimate of 2.771×108m3. Tables 16 and 17 exhibit the calculation results when other initial values of iteration are considered.

Table 16

Estimation of GIP by DMBM-APGR with G0=0.15×108m3 for the field case.

No.G0×108m3b×101MPa/104m3/dmmbe×103MPa/104m3Ge×108m3Ed (%)
10.1511.718223961.6518911563.3440621392129
23.3440621398.8495599271.9192606352.878205582-13.931
32.8782055828.2520711771.9808464952.788720219-3.109
42.7887202198.1168532421.9951128952.768778993-0.715
52.7687789938.0856998221.9984200962.764196919-0.165
No.G0×108m3b×101MPa/104m3/dmmbe×103MPa/104m3Ge×108m3Ed (%)
10.1511.718223961.6518911563.3440621392129
23.3440621398.8495599271.9192606352.878205582-13.931
32.8782055828.2520711771.9808464952.788720219-3.109
42.7887202198.1168532421.9951128952.768778993-0.715
52.7687789938.0856998221.9984200962.764196919-0.165
Table 17

Estimation of GIP by DMBM-APGR with G0 =15 × 108 m3 for the field case.

No.G0×108m3b×101MPa/104m3/dmmbe×103MPa/104m3Ge×108m3Ed (%)
11511.304763431.6877220633.273066575-78.180
23.2730665758.7688585341.9274492212.865977797-12.438
32.8659777978.2340259351.9827546672.786036399-2.789
42.7860363998.1126825801.9955590022.768160033-0.642
52.7681600338.0847267241.9985237222.764053592-0.148
No.G0×108m3b×101MPa/104m3/dmmbe×103MPa/104m3Ge×108m3Ed (%)
11511.304763431.6877220633.273066575-78.180
23.2730665758.7688585341.9274492212.865977797-12.438
32.8659777978.2340259351.9827546672.786036399-2.789
42.7860363998.1126825801.9955590022.768160033-0.642
52.7681600338.0847267241.9985237222.764053592-0.148

As shown in Figure 12, the estimated values will eventually stay at 2.762819832×108m3 if the above iteration process continues. The changes in initial iteration values have little effect on the calculation accuracy of DMBM-APGR for the field case. The calculation process, in fact, can be terminated after 4 iterations because Ed is less than 1% at that time. The above calculation results for the field case match well with the predicted value of 2.775050966×108m3 (or 9.8Bcf) presented by Ibrahim et al. (2003) and the estimate of 2.690100426×108m3 (or 9.5Bcf) by Wang and Ayala (2020), which further demonstrates the feasibility and applicability of DMBM-APGR during boundary-dominated flow though production rate and bottomhole pressure may fluctuate slightly.

Figure 12

Iteration processes of DMBM-APGR with different initial values for the field case.

Figure 12

Iteration processes of DMBM-APGR with different initial values for the field case.

5.3. Discussion

There are many methods for determination of GIP, the majority of which, nevertheless, fail to calculate GIP of abnormally pressured gas reservoirs with low cost and acceptable accuracy. For example, the conventional material balance method needs average reservoir pressure data by shutting the gas well or conducting a well testing destined for static pressure. Most methods for production decline analysis (or rate transient analysis), however, do not consider the elastic effects of rock pore and irreducible water. As a matter of fact, both formation and water compressibility effects are supposed to be considered in the production performance analysis.

Therefore, this study is aimed at offering a reliable approach for determining gas reserves applicable to abnormally high-pressure gas reservoirs with desirable accuracy and low cost. The novelties of the work lie in the following. (1) The iterative method DMBM-APGR presented in the paper reckons with the influence of rock and water compressibility, and avoids well shutdown or well testing since it only requires production data (q, pwf, Gp, etc.) available to determine G, which boasts low cost and high convenience. (2) The dynamic material balance equation is strictly derived from the gas flow model for a closed system with abnormally high pressure. (3) The method combines the static material balance equation with the dynamic material balance equation for overpressured gas reservoirs to design the iterative steps for GIP estimation. (4) A practical selection strategy of an initial iteration value of GIP is proposed though the method is not sensitive to choices of the initial value.

Three numerical cases and a field case validate the applicability of DMBM-APGR to overpressured gas reservoirs with constant production rate, constant bottomhole pressure, and variable rate/variable pressure conditions. The changes in initial iteration values have little effect on the calculation accuracy of DMBM-APGR for all the cases because the estimated GIP by the method always ultimately stabilizes to a fixed value independent of initial iteration values, which shows its robustness and reliability.

6. Conclusions

  • (1)

    A new definition of total compressibility, pertinent to the compressibilities of gas, rock, and irreducible water, is developed in this paper, and the dynamic material balance equation for abnormally pressured gas reservoirs is rigorously proved using the gas flow theory with a variable flow rate, too

  • (2)

    Based on the relationship between the quasipressure corresponding to bottomhole pressure and the one corresponding to average reservoir pressure (i.e., dynamic material balance equation) and the material balance principle of overpressured gas reservoirs, the authors present the dynamic material balance method for abnormally pressured gas reservoir (DMBM-APGR), an iterative method for calculating GIP of gas reservoirs without water drive by production data. The method not only considers the pore reduction in reservoir rock and the volume expansion of bound water, but also embodies the fluctuations in production rate and changes in gas properties

  • (3)

    Simulated cases under various production scenarios and a field case are used to verify the applicability and effectiveness of the method. The paper creatively proposes a selection strategy of an initial iteration value of GIP where initial value G0 can be set as the estimated value calculated by the slope of pipwf/q vs. t straight line. Generally, the accuracy requirement can be fulfilled after 2 to 4 iterations using DMBM-APGR, and the calculation error of gas in place is less than 5%, acceptable and desirable for reservoir engineering application

  • (4)

    DMBM-APGR is not sensitive to the selection of the initial iteration value of GIP, and the estimated gas in place by the method will always ultimately stabilize to a fixed value, which shows that it is robust and reliable. It is recommended to adopt the selection strategy presented here for initial values of iteration; otherwise, the number of iterations to satisfy the stopping criterion may rise

Data Availability

The gas properties data used for numerical simulation and filed application are generated by Londono et al. [37, 38] and Hall-Yarborough [39] correlations (provided in Figures 3, 4, 9, and 10 in the manuscript). The simulated production data (such as production rate, bottomhole pressure, and cumulative gas production) under various scenarios are derived from Eclipse, a commercial reservoir simulator. The field case is derived from Wang and Ayala (2020) [17] and Ibrahim et al. (2003) [40].

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work has been supported by the Department of Asia-Pacific E&P and Department of Middle East E&P, Research Institute of Petroleum Exploration and Development, PetroChina. We gratefully acknowledge the financial support from the National Science and Technology Major Project of China (Grant No. 2017ZX05030-003). Our thanks also go to Jinyu Guo, Lang Li and Ping Chen for revision suggestions.

1.
Satter
A.
Iqbal
G. M.
Reservoir Engineering: the Fundamentals, Simulation, and Management of Conventional and Unconventional Recoveries
 , 
2016
Waltham
Gulf Professional Publishing
2.
Gonzalez
F. E.
Blasingame
T. A.
A quadratic cumulative production model for the material balance of an abnormally pressured gas reservoir
SPE Western Regional and Pacific Section AAPG Joint Meeting
2008
Bakersfield, California
(pg. 
1
-
33
)
3.
Liu
N. Q.
Practical Modern Well Test Interpretation
 , 
2008
5th
Beijing
Petroleum Industry Press
4.
Wu
Z.
Cui
C.
Lv
G.
Bing
S.
Cao
G.
A multi-linear transient pressure model for multistage fractured horizontal well in tight oil reservoirs with considering threshold pressure gradient and stress sensitivity
Journal of Petroleum Science and Engineering
 , 
2019
, vol. 
172
 (pg. 
839
-
854
)
[PubMed]
5.
Mattar
L.
Anderson
D. M.
A systematic and comprehensive methodology for advanced analysis of production data
SPE Annual Technical Conference and Exhibition
2003
Denver, Colorado
(pg. 
1
-
14
)
6.
Sun
H. D.
Advanced Production Decline Analysis and Application
 , 
2015
Waltham, MA
Gulf Professional Publishing
7.
Wu
Z.
Dong
L.
Cui
C.
Cheng
X.
Wang
Z.
A numerical model for fractured horizontal well and production characteristics: comprehensive consideration of the fracturing fluid injection and flowback
Journal of Petroleum Science and Engineering
 , 
2020
, vol. 
187, article 106765
 
8.
Mattar
L.
McNeil
R.
The “flowing” material balance procedure
46th Annual Technical Meeting of The Petroleum Society of CIM
1995
Banff in Alberta
(pg. 
1
-
6
)
9.
Mattar
L.
McNeil
R.
The “flowing” gas material balance
Journal of Canadian Petroleum Technology
 , 
1998
, vol. 
37
 
2
(pg. 
52
-
55
)
10.
Blasingame
T. A.
Lee
W. J.
Variable-rate reservoir limits testing of gas wells
SPE Gas Technology Symposium
1988
Dallas, Texas
(pg. 
43
-
56
)
11.
Mattar
L.
Anderson
D.
Dynamic material balance (oil or gas-in-place without shut-ins)
Petroleum Society’s 6th Canadian International Petroleum Conference
2005
Calgary
(pg. 
1
-
9
)
12.
Mattar
L.
Anderson
D.
Dynamic material balance—oil-or gas-in-place without shut-ins
Journal of Canadian Petroleum Technology
 , 
2006
, vol. 
45
 
11
(pg. 
7
-
10
)
13.
Ismadi
D.
Kabir
C. S.
Hasan
A. R.
The use of combined static and dynamic material-balance methods in gas reservoirs
SPE Asia Pacific Oil and Gas Conference and Exhibition
2011
Jakarta, Indonesia
(pg. 
1
-
16
)
14.
Ismadi
D.
Kabir
C. S.
Hasan
A. R.
The use of combined static- and dynamic-material-balance methods with real-time surveillance data in volumetric gas reservoirs
SPE Reservoir Evaluation & Engineering
 , 
2012
, vol. 
15
 
3
(pg. 
351
-
360
)
[PubMed]
15.
Zhang
M.
Ayala
L. F.
Gas-production-data analysis of variable-pressure-drawdown/variable-rate systems: a density-based approach
SPE Reservoir Evaluation & Engineering
 , 
2014
, vol. 
17
 
4
(pg. 
520
-
529
)
[PubMed]
16.
Stumpf
T. N.
Ayala
L. F.
Rigorous and explicit determination of reserves and hyperbolic exponents in gas-well decline analysis
SPE Journal
 , 
2016
, vol. 
21
 
5
(pg. 
1843
-
1857
)
[PubMed]
17.
Wang
Y.
Ayala
L. F.
Explicit determination of reserves for variable-bottomhole-pressure conditions in gas rate-transient analysis
SPE Journal
 , 
2019
, vol. 
25
 
1
(pg. 
369
-
390
)
18.
Qin
T. L.
Li
D.
Chen
Y. Q.
Practical Reservoir Engineering Methods
 , 
1989
Beijing: Petroleum Industry Press
19.
Chen
Y. Q.
Li
D.
Modern Petroleum Reservoir Engineering
 , 
2001
Beijing: Petroleum Industry Press
20.
Al-Hussainy
R.
Ramey
H. J.
Crawford
P. B.
The flow of real gases through porous media
Journal of Petroleum Technology
 , 
1966
, vol. 
18
 
5
(pg. 
624
-
636
)
21.
Al-Hussainy
R.
Ramey
H. J.
Application of real gas flow theory to well testing and deliverability forecasting
Journal of Petroleum Technology
 , 
1966
, vol. 
18
 
5
(pg. 
637
-
642
)
22.
Russell
D. G.
Goodrich
J. H.
Perry
G. E.
Bruskotter
J. F.
Methods for predicting gas well performance
Journal of Petroleum Technology
 , 
1966
, vol. 
18
 
1
(pg. 
99
-
108
)
23.
Agarwal
R. G.
Real gas pseudo-time - a new function for pressure buildup analysis of MHF gas wells
SPE Annual Technical Conference and Exhibition
1979
Las Vegas, Nevada
(pg. 
1
-
12
)
24.
Lee
W. J.
Holditch
S. A.
Application of pseudotime to buildup test analysis of low-permeability gas wells with long-duration wellbore storage distortion
Journal of Petroleum Technology
 , 
1982
, vol. 
34
 
12
(pg. 
2877
-
2887
)
[PubMed]
25.
Fraim
M. L.
Wattenbarger
R. A.
Gas reservoir decline-curve analysis using type curves with real gas pseudopressure and normalized time
SPE Formation Evaluation
 , 
1987
, vol. 
2
 
4
(pg. 
671
-
682
)
[PubMed]
26.
Meunier
D. F.
Kabir
C. S.
Wittmann
M. J.
Gas well test analysis: use of normalized pressure and time functions
SPE Annual Technical Conference and Exhibition
1984
Houston, Texas
(pg. 
1
-
16
)
27.
Meunier
D. F.
Kabir
C. S.
Wittmann
M. J.
Gas well test analysis: use of normalized pseudovariables
SPE Formation Evaluation
 , 
1987
, vol. 
2
 
4
(pg. 
629
-
636
)
[PubMed]
28.
Ye
P.
Ayala
H. L. F.
A density-diffusivity approach for the unsteady state analysis of natural gas reservoirs
Journal of Natural Gas Science and Engineering
 , 
2012
, vol. 
7
 (pg. 
22
-
34
)
[PubMed]
29.
Vardcharragosad
P.
Ayala
L. F.
Rate-time forecasting of gas reservoirs with significant transient flow: a density-based method
Journal of Unconventional Oil and Gas Resources
 , 
2015
, vol. 
11
 (pg. 
111
-
126
)
[PubMed]
30.
Mohr
P. J.
Newell
D. B.
Taylor
B. N.
CODATA recommended values of the fundamental physical constants: 2014
Journal of Physical and Chemical Reference Data
 , 
2016
, vol. 
45
 
4, article 043102
[PubMed]
31.
Mohr
P. J.
Newell
D. B.
Taylor
B. N.
CODATA recommended values of the fundamental physical constants: 2014
Reviews of Modern Physics
 , 
2016
, vol. 
88
 
3, article 035009
[PubMed]
32.
Kong
X. Y.
Advanced Fluid Mechanics in Porous Medium
 , 
2010
Hefei
University of Science and Technology of China Press
33.
Blasingame
T. A.
Lee
W. J.
Variable-rate reservoir limits testing
Permian Basin Oil and Gas Recovery Conference
1986
Midland, Texas
(pg. 
361
-
369
)
34.
Anderson
D. M.
Mattar
L.
An improved pseudo-time for gas reservoirs with significant transient flow
Canadian International Petroleum Conference
2005
Calgary, Alberta
(pg. 
1
-
11
)
35.
Anderson
D. M.
Mattar
L.
An improved pseudo-time for gas reservoirs with significant transient flow
Journal of Canadian Petroleum Technology
 , 
2007
, vol. 
46
 
7
(pg. 
49
-
54
)
36.
Ye
P.
Ayala
H. L. F.
Straightline analysis of flow rate vs. cumulative-production data for the explicit determination of gas reserves
Journal of Canadian Petroleum Technology
 , 
2013
, vol. 
52
 
4
(pg. 
296
-
305
)
[PubMed]
37.
Londono
F. E.
Archer
R. A.
Blasingame
T. A.
Simplified correlations for hydrocarbon gas viscosity and gas density - validation and correlation of behavior using a large-scale database
SPE Gas Technology Symposium
2002
Calgary
(pg. 
1
-
16
)
38.
Londono
F. E.
Archer
R. A.
Blasingame
T. A.
Correlations for hydrocarbon gas viscosity and gas density - validation and correlation of behavior using a large-scale database
SPE Reservoir Evaluation & Engineering
 , 
2005
, vol. 
8
 
6
(pg. 
561
-
572
)
39.
Hall
K. R.
Yarborough
L.
A new equation of state for Z-factor calculations
Oil and Gas Journal
 , 
1973
, vol. 
71
 
7
(pg. 
82
-
85, 90, 92
)
40.
Ibrahim
M.
Wattenbarger
R. A.
Helmy
W.
Determination of OGIP for wells in pseudosteady-state-old techniques, new approaches
SPE Annual Technical Conference and Exhibition
2003
Denver, Colorado
(pg. 
1
-
12
)
This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.