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.
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 , material balance method , transient well test analysis [3, 4], and production decline analysis [5–7], 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 -factor) with pressure or time are not considered.
Blasingame and Lee (1988)  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 . 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 , 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 [13–17] 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, -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)  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)  proposed the concept of pseudotime, representing a new stage for solution of gas flow. Lee and Holditch (1982)  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) , 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
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
From the definitions of quasivariables, it can be easily noted that quasitime is a function of time , while the gas viscosity and total compressibility are the functions of pressure dependent on location and time; so is a function of location and time . Similarly, quasipressure is only related to pressure; therefore, it is also a function of and . Strictly speaking, and correspond to the same pressure , 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).
2.3. Variable Rate Solution
It is worth noting that once the vs. graph is drawn and the slope 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 vs. curve cannot be obtained in advance, the authors innovatively use the slope of vs. curve to replace the initial value of , and GIP estimated by this slope is taken as the initial value of iteration .
2.4. Dynamic Material Balance Equation for Abnormally Pressured Gas Reservoirs
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 , but values are often difficult to measure. Therefore, the iterative method, where an initial GIP value is postulated to calculate the and data and then GIP is updated by these data time and again, can be employed to determine the gas reserves. The average reservoir pressure 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
The material balance equation for abnormally pressured gas reservoirs, i.e., Equation (70), expresses a relationship between the -dependent pressure function in the overpressured gas reservoir and the cumulative gas production , from which it can be observed that a straight-line section should appear on the vs. relationship curve with its intercept and the slope or , helpful to determine the value of . If the elastic effect of rock and bound water is insignificant, that is , then can convert into , which is applicable to volumetric gas reservoirs. For overpressured gas reservoirs, nevertheless, the function is supposed to be used; that is to say, it seems unreasonable to rashly apply relation to 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 , average reservoir pressure , and material balance quasitime by the iterative method shown in Figure 2. Its detailed steps are as follows:
Collect and sort out the known data available, such as gas production rate , cumulative gas production , and bottomhole pressure
Determine the intercept of vs. curve according to Equation (34)
Calculate by Equation (71), draw vs. relationship curve, and set the intercept of the straight line as , then use the slope to determine the gas in place
Compare calculated in Step 7 with the initial value in Step 2. If their difference meets accuracy requirement, then the estimated gas in place is ; otherwise, assign to and repeat Steps 2 to 8
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 (, , , etc.) available to determine gas in place , 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 ; the original reservoir pressure is 55 MPa, the temperature is 108°C, the initial rock porosity is 0.2, and the initial irreducible water saturation 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 and actual gas in place obtained by initialization are 565864 m3 and , respectively. Other parameters about gas reservoir and natural gas are shown in Table 1.
In Table 1, denotes the radial length of the grid, denotes the angular size of each grid block, and denotes the grid height; , , and are the radial permeability, azimuthal permeability, and vertical permeability, respectively; represents the pseudocritical temperature of gas mixture and 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  is used to determine the -factor and gas compressibility . The resulting relation between gas formation volume factor and pressure and that between gas viscosity and are shown in Figure 3. In line with Equation (13), Figure 4 reveals the relationship between quasipressure function and pressure .
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 . The production schedules are shown in Table 2 where denotes time step and represents time interval.
The time-dependent data such as production rate , bottomhole flowing pressure, and cumulative gas production , 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 in Equation (34) can be also determined. Now, we only need to iteratively calculate the values of the abscissa plotting function .
To determine the material-balance quasitime requires a given GIP value . This initial value of iteration can be determined by the slope of vs. curve, that is, this slope can be used to replace in Equation (38). As shown in Figure 5, when the quasipressure and quasitime are not used, the slope of vs. curve is , 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.
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 vs. 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 , the estimation by DMBM-APGR will tend to be the same result (see Tables 4 and 5).
5.1.2. Constant Bottomhole Pressure Simulation
Setting the production schedule of gas well as without changes in other parameters generates the constant pressure production case as shown in Table 6. According to the simulation results, the slope of vs. curve is ·(104 m3)-1 as shown in Figure 6, so the GIP estimation from this slope is 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 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.
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 7–9 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 vs. 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.
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 . 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.
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 11–13, 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)  and Wang and Ayala (2020) . The reservoir and gas properties for the well are listed in Table 14. Similar to the previous numerical cases, the HY (1973) method  is employed to estimate the -factor and gas compressibility , and viscosity correlations developed by Londono et al. (2002, 2005) [37, 38] are used for gas viscosity calculations. Figure 9 shows the relationship between and and that between and . Figure 10 exhibits the relationship between and , depicted by Equation (13).
It can be seen from Table 15 that for a given initial iteration value of , the variation in estimated GIP after four iterations is less than 1% with an estimate of . Tables 16 and 17 exhibit the calculation results when other initial values of iteration are considered.
As shown in Figure 12, the estimated values will eventually stay at 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 is less than 1% at that time. The above calculation results for the field case match well with the predicted value of (or ) presented by Ibrahim et al. (2003) and the estimate of (or ) 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.
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 (, , , etc.) available to determine , 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.
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
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
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 can be set as the estimated value calculated by the slope of vs. 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
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
The gas properties data used for numerical simulation and filed application are generated by Londono et al. [37, 38] and Hall-Yarborough  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)  and Ibrahim et al. (2003) .
Conflicts of Interest
The authors declare that they have no conflicts of interest.
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.