A simple approach to dynamic material balance for a dry-gas reservoir

The dynamic material balance methodology can be used to estimate gas initially-in-place using only production and PVT data. With this methodology, reservoir pressure is obtained without requiring the well to be shut in; it is therefore superior to the static material balance method since there is no loss of production. However, the technique requires iterative calculations and numerical integration of gas pseudotime and is quite complex to implement in practice. A simpler and equally accurate methodology is proposed in this study. It requires only production and PVT data and also does not rely on a shut-in pressure survey. In addition, it requires neither iterative calculations nor numerical integration of gas pseudotime. The results of the analysis include gas initially-in-place and gas productivity index, and can easily be extended to production forecasting. Gas initially-in-place is evaluated with a high degree of reliability. The methodology is successfully tested with two simulated cases and one field case, giving high-accuracy results.


Introduction
One of the most reliable methods for hydrocarbon-in-place estimation is material balance (MB). It requires production data and average reservoir pressures at different points in time, which are obtained by running pressure surveys with the well shut-in. The resulting production loss is the main drawback of this method and is especially critical for tight reservoirs in which lengthy shut-in pressure surveys would be needed to obtain accurate reservoir pressure estimates.
Moving beyond the traditional MB approach, Mattar and McNeil (1998) proposed the "flowing material balance" methodology. This method requires only production data (flow rate and wellbore flowing pressure) and PVT data. MB calculations can be performed without average reservoir pressure, which means that the well does not need to be shut-in. However, this method is valid only for a well with constant flow rate and breaks down when the flow rate varies with time.
Subsequently, Mattar and Anderson (2005) presented the "dynamic material balance (DMB)" methodology which is an extension of flowing material balance. The method can be applied for both gas and oil reservoirs, with either constant flow rate or variable flow rate. It requires profiles of flow rates and wellbore flowing pressures, and PVT data, and allows flowing pressure at any point in time to be converted to the average reservoir pressure. Once the average reservoir pressure is available, without shutting-in the well, the classical material balance calculations are applicable. The method has been successfully applied to many field studies (Ojo et al. 2004;Osisanya et al. 2006;Zhao and Tian 2006;Mazloom et al. 2007;Sureshjani et al. 2014;Hussain et al. 2016). However, the drawbacks of the method are the following: • It needs iterative calculations: using an initial guess for the value of gas initially-in-place (G), iterations are performed until the values of G converge. • It requires calculation of pseudotime for gas (t ca ), initially introduced by Agarwal (1979). For varying gas flow rates, t ca is approximated using numerical integration.
Note that the pseudotime for gas, Eq. (C.11) in Agarwal's paper, does not use q g in the integration since q g is assumed to be constant. t ca is similar to the pseudo-equivalent time (t a ) which was developed by Palacio and Blasingame (1993) by extending the earlier work of Fraim and Wattenbarger (1987).
The above techniques require the calculation of pseudopressure and pseudotime. Ye and Ayala (2013) and Ye (2013a, 2013b) proposed a density-based approach. They used Carter's λ and its time-average value (β) to account for nonlinearities. For a volumetric reservoir, a gas rate can be expressed in terms of the well-known liquid exponential decline, using λ and β as rescaling parameters. They showed that the gas rate declines hyperbolically within the hyperbolic window. This window occurs during the early boundary-dominated flow (BDF) period for which there is insignificant pressure depletion (p ≅ p i ) and the value of b remains constant. Stumpf and Ayala (2016) showed that the application of type-curve and straight-line analysis techniques can be applied to determine G and reservoir properties. However, when pressure depletion becomes significant, the gas rate profile falls outside the hyperbolic window, and b varies with cumulative production. Then, the data outside the hyperbolic window cannot be used for interpretation using this methodology.
The objective of this study is to propose a simpler DMB method for a dry-gas reservoir that needs neither iterative procedures nor numerical integration of the pseudotime function. In addition, it does not require the constant-b assumption and a wider range of data can be used for an interpretation. These methodological simplifications are not traded off against accuracy.

Gas initially-in-place (G)
After a well is opened for production, three key pressures characterize the system: • Initial reservoir pressure (p i ) • Average reservoir pressure (p) • Wellbore-flowing pressure (p wf ) The pressure loss due to depletion can be estimated using the gas material balance equation which relates p i and p with recovery factor (G p /G) (Dake 1978).
The pressure loss due to inflow can be calculated using a stabilized gas flow rate (q g ) which is a function of gas productivity index (J g ) and pseudopressure drawdown [m(p)−m(p wf )] (Al-Hussainy and Ramey 1966).
where pseudopressure [m(p)] is defined as The key assumptions for these two equations are: • A volumetric dry gas reservoir • Negligible effects of water and formation compressibility • Boundary-dominated flow (BDF) The absence of non-Darcy flow effects.
Normally, the values of p i , G p , q g , and p wf are available. There are three unknowns; p, G, and J g . We cannot use this system of equations to directly solve for the unknowns since there are fewer equations than there are unknowns. This is the reason why the available methods (Palacio et al. 1993;Agarwal et al. 1999) in the literature require iterative calculations. For a given value of G, we can calculate J g at each production data point and standard deviation of J g , SD(J g ), for the whole production profile. For different values of G, we can construct a plot of SD(J g ) versus G. Theoretically, J g is constant if and only if the value of G is correct. The correct value of G is the point where SD(J g ) is the lowest.
Step-by-step procedures are as follows: 1. Assume a value of G. 2. Use Eq. (3) to calculate J g at each production data point. 3. Calculate SD(J g ) for the whole production profile. 4. Repeat the calculations in steps 1-3 for different values of G. 5. By plotting SD(J g ) versus G, find the value of G where SD(J g ) is minimum.
Alternatively, based on Eq.
(3), a plot of q g versus [m(p)−m(p wf )] yields a straight line with a y-intercept of zero if a value of G is correctly assumed. The slope of this straight line is J g . This approach requires varying flow rates within the production history for production with a constant flow rate, and the plot yields a straight line with zero slope. If the value of G is correctly assumed, [m(p)−m(p wf )] should have a single value. For noisy production data, the value of SD[m(p)−m(p wf )] is a minimum when a value of G is correctly assumed.

Decline rate (D)
Decline rate is defined as the fractional change in flow rate per unit change in time.
Combining Eqs. (2-5) leads to an expression for decline rate as follows: where the initial decline rate (D i ) and the viscosity-compressibility ratio (λ) are defined, respectively, as: The details of the derivation are set out in AppendixA. D is composed of two parts which are attributable to reservoir pressure depletion and to changes in p wf . The part of the decline rate due to reservoir pressure depletion is a function of reservoir rock and fluid properties. Note that D i is the initial decline rate for a given constant p wf .

Decline exponent (b)
The decline exponent is defined as the change in the reciprocal of decline rate per unit change in time.
Because, for a constant p wf , the decline rate is caused only by reservoir pressure depletion, Eqs. (6) and (9) can be combined to yield The details of the derivation are set out in AppendixB. The value of the decline exponent depends only on fluid properties and on the specified well operating conditions. Therefore, it can be predicted without any field production or test data.

Identification of the BDF period
The proposed methodology in this study is applicable only for data during the BDF period. Therefore, identification of the beginning of the BDF period is a critical issue. During the transient flow period, prior to BDF, the productivity index is not constant but varies with time. Araya and Ozkan (2002) defined the transient productivity index for a gas reservoir as where the material-balance pseudotime (t a ) is defined as (Palacio et al. 1993;Agarwal et al. 1999;Raghavan 1993) Agarwal et al. (1999) presented the productivity indices as a function of t a for different reservoir pore volumes. During the transient flow period, there is insignificant pressure depletion (p ≈ p i ), and p wf is decreasing with production time. Therefore, J g is decreasing during the transient flow period. During the BDF period, J g is constant for the correct reservoir pore volume. Too high and too low estimations of reservoir pore volume cause, respectively, downward and upward bends in the profile of J g versus the appropriate time function. These effects of the estimated values of G on the profiles of J g are shown in Fig. 1. Araya and Ozkan (2002) presented similar results. A reservoir of limited volume has a constant J g during the BDF period. We can, therefore, identify the beginning of the BDF period as the time when J g stops decreasing and starts to approach a constant value.

Production forecast
Application of the proposed methodology to analyze the production data yields G and J g . Let us assume that we have production data (q g,0 and G p,0 ) up to time t 0, and the future profile of p wf . We can forecast the production time t 1 when the gas flow rate declines to the value q g,1 using the following steps.
By moving to the next time step and repeating the above calculations, we can generate a production forecast profile.

Applications
Example 1: A gas well producing at two different wellbore-flowing pressures The production data are generated using Saphir software. The effects of wellbore-storage and the transient flow periods are incorporated into the production data. Table 1 shows the reservoir and fluid properties and operating conditions. The gas production profile is shown in Fig. 2. The wellbore flowing pressures (p wf ) are stepwise constant, and their values during the first and second years are 3,500 and 2,000 psia, respectively. Figure 3 shows the profiles of J g for three different values of G during the first year of production. During the transient flow period (t < 10 days), the values of J g are independent of the values of G. The J g profiles based on these three different values of G are exactly the same during transient flow, and the values of J g are on a decreasing trend. These profiles are significantly different during the BDF period. With an approximately correct value of G (62.3 Bscf), J g becomes constant, at 1.17*10 -7 MMscf/d/(psi 2 /cp). The values of J g are on increasing and decreasing trends when G is too small and too big, respectively. The results are consistent with the previous studies (Agarwal et al. 1999;Araya and Ozkan 2002;Medeiros et al. 2010;Ismadi et al. 2011).
Since the proposed methodology is only valid for the BDF flow periods, production data during transient flow periods (the first 10 production days of the first and second years) are removed from the analysis. The versus-time profiles of J g for three different values of G are shown in Fig. 4. The values of J g are overestimated for G < 62.3 Bscf and are underestimated for G > 62.3 Bscf. For a given value of G, the profiles of J g during the first and second years are similar. Figure 5 shows the minimum, average, and maximum values of J g plotted as a function of the  assumed value of G. The minimum value of J g approaches the correct value when G is underestimated, and the maximum value of J g approaches the correct value when G is overestimated. The range of J g for G < 62.3 Bscf is significantly larger than for G > 62.3 Bscf. Figure 6 shows the plot of SD(J g ) versus G. The value of SD(J g ) is a minimum for a value of G of 62.3 Bscf. This value of G is about 0.8% larger than the input value in the model (61.78 Bscf). This is due to one of the fundamental assumptions in gas material balance equation. In Eq.
(2), the effects of connate water and formation compressibilities are ignored.
To account for these two factors, the gas material balance equation becomes (Dake 1978) The analysis with this new gas material balance equation yields a value of G of 61.6 Bscf, about 0.3% underestimated. The residual is attributed to numerical approximation. It should be noted that the values of SD(J g ) on the left (smaller) side of the correct value of G are more sensitive to the value of G than those on the right (larger) side. Figure 7 shows a plot of q g versus [m(p)−m(p wf )] for different values of G. These three profiles are straight lines which differ in the values of their y-intercepts. When the value of G is approximately correct, 62.3 Bscf, the y-intercept is close to zero. When the value of G is underestimated (overestimated), the sign of y-intercept is positive (negative). The value of J g , from the slope of the correct straight line, is 1.165*10 -7 MMscf/d/(psi 2 /cp). The values of b during BDF period are estimated from Eq. (10) and are plotted in Fig. 10. The value of b is declining during Year 1 when p wf is held constant at 3,500 psia. Then, p wf is reduced to 2,000 psi and is kept constant during Year 2. The value of b initially increases sharply and then declines. The results clearly show that the value of b is not constant, but rather is changing with pressure depletion. Therefore, application of Arps' model (Arps 1945) to analyze gas production data should be performed with caution since the model assumes that the value of b is constant, and that method can therefore yield optimistic results.
Next, the effect of including transient data points on the accuracy of the method is investigated. The results are plotted in Fig. 11. The parameter on the x-axis is defined as the number of data points in the transient flow period divided by the total number of data points in the analysis. Productivity 0.0E+00 1 .0E+08 2.0E+08 3 .0E+08 4.0E+08 5 .0E+08 6.0E+08 indices during the transient flow period are higher than those during the BDF period. As shown in the previous section, data points with higher productivity index correspond to lower values of G, and vice versa. This is the reason why including more data from the transient flow period will lower the estimated value of G. The values of G are 0.4%, 2.3%, and 6.2% underestimated when there are about 10%, 20%, and 30%, respectively, of transient data points in the analysis. The error profile is concave downwards. As more transient data are added into the analysis, the marginal effect on error in G estimation becomes more significant. Figure 11 suggests that small errors in estimating the time of the start of the BDF period should not cause significant error in the estimation of G.
To generate the production forecast accuracy, we use data from the first year for history matching, and forecast the second year of production during BDF period. There is a good match between production data and production forecast, Fig. 12. Note that the proposed model cannot forecast the production profile during the transient flow period. Now the dynamic material balance method (DMB) of Mattar and Anderson (2005) is applied to analyze the production ProducƟon data data in Example 1. For varying gas flow rates, the pseudoequivalent time (t a ) of Palacio and Blasingame (1993) is used to replace the pseudotime for gas (t ca ). Then, the relevant equation in the DMB method becomes.
Where (16) The value of b pss can be read at the y-intercept on the plot of [m(p i )−m(p wf )]/q versus t a . The above equation can be rearranged as Agarwal et al. (1999) showed that plotting [q/Δm(p)] vs [2p/(μc g z) i ]*[ q/Δm(p)]*t a yields a straight line with an x-intercept of G. The method proposed in this study yields a more accurate result, less than 1% error. In addition, the procedure and calculations in the proposed method are not iterative and are much easier to understand and apply than those of the other methods. The remaining limitations of the proposed method will be discussed later, in the discussion section.

Example 2: A gas well producing at two different flow rates
The reservoir is exactly the same as one in Example 1. The difference is the operating conditions. The production profile is shown in Fig. 13. Once the BDF period is reached, with a correct value of G, the value of [m(p)−m(p wf )] becomes stable or shows minimal variation. When the value of G is underestimated (overestimated), the value of [m(p)−m(p wf )] decreases (increases) with production time during the BDF period. Figure 16 shows profiles of J g as a function of t a for constant-p wf and constant-q g productions during the first year of production. There is good agreement between the two cases. It means that the profile of J g versus t a is independent of well operating conditions. The results are consistent with the findings of Medelros et al. (2010). As before, production data during the first year are used for history matching, and the production forecast is generated for the second year of production during the BDF period. There is a good match between production data and production forecast, Fig. 17.

Example 3: West Virginia hydraulically fractured gas well A
Well A is a gas well completed in the Onondaga chert in West Virginia. The reservoir and fluid properties are shown in Table 2. Hydraulic fracturing is essential for production since this reservoir has very low porosity and productivity. The well was on production for 200 days and shut in for a 106-day pressure build-up. The wellbore flowing pressure was kept at 500 psia. The well suffered from water load up at low gas rates. Production data and fluid properties are from Fraim and Wattenbarger (1987). Figure 18 shows the production profile. The profiles of J g for three different values of G are shown in Fig. 19. Figure 20 shows the plot of SD(J g ) versus G. The values of J g are clearly more sensitive when too-low values of G are input than when too-high values are input. The value of SD(J g ) is a minimum at a value of G of 2.7 Bscf. Figure 21 shows the plot of q g versus [m(p)−m(p wf )] for a value of G of 2.7 Bscf, with the best-fit straight line. The y-intercept is very close to zero. The value of J g , slope of the straight line, is 1.84*10 -6 Mscf/d/(psi 2 /cp). The results of other studies performed on this data set using alternative methods are reported in Table 3. The result of this study is consistent with the results from the previous studies. We perform history matching using production data during the first year, and generate the production forecast for the second-eighth year. As shown in Fig. 22, the forecast is consistent with actual production data.

Discussion
• The value of J g is assumed to remain constant throughout the production life (except during transient periods). However, reservoir productivity may decline with pres-  (2021) proposed an analytical model and a solution for the drawdowndependent well productivity index by accounting for the depletion-induced permeability variation for the steadystate regime of reservoir flow. The nth-order well productivity index can be expressed as where p * n (r w ) is the nth-order perturbation solution for pore pressure, expressed in Eq. (13) of their paper. To apply their model to a gas reservoir, the pressure must be replaced by pseudopressure. With these modifications, the proposed methodology can be extended to analyze production data from a gas system with varying productivity index.
• For a gas-condensate reservoir, productivity index varies with production after the reservoir pressure is depleted below the dew-point pressure. Heidari Sureshjani and Gerami (2011) proposed a simple model for a gas-condensate reservoir using flowing well data. The model can be applied to account for the variation in productivity index. The productivity index at the initial reservoir condition can be expressed as where two-phase pseudopressure (m tp ), two-phase material balance pseudotime (t acr,tp ), and two-phase pseudotime (t a,tp ) are defined as  For a gas-condensate system, Hagoort (1988) and Vo et al. (1990) recommended the following two-phase material balance equation.
Produced moles (n p ) and moles at initial condition (n i ) are defined in terms of the molar density of the gas component at standard conditions (ρ gsc ), the molar density of the oil With these modifications, the proposed methodology can be extended to analyze production data from a gas-condensate system.
• For unconventional oil/gas reservoirs, productivity index is decreasing during the transient-flow period and become a constant during BDF period (Al-Rbeawi, 2018). Therefore, the proposed method is still applicable but only for data during BDF period. However, these kinds of reservoirs normally have an extremely long transient-flow period depending on the size of simulated reservoir volume (SRV), connectivity of fracture net-  works, and matrix permeability. Therefore, the application both of the proposed method and of Arps' model is limited for such cases. As an alternative, Gupta et al. (2018) modified the Arps exponential-decline equation, where the constant decline rate is replaced by a powerlaw-function variable decline rate, and proposed the following variable-decline-modified Arps (VDMA) model: where Di is the initial decline rate and a is the decline-rate exponent. The relationship between instantaneous decline rate (D) and time (t), which can be used to estimate the two parameters, is They applied this model to interpret field data and found that the average error in the prediction of cumulativeproduction is less than 3%. This method can account for changes in fracture conductivity with time and for different flow regimes and is very robust. For a given abandonment pressure, coupling this model with the gas material balance equation allows a value of G to be obtained.
• The proposed methodology is applicable only for data in the BDF period. Therefore, the transient flow period must be identified and its data must be removed from the analysis. Including transient data tends to lead to underestimation of G. However, the profile of J g versus t a is independent of well operating conditions, as illustrated in Fig. 16. Therefore, we can use t a to link transient data for any operating condition to the equivalent data for a constant-flow-rate condition whose analytical solutions are well established in pressure transient analysis. The results are reservoir fluid and rock properties, not G.

Conclusions
• This study proposes a new simple and accurate methodology to estimate gas initially-in-place and gas productivity index by coupling a gas material balance equation and a stabilized gas flow equation. It requires only the flowing production profile and PVT data. The correct value of gas initially-in-place (G) is the one which yields the least variation in gas productivity index (J g ), and a plot of q g versus [m(p)−m(p wf )] yields a straight line with a y-intercept of zero. o Applicable only for data during the BDF period o Constant-J g assumption o Non-Darcy flow is not considered • The method has been successfully tested with 2 synthetic data sets and 1 field data set, with high accuracy results. • The method can potentially be extended to a reservoir with varying productivity index, a gas-condensate reservoir, and an unconventional reservoir by applying established methodologies to account for those conditions.

And substituting Eq. (B.5) in (B.3) yields
Funding The author(s) received no specific 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, 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/.