A novel material-balance approach for estimating in-place volumes of gas and water in gas reservoirs with aquifer support

Applications of the traditional static material-balance method in gas reservoirs become a challenge with production maturity due to variability in aquifer influx, infill drilling, and production-operational changes in offset wells, among others. Besides, some existing modeling approaches involve a trial-and-error method, making the solution outcomes nontrivial. This study proposes a new methodology for analyzing production data involving water-drive gas reservoirs. The main findings of this study include the following: (1) A straight-line plot that yields gas and water in-place volumes, (2) A modified-(pav/z)* plot exhibits a straight-line with an x-intercept of gas initially-in-place, similar to that in a conventional-(pav/z) plot, (3) A new definition of degree of aquifer support that is quantifiable using production data. Synthetic data verified the proposed modeling approach, whereas a field dataset provided validation.


p *
Average reservoir pressure in a gas reservoir without aquifer, psia S gr Residual gas saturation, % S wi Irreducible water saturation, % V ba Bulk volume of an aquifer, RB V bap Predicted bulk volume of an aquifer, RB V bg Bulk volume of a gas reservoir, RB V bgp Predicted bulk volume of a gas reservoir, RB V cw Volume of connate water, RB V cwf Volume of connate water in a water-influx zone, RB V cwg Volume of connate water in a gas reservoir, RB V f Bulk volume of a water-influx zone, RB V g Volume of gas in a gas reservoir, RB V rg Volume of residual gas in a water-influx zone, RB V wa Volume of water in an aquifer, RB V wf Volume of aquifer water in a water-influx zone, RB 1 3

Introduction
Reservoir engineers have relied on the material-balance equation or MBE to estimate gas initially-in-place (G). Some of the landmark studies, particularly those of Havlena and Odeh (1963), Cole (1969), and Campbell and Campbell (1978), developed tools for a water-drive gas reservoir. The presence of aquifer support increases the number of unknowns. This challenge of solving the gas MBEs appears in Bruns et al. (1965), Chierici et al. (1967), Dake (1983), Tehrani (1985), and Vega and Wattenbarger (2000), among others. Bruns et al. (1965) and Agarwal et al. (1965) have documented the effects of water influx on the (p av /z) plot.
Specifically, Bruns et al. (1965) used various aquifer models, those of Schilthuis (1936), Hurst simplified (1943), and van Everdingen and Hurst (1949). Decades later, Elahmady et al. (2007) showed a non-uniqueness problem, leading to an overestimation of in-place gas volume, G. This issue arises due to the time-dependent nature of the aquifer influx model. To that end, Izgec and Kabir (2010) showed that a simplified form of the van Everdingen and Hurst aquifer-influx model could discern the variable nature of aquifer influx. Yildiz and Khosravi (2007) showed the performance of their analytical tools in edge-water and bottom-waterdrive systems. In a follow-up study, Yildiz (2008) offered a hybrid approach involving a modified version of the Roach (1981) plot and the McEwen (1962) method to ascertain in-place gas volume with production data in water-drive reservoirs. The illustrative examples supported the goodness of this methodology, particularly in strong waterdrive systems. Some of the studies in the modern era include those of Moghadam et al. (2011), exploring the application of the p/z method in coalbed-methane, overpressured, and water-drive reservoir systems. Wang and Teasdale (1987) established the material balance equation (MBE) for a fractured vuggy gas reservoir with a bottomwater drive, including stress sensitivity and gravity segregation. The proposed MBE can be applied to estimate the original gas-in-place and the distribution of the reserve in the matrix, fracture, and cavity.
More recently, Yu et al. (2019) showed that the Blasingame method and the flowing-material balance (FMB) approach could lead to visual identification and analysis of the early aquifer-influx period data. Yu et al. (2019) proposed no-aquifer influx period-type curves and early aquifer-influx period-type curves, which can be applied to determine reservoir and aquifer parameters. Well history production data are divided into three influx periods: noaquifer, early aquifer, and middle-late aquifer. It yields the dimensionless drainage radius, initial gas-in-place, other reservoir parameters, and aquifer energy. Also, Kazemi and Ghaedi (2020) presented a semi-analytical method for the FMB for edge-water-drive systems with reservoir pressure exceeding 3000 psi. For lower pressures, it requires a correction factor. Zaremoayedi et al. (2022) proposed a new production data analysis approach by applying a modified compressibility parameter to the conventional FMB equations to estimate the initial gas in place accurately and average reservoir pressure in a non-volumetric naturally fractured gas condensate reservoir. Without using complex multi-phase reservoir modeling, they modified the isothermal rock compressibility factor to model the aquifer's influence on reservoir performance.
Despite the plethora of studies in literature, aquifer volume (W e ) estimation becomes a challenge due to the uncertainty affecting the microscopic displacement efficiency and the conformance factor, as Chierici et al. (1967) pointed out. None of the previous works show a straight-line analysis. This study aims to develop a new methodology for a waterdrive gas reservoir to estimate the aquifer volume during the boundary-dominated flow period. Our ultimate objective is to create a holistic modeling approach, given the challenges in the modern era involving complex fluid PVT, overpressured reservoirs, and aquifer influx. This investigation develops a modified-(p av /z)* plot handling the variable water-drive situation and quantifies various degrees of aquifer support in a conventional reservoir. Estimations of gas and aquifer in-place volumes became the underlying motivation, given that this item has remained unexplored. Of course, the knowledge of aquifer strength associated with each well helps develop a reservoir's depletion strategy.

Material balance equation in a water-drive gas reservoir
Let us write the material balance equation of a water-drive gas reservoir as where the fluid withdrawal (F), the cumulative expansion of gas (E g ), the incremental expansion of the connate water and reduction in the reservoir-pore volume (E fw ), and the cumulative expansion of water (E a ) are defined, respectively, as the following: W is water initially-in-place. While the classical gas MBE uses W e , Eq. 1 uses W. Dividing Eq. 1 by (E g + E fw ) and rearranging yields the following expression: A dimensionless plot of F/[E g + E fw ] vs. E a /[E g + E fw ], as illustrated in Fig. 1, yields a straight line with a slope of W and G from the intercept. The larger the aquifer, the steeper the slope of a straight line. However, the y-intercept is independent of an aquifer's size. While both a straight-line plot and Cole plot yield the same y-intercept of G, this study offers straightline profiles. Given that reservoir connectivity is a critical parameter not explicitly considered, Fig. 1 cannot assess the correct aquifer size. So, this index implies the relative contribution of water influx in the context of the drive mechanism. Once G and W are known, the bulk volumes of a gas reservoir (V bg ), an aquifer (V ba ), and a water-influx zone (V f ) can be calculated, respectively, as

Correction for the (p av /z) plot
A p av /z plot of a water-drive-gas reservoir deviates from a straight line because of the presence of an aquifer. A part of a gas reservoir with a volume of V wf represents the aquifer water. Without an aquifer, this volume gets filled with gas. Then, gas-pore volumes for models with and without aquifer are (V g + V rg ) and (V g + V rg + V wf ), respectively. V g is the volume of gas in a gas reservoir, and V rg is the residual gas volume in a water-influx zone. Therefore, the average reservoir pressure of a model without an aquifer (p * ) should be lower than that for an aquifer (p). The value of p is corrected to p * using the real-gas law. Then, the (p av /z) * relation becomes The values of V g , V rg , and V wf can be estimated using the following equations: Similar to the conventional (p av /z) plot, a plot of (p av /z) * vs. G p yields a straight line with y-intercept = (p/z) i and x-intercept of G. Note that this plot requires prior knowledge of gas in place which can be estimated from Fig. 1. The real objective of this section is to verify whether the value of G is correctly estimated.

Drive indices
Following Pletcher's pot aquifer (2002) definition of drive indices, let us rewrite Eq. 1 as where gas-drive index (I GD ), compaction-drive index (I CD ), and water-drive index (I WD ) are defined, respectively, as While Pletcher (2002) used W e in I WD , this study uses (WE a ). Quality of production data, reservoir complexity, infill drilling, and operational activities may pose challenges for Eq. 14. So, Kabir et al. (2016) recommended that I WD gets evaluated using the following expression:

Degree of aquifer support
Two dimensionless variables assisted in quantifying the degree of aquifer support, D as . Let us define D as within a range of 0 and 1; zero implies no energy support, whereas (14)

Fig. 3 Flowchart of the proposed methodology
• For a gas reservoir with weak aquifer support, it has an x-intercept @ (p av /p i ) = 0; 0 < D as < 0.5. When the reservoir gets fully depleted, no water breakthrough occurs. • For a gas reservoir with strong aquifer support, it has a y-intercept @ (V f /V bg ) = 100%; 0.5 < D as < 1. When the gas recovery factor (RF) is maximum, there is still energy available in the system. • For a gas reservoir with infinite aquifer support, it is a horizontal line @ (p av /p i ) = 100%; D as = 1. In other words, no pressure depletion occurs in this case.

Aquifer expansion & production (RB/D) t (D)
Aquifer expansion Aquifer production One can estimate the value of D as as where Figure 3 shows the flowchart of the proposed methodology, which provides a visual representation of the step-bystep process and will help to make the calculation procedure more transparent and accessible.

Model validation and verification
Two synthetic examples validate the modeling approach, and one field example from the literature verifies the proposed method.
Example 1: Synthetic data The data in this example are from Run #9 of Kabir et al. (2016). This layered Eclipse model has about 198,000 active cells with a cell size of 100 ft × 100 ft × 3 ft. For simplicity and transparency, the model has uniform rock properties. The authors used numerical modeling to generate production data using the following rock and fluid properties in Table 1.
Some of the other relevant attributes of this modeling approach appear in the context of presenting the following example. Production profiles of gas and water (q g and q w ) and flowing bottom-hole pressure (p wf ) appear in Fig. 4. The gas rate remained constant at 50 MMscf/D for about five years without water production. When the well started producing water, the gas rate began to decline. However, the p wf remained stable until about 2000 days before its decline started and reached a minimum value of 1000 psia. Figure 5a shows the classical material-balance straightline plot. As expected, the late-time data in the boundarydominated flow or BDF period appear as a straight line, whereas the early-time transient data do not. The straight line yields the following parameters: G = 115 Bscf; from y-interceptV bg = 275 MMRB W = 417 MMSTB; from slopeV ba = 2090 MMRB The aquifer is about 7.6 times larger than the gas reservoir relative to the model input value of 10. During transient flow, the values on the left-hand side (LHS) of Eq. 1 are significantly less than those on the right-hand side (RHS). The ratios of (LHS) and (RHS) appear in Fig. 5b. The ratio is less than one during the transient flow period and reaches approximately one during the BDF period. In Eq. 1, using G = 115 Bscf, one can estimate the water initially-in-place, W a . The aquifer volume actively contributes to production    1 3 beyond the transient period. W a is much less than W during transient flow but increases with time. As expected, W a equals W during the BDF period when the entire aquifer actively contributes to the production. Figure 6a shows two bulk volume profiles of the waterinflux zone (V f ) in a dimensionless form. The results from these two equations are slightly different. Without any water production, V f increases monotonically with producing time.
Water production begins after about 1800 days, leading to decreasing slope of the profiles. The profile slope is positive if an instantaneous aquifer production rate (q w B w ) is less than the aquifer expansion rate. As Fig. 6b illustrates, after 2830 days of production, q w B w becomes more significant than an instantaneous aquifer expansion, causing a declining V f trend.
Profiles of V g , V rg , and V wf appear in Fig. 7a. The slope of V f can explain the individual profiles. While changes in V rg and V wf are in the same direction as V f , the V g trend occurs in the opposite direction. Then, one can calculate a profile of [(V g + V rg )/(V g + V rg + V wf )] for estimating (p av /z)* in Eq. 10. In this case, it declined from 1.0 to around 0.2. Figure 7b shows the profiles (p av /z) and (p av /z)* vs. G p . A (p av /z) plot is concave downward and yields an overestimated G, especially during early production. At RF = 50%, the (p av /z) plot produces a straight line with R 2 of 0.992, but a G of 375 Bscf, reflecting a 226% overestimation. In other words, the early portion of the (p av /z) plot becomes unsuitable for G estimation. However, with correction, the (p av /z)* resulted in a G of 116 Bscf, close to the input value of 115 Bscf. Figure 8a shows different drive indices. I WD dominates production at early times, while I GD dominates during the late production period. As expected, I CD is insignificant. Figure 8b displays the production contributions from individual energy sources, and Fig. 8c shows the relative   expansion factors. Note that one for the gas drive increases significantly faster than others when depleting the reservoir. The D as profile in Fig. 8d shows that it increases before reaching the maximum during the early production period. When the well starts producing water, D as declines to 0.5. Note that D as is estimated using historical production data only during the BDF period. Figure 9a shows the results of Cole and modified-Cole plots. The profiles increase very early times and reach a peak before declining. The minimum G p B g /(B g − B gi ) value is 132 Bscf, and G p is 111 Bscf. Therefore, the range of G is 111 -132 Bscf. The result of the Roach (1981) plot appears in Fig. 9b. The value of G is 109 Bscf, derived from a reciprocal of a slope of the straight line, which is reasonably close to the input value of 115 Bscf.
Example 2: Synthetic data The data in this example is from Run #15 of Kabir et al. (2016), with an initial pressure of 13,248 psia. The gas rate is kept constant at 30 MMscf/D with declining p wf , as illustrated in Fig. 10a. When p wf reaches 1000 psia, the gas rate starts to decrease. No water production occurred throughout history. The aquifer's size is about one-half of the gas reservoir, as shown in Fig. 10b  involving 200 md permeability with 20% porosity and water saturation of 20% above the gas/water contact. Figure 11a shows that two straight lines emerged, given a decreased aquifer volume, and the well's proximity to the aquifer delayed the total system response. The first straight line leads to the following results: G = 87.2 Bscf; W = 377 MMSTB, whereas the second straight line yields G = 101 Bscf; V bg = 316 MMRB; W = 408 MMSTB; V ba = 2,040 MMRB The second straight line corresponds to a more extensive gas-water system. The aquifer size is about 6.5 times larger than the gas volume. The ratios of the LHS to the RHS of Eq. 1 suggest using values of G and W from the second straight line. The well produces from a smaller reservoir system during BDF-1 than BDF-2.
In the absence of water production, V f increases monotonically with production time, as Fig. 12a shows. When p wf reaches the minimum, the gas rate starts to decline. V f reaches the maximum during this period as the aquifer expansion declines, as Fig. 12b illustrates. Figure 6b shows an increasing trend of aquifer expansion during the early flow period, while Fig. 12b shows a decreasing trend. Let us attribute this difference to the two different production periods. A higher production rate causes an increasing pressure depletion rate, leading to a growing aquifer expansion rate. This reality explains the two BDF periods in Fig. 11. Figure 13a presents the profiles of V g , V rg , and V wf . They are all monotonic functions in the absence of water production. The (p av /z) correction factor decreases from 1.0 down to 0.3. On the (p av /z) plot in Fig. 13b, two trend lines emerged, corresponding to two straight lines in Fig. 11a. With the proper correction, (p av /z)* yields a straight line, with G = 102 Bscf. Figure 14a shows different drive indices. Their behaviors are similar to the ones in Example 1. Figure 14b displays an increasing D as trend at early times, then reaches the maximum, followed by a decline before stabilizing at a value of 0.38. Figure 15a shows the results of Cole and modified-Cole plots. The profiles on plots decrease with time, indicating weak aquifer support. The minimum G p B g /(B g -B gi ) value is 114 Bscf, and G p is 97 Bscf; therefore, the G range spans 97 to 114 Bscf. The result of the Roach (1981) plot appears in Fig. 15b. The value of G is 102 Bscf, derived from a reciprocal slope of the straight line.
Example 3: Field data This data from Reservoir A appeared in the study of Rossen (1975). The pertinent reservoir properties appear in Table 2. Figure 16 exhibits a cumulative gas production of 200 Bscf and a pressure depletion of 3938 psig. No water production appeared in the study. Figure 17 presents the results of our analysis. The estimated values of G are 417 Bscf from a straight-line plot in Fig. 17a, and 425 Bscf from the proposed (p av /z)* graph in Fig. 17b and W is 1640 MMSTB. The aquifer volume is about five times the gas reservoir. Figure 18a shows that I GD dominates production for the entire period, given that the aquifer size is insignificant. D as becomes stable after 700 days of production and is slightly less than 0.2, as shown in Fig. 18b. Figure 19a presents the Cole and modified-Cole plots. Both profiles indicate moderate water drive. The Roach plot in Fig. 19b yields a G value of 431 Bscf. Given the low-D as value of about 0.2, meaning a weak aquifer influx, this G value is close to the other methods shown in Fig. 17.

Discussion
The proposed method for a water-drive gas reservoir overcomes the obstacles that Cole and modified-Cole methods experience. While the Roach method is valid for weak and medium aquifer support, it yields erroneous results when the aquifer support gets stronger. In contrast, the proposed method works well for all degrees of aquifer support, given that it captures the critical mechanisms that are in play in the model. However, the proposed methodology has the following requirements and limitations.
• Besides flow rate, periodic average-reservoir pressure, rock, and fluid properties are essential items. • The BDF flow period is a requirement. • Tight reservoirs may have long transient-flow periods and require long shut-in periods to represent the average system pressure. • The approach may not work for a system with partial communication between the gas reservoir and the aquifer if they operate at different average pressures.  1 3 The proposed method shows that the early production period does not lend itself to a meaningful analysis due to the absence of the boundary-dominated flow (BDF) period. In other words, the whole system does not actively contribute to various production mechanisms during the transient flow period. In this context, the larger the aquifer than the inplace gas volume (V ba /V bg ), the longer it takes for the entire production system to reach the BDF due to the increased water-drive index I WD . The degree of aquifer support, D as , as introduced here, measures an aquifer's efficacy. The two field examples showed a range of D as from about 0.2 to 0.4, whereas the two synthetic cases had 0.4 to 0.6, capturing the evolving nature of its trend. In reservoir management, one needs to learn both the aquifer strength (D as ) associated with each well and the interwell connectivity through primary capacitance-resistance modeling or PCRM, as shown by Izgec and Kabir (2012). Then, one can adjust withdrawal rates in each reservoir compartment to ensure that premature water breakthrough does not occur. This strategy can assist in maximizing gas recovery.
The increasing aquifer size or the V ba /V bg ratio also adversely impacts the (p av /z) trend on the straight-line plot; the profile with a negative slope signifies the aquifer support. This trend implies that evaluating G from the conventional (p av /z) method will yield an overestimation. In theory, the values of G from the traditional (p av /z) straightline plot and the modified-(p av /z)* plot, as introduced here, should be the same in a volumetric system. Application of the proposed method to analyze production data may yield a difference in these two values. In other words, their difference reflects the uncertainties in G and W. These differences may indicate data quality relating to production data, reservoir pressure, and rock and fluid properties. However, the degree of aquifer influx may be the principal underlying reason. In this context, Brigham's (1997) simplified aquifer model suggests that the pressure change is independent of mobility (k/μ) at each timestep during the BDF period. This reality reassures the estimation of the aquifer volume as presented here.
Studies abound with various forms of rate-transient analysis for estimating in-place gas volume. These methods involve conventional rate-transient analysis (Palacio and Blasingame 1993), dynamic material balance (Mattar et al. 2006), transient-PI (Medeiros et al. 2010), and combined static and dynamic material balance (Ismadi et al. 2012). Those and related studies make serious value propositions in estimating in-place gas volume. Nevertheless, static, and dynamic material-balance methods must objectively arrive at realistic in-place solutions in a waterdrive system.
In that regard, our underlying intention focused on exploring the value proposition of estimating both the aquifer and gas in-place volumes to assist the development of numerical models to study the operational strategy in water-drive reservoirs to maximize recovery. That is primarily because all producers experience different degrees of aquifer support, given the variability in the formation permeability and potential geologic baffles leading to the aquifers.

Conclusions
The following conclusions appear pertinent here: • This study proposes a new methodology to interpret production performance in water-drive-gas reservoirs involving various degrees of aquifer support. Two synthetic datasets verified this method, and one field example provided validation. • The proposed definition of the degree of aquifer support (D as ) quantifies the aquifer strength; it ranges from zero (depletion drive) to one (infinite aquifer support). In this context, the D as values of the three examples spanned from about 0.2 to 0.6. Also, the novel material-balance plot yields a straight line with a slope of W (connected aquifer volume) and a y-Intercept of G (in-place gas). Production data during the early period might not be on a straight line and may have a negative slope, identifying the transient flow period. • The proposed new plot of (p av /z)* vs. G p yields a straight line in a water-drive gas reservoir, similar to a volumetric system's conventional (p av /z) plot. The proposed method estimates the in-place gas volume more accurately than Roach and Cole.
Funding This is no funding for this work.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, 1 3 provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.