Integration of Multi-Region Material Balance Equation with Binomial Productivity Equation for Performance Prediction

In order to overcome the limitations of conventional material balance methods, a new multi-region material balance method, considering the strong heterogeneity and low-velocity non-Darcy ﬂ ow characteristics of bio-reef gas reservoir, is proposed to estimate the development performance of reef reservoirs. Firstly, the reef reservoir is divided into two regions, the reef cover high-permeability development area and the reef edge low-permeability supplying area. Then, based on one-dimensional non-linear seepage mechanism and mass conservation principle, the multi-region supplying material balance equation is established. Finally, combined with the binomial productivity equation, the development performance can be easily estimated at di ﬀ erent P / Z curves. The application results indicate that the low-permeability replenishment area plays an important role, especially after entering the decline stage, on the development performance of high-permeability region. Compared with the traditional material balance method, the new method can help for better understanding of unbalanced development in di ﬀ erent regions and the replenishment e ﬀ ect of low-permeability region.


Introduction
A typical gas well needs to go through three periods, including constant gas rate production, production decline, and pressurization production [1]. At the stage of constant rate production, traditional material balance, production dynamic fitting, and numerical simulation are three main methods, which are usually adopted to predict the dynamic of gas reservoirs [2]. However, these technologies are not applicable to bio-reef gas reservoirs for the heterogeneity and the low-velocity non-Darcy flow characteristics.
(1) Traditional material balance is a simple and convenient method to evaluate the dynamic reserves and only can carry out simple development performance prediction [3,4]. However, the application indicates that the traditional material balance cannot accomplish to accurately describe the production characteristics of different types of reservoirs (2) Numerical simulations are common methods used to analyze the dynamic characteristics [5,6]. However, long well spacing and sparse well pattern make it difficult to build accurate geological models, and it is inconvenient to be applied to evaluate the production characteristics of low-permeability reservoirs (3) Production performance fitting methods are another method to reveal the low-permeability reservoir seepage law [7,8]. It requires more detailed geological information and reservoir fluid parameters, which cannot be adopted to quickly evaluate the production characteristics of heterogeneous reservoirs based on constant rate production period data, but lacks a complete theoretical basis. Based on the concept of multiregion replenishment, Shiqing Cheng et al. [12][13][14][15][16] conducted an in-depth dynamic evaluation of the low-permeability gas reservoirs, but early predicting the dynamic development characteristics of a single well cannot realize. In addition, the MMBE, considering the low-velocity non-Darcy flow of lowpermeability reservoirs, was used to evaluate the elastic flooding and dissolved gas flooding recovery efficiency [17][18][19][20]. However, the reservoir is simplified to a number of annular bands which made the theoretical foundation imperfect.
In this work, a new MMBE method is developed to carry out gas well performance prediction. Compared with the previous work [9][10][11][12][13][14][15][16][17][18][19][20][21][22], this new method initially introduced the effects of low-velocity non-Darcy seepage and combined the MMBE with binomial productivity equation for development performance prediction. The specific steps of this method are as follows: firstly, based on the architectural model of reef body and low-velocity non-Darcy seepage law in low-permeability region, the heterogeneous bio-reef reservoirs could be simplified to 1D reservoir model with two regions. Then, considering the low-velocity non-Darcy seepage in low-permeability region, a new mathematic method derived from heterogeneous multi-region model is proposed to estimate the reservoir parameter. Subsequently, integrating the multi-region material balance equation with binomial productivity equation, an integrated approach assessing the gas well development performance is derived. Finally, practical application was conducted, and some important conclusions are obtained. The flow chart of using this method to carry out development performance prediction can be seen in Figure 1.

Reservoir Characteristics of Typical Bio-Reef Gas
Reservoirs. The bio-reef gas reservoir presents the characteristics of multi-stage sedimentation; the reef core and reef front are impervious regions. The reef cover with highpermeability is a development area, and the reef edge with low-permeability is a replenishment area, shown in Figure 2. The low-permeability region exhibits strong non-Darcy flow characteristics. As shown in Figure 3, the fluid seepage of lowpermeability reservoir is divided into two stages: (1) Affected by gas-solid interaction, the first stage shows low-velocity non-Darcy seepage characteristics, and the gas seepage velocity closes to 0. For the convenience of research, we assume that the seepage velocity is equal to 0, when the pressure gradient is lower than G, non-Darcy seepage factor, shown in Figure 3 and Eq. (1) (2) When the pressure gradient is higher than G, the second stage behaves linear Darcy flow. As shown in Figure 3 and Eq. (1), there is a positive correlation between seepage velocity and pressure gradient The reef configuration shows that the gas reservoir is divided into two parts: reef cover highpermeability development area and reef edge lowpermeability supply area, as shown in Figure 4. For the convenience of research, a typical one-dimensional two-zone seepage model is used to approximately describe the configuration characteristics of the typical reef. The assumptions of the model are as follows: (1) The reservoirs in high-permeability area and lowpermeability area are horizontally homogeneous with equal thickness, ignoring the influence of gravity and water (2) The reservoir in the low-permeability area is dense and has strong gas-solid interaction, showing the characteristics of low-velocity non-Darcy seepage, while the reservoir in the reef cover has highpermeability and the seepage meets the characteristics of Darcy seepage (3) The supply of the two sections is controlled by the pressure difference between the two regions, the supply gas volume is directly proportional to the pressure difference between the two areas, and the supply flow distance is equal to ðL 1 + L 2 Þ/2, shown in Figure 4 (4) The length of high-permeability area and lowpermeability area is L 1 and L 2 , respectively, the width is L, and the formation thickness is h (5) In the high-permeability area and the lowpermeability area, the geological reserves are V 1 and V 2 , respectively, and the permeability is K 1 and K 2 , respectively (6) The gas production rate in the high-permeability area is q g 2.2.2. Mathematical Model and Solution. Reservoir heterogeneity has a greater impact on the development of bio-reef gas reservoirs and fluid migration between different regions. The low-permeability area shows obvious low-velocity non-Darcy seepage characteristics, which makes the pressure drop in the reef edge low-permeability area slower than that in the reef cover high-permeability area, and the formation pressure difference between the two areas is obvious. Affected by the physical property difference between the two areas, the gas well production, in the high-permeability area, shows two obvious development stages. In the initial stage, only the fluid in the reef cover high-permeability development area flows to the near well region. When the pressure gradient between the two areas is greater than the 2 Geofluids low-velocity non-Darcy flow factor G, the fluid in the reef edge low-permeability area starts to replenish to the reef cover high-permeability area.
(1) Early seepage stage in high-permeability region Affected by the low-velocity non-Darcy seepage characteristics in the low-permeability area, and when the pressure gradient between the two areas is less than the low-velocity non-Darcy flow factor G. Therefore, only the fluid in the reef covers high-permeability area seepages, and the fluid in the low-permeability area does not flow.
Based on the principle of mass balance, the initial material balance equation in the early seepage stage can be obtained. For the specific formula derivation process, see Appendix A.
(2) Late replenishment stage in low-permeability region When the cumulative gas production volume N P is greater than N P0 , see Figure 3, the pressure gradient between the two areas is greater than the low-velocity non-Darcy flow factor G, and the fluid in the low-permeability area begins to supply to the high-permeability area.
Based on the principle of zonal replenishment, the material balance equations of different regions can be obtained. For the specific formula derivation process, see Appendix B.    Figure 3: The relationship curve between seepage velocity and pressure gradient.

Geofluids
2.3. Dynamic Prediction Method. Assuming that the development is divided into two development stages, constant rate production and decreasing rate production, and there is only one well in the high-permeability area. After productivity testing, we can get the binomial productivity equation of the gas well: where P wf is bottom hole flowing pressure, MPa; a and b are the binomial capacity equation coefficients. By synthesizing binomial production equation and multi-region material balance equation, it can be able to predict the development performance of gas wells.
(1) Calculation of constant rate production period (1) Calculation of high-permeability development zone Suppose a single well is initially produced at a constant rate q g until the bottom hole pressure drops to the lowest value P wL . According to Eq. (5), the reservoir pressure P s in the high-permeability development zone at the end of the constant rate production period is determined by the following equation: The reservoir pressure interval [P 0 , P S ] is divided into uniform m 1 intervals with node P j (j =0, 1, 2, ..., m 1 ). The production rate q g and pressure P j , in the constant rate production period, are brought into Eq. (2) or Eq. (3) to obtain the cumulative production time t, the cumulative gas production volume G pj , and the apparent formation pressure (P 1 , j /Z 1 , j ) of high-permeability area.
And the remaining reserves in the high-permeability development zone can be obtained through the following formula: (2) Calculation of low-permeability recharge area Furthermore, the cumulative gas production G pj , the cumulative production time t, and the gas production rate q g , in the constant rate production period, are brought into Eq. (4). And we can get the apparent formation pressure (P 2 , j /Z 2 , j ) of low-permeability area.
According to the following Eq. (8), the cumulative gas volume supply from the low-permeability area to the highpermeability development area is: In addition, the remaining reserves in the lowpermeability development zone can be obtained by the following formula: (2) Calculation of production rate decline period (1) Calculation of high-permeability development zone Assuming that the economic limit production is q L , and the abandoned reservoir pressure P L can be obtained by the following formula: Furthermore, the reservoir pressure and abandoned formation pressure interval [P s , P L ] are divided into any m 2 equal cells, and the node pressure is P j (j = m 1 + 1, m 1 + 2 ⋯ , m 1 + m 2 ). Similar to the calculation of constant rate production period, the average reservoir pressure and average production in interval j can be obtained, and the relationship between those parameters is determined by Eq. (5), that is.
Solving Eq. (11) can be obtained: Incorporating the average production rate q g and node reservoir pressure P j into Eq. (2), we can get the cumulative production time t, the cumulative gas production volume G pj , and the apparent formation pressure (P 1 , j /Z 1 , j ) of high-permeability area. And the remaining reserves in the high-permeability development zone can be obtained through the following formula: 4 Geofluids (2) Calculation of low-permeability recharge area Furthermore, cumulative gas production volume G pj , cumulative production time t, and average production rate q g are taken into Eq. (3), and the apparent formation pressure (P 2 , j /Z 2 , j ) of the low-permeability replenishment area can be obtained. The cumulative gas replenishment volume from the low-permeability replenishment area to the highpermeability development area can be obtained by the following equation: In addition, the remaining reserves in the lowpermeability development zone can be obtained by the following formula: 2.4. The Limitations of This Method. Take into account that the physical and mathematical models have many presuppositions, the proposed method has the following three limitations: (1) Reservoirs in the high-and low-permeability zones are simplified into two independent homogeneous whole, which cannot reflect the internal seepage characteristics of the two areas. And the evaluation parameters can only reflect the overall characteristics of each seepage unit (2) This method can only evaluate the reserves in the low-permeability area and the average distance from the low-permeability area to the high-permeability area. And the distribution characteristics of lowpermeability reserves can only be determined by combining geological analysis (3) As the binomial productivity equation will change with the formation pressure drop, introducing the initial binomial productivity equation to predict the development performance has relatively low accuracy

P/ZCurve
Characteristics. Based on the above proposed new multi-region material balance, the relationship curve between the apparent reservoir pressure ðP/ZÞ and the cumulative gas production volume G P in the different areas can be obtained, shown in Figures 5 and 6. Those curves illustrate obvious two-stage segmentation characteristics.
(1) In the first stage, only the fluid, in the highpermeability zone, seepage to the bottom of the gas well. In the high-permeability development zone, the first decline characteristic of the curve is mainly affected by the reserves in the high-permeability development zone, and the decline rate is equal to ðP/ZÞ 1 /V 1 . In the low-permeability development zone, the P/Z curve value is constant as the initial P/Z value (2) At the second stage, the reserves in the lowpermeability area began to be recharged to the high-permeability area. Affected by the reserves in the low-permeability recharge zone, the decline rate closes to ðP/ZÞ 1 /ðV 1 + V 2 Þ in the high-permeability zone. Compared with the traditional material balance method, the recoverable reserves evaluation result of this model is significantly lower, and the difference is ΔG P , as shown in Figure 5. And the decline rate of low-permeability zone is equal to ðP/ZÞ 2 /V 2

Practical Application.
Taking well W of a bio-reef gas reservoir as an example, the initial reservoir pressure is 68.69 MPa, and the average daily production rate is Abandoned reservoir pressure G P P ab v 2 P 1 i Z Figure 5: P/Z curve of typical multi-region material balance equation (high-permeability zone).
Abandoned reservoir pressure P ab 0 Figure 6: P/Z curve of typical multi-region material balance equation (low-permeability zone). 5 Geofluids 470,000 m 3 /d. After productivity test, the open flow rate is 1,990,000m 3 /d, and the binomial productivity equation is: (1) P/Z curve fitting Figure 7 illustrates the fitting results of this well and shows the relationship curve between P/Z and cumulative gas production volume in the high-permeability development area of this well. Applying the new multi-region material balance method to fit the P/Z values of different times, we can get relevant reservoir parameters of different regions. The evaluation results show that the reserves in the highpermeability development zone are equal to 30 × 10 8 m 3 , the reserves in the low-permeability recharge zone are 10 × 10 8 m 3 , the low-velocity non-Darcy seepage factor G in the low-permeability recharge zone is 0.008Mpa/m, and the seepage resistance coefficient A is equal to 0.2.
As shown in Figure 8, the cumulative gas volume supplied from low-permeability region behaves a linear increase trend. When the cumulative gas production reaches 1:94 × 10 8 m 3 , the fluid in the low-permeability zone begins to supply to high-permeability. When the gas well produced 1 × 10 8 m 3 natural gas, the low-permeability area could supply 0:26 × 10 8 m 3 natural gas to the high-permeability area.
(2) Gas well performance prediction The lowest bottom hole flow pressure P wf of this well is 16 MPa, the economic limit production q L is 2 × 10 4 m 3 /d, and the gas rate of constant rate production period is 47 × 10 4 m 3 /d. Based on the new multi-region material balance and traditional material balance methods, we can get the dynamic forecast curve of this well, as shown in Figure 9. Using the new multi-region material balance equation to predict the development performance, the constant rate production time is 10.3 years, which is 0.65 years lower than prediction results of the traditional material balance method.
According to the new multi-region material balance method proposed in this work, the change law of average formation pressure and supplying gas rate can be obtained, as shown in Figure 10. The curve characteristics can be divided into three stages: (1) First is the low-velocity non-Darcy seepage effect stage. The reservoir pressure P 1 and bottom hole flowing pressure P wf , in the high-permeability development area, show rapid decline. However, affected by the low-velocity non-Darcy seepage law, the reservoir pressure (P 2 ) in the low-permeability recharge area remains constant, and the gas rate supplied from low-permeability region is equal to 0 (2) Second is the later stage of constant rate production period. In this period, the pressure gradient between the low-permeability and the high-permeability area is greater than the low-velocity non-Darcy seepage factor G. With the cumulative gas volume increasing, the reservoir pressure (P 1 ) in the high-permeability area, the reservoir pressure (P 2 ) in the lowpermeability area, the bottom hole pressure, and the cumulative gas volume (Q C ), supplying from the low-permeability zone to the high-permeability development zone, illustrates linear rapidly changing characteristics (3) Third is the stage of decline in production. During the production decline stage, the gas production rate decreases rapidly. The reservoir pressure in highpermeability area, the reservoir pressure in lowpermeability area, and the bottom hole flowing pressure show a slow decline. At the same time, the cumulative gas volume supplying from the lowpermeability zone to the high-permeability development zone shows a slow increase In summary, as production progresses, the lowpermeability replenishment area will gradually play an important role. Especially after entering the decline stage, Actual data Fit data Figure 7: P/Z data fitting curve of W gas well.  6 Geofluids the proportion of the supplied gas in the cumulative production gas volume is even greater. In addition, the total cumulative gas production volume is estimated to be 25:3 × 10 8 m 3 , and the final recovery degree is 58.8%.

Summary and Conclusions
In this work, a new multi-region material balance method for dynamic prediction of bio-reef gas reservoir is presented, which comprehensively considers the heterogeneity of reefs, low-permeability reservoir supply, and initially introduced the effects of low-velocity non-Darcy seepage. And this method further combined the MMBE with binomial productivity equation for development performance prediction. In addition, practical application was also conducted. Then, some important conclusions are drawn as following: (1) Through fitting the P/Z curve of typical bio-reef gas well, dynamic reserves in high/low-permeability zones, the resistance coefficient, and the non-Darcy seepage factor in the low-permeability zone can be obtained. Compared with the traditional material balance method, this method can more precisely describe the development dynamics of the low-permeability zone (2) Combined the new multi-region material balance with productivity equations, we can achieve early dynamic prediction. The high/low-permeability zone reservoir pressure, bottom hole flowing pressure, and the recharge gas volume can be achieved by P/Z curves features in the constant rate production period   Gas production rate (ten thousand cubic meters/day) cumulative gas production G P (100 million squares) Figure 10: Dynamic prediction curve of multi-region material balance method. 7 Geofluids time is relatively low. It can more accurately predict the reserves of the low-permeability zone, providing a theoretical basis for the determination of the production well working system, and the deployment/adjustment of the development well pattern Appendix A. Early seepage stage in highpermeability region At the early seepage stage in high-permeability region, the material balance equation of the high-permeability development zone is: Now with the definition of gas compressibility: The chain rule for derivatives is: For the spatial derivative of ρ across the physical model of average thickness L, Eq. (2) can be rewritten as: Initial conditions: The results for the seepage in the high-permeability are: And the definition of cumulative gas produced: N P = q sc t: ðA:8Þ Eq. (A.6) can be written in the form: Where: where M is molecular weight; R is universal gas constant, MPa.m 3 /mol.K; P is reservoir pressure, MPa; P i is initial reservoir pressure, MPa; Z is gas deviation factor, decimal; P sc is pressure under standard conditions, MPa; T is the reservoir temperature, K; T sc is the temperature under standard conditions, K; ρ 1 is the fluid density in the high-permeability zone, kg/m 3 ; N P0 is the cumulative gas production volume starting to supply from the low-permeability area to the high-permeability area, 10 8 m 3 .

B. Late replenishment stage in lowpermeability region
Based on the one-dimensional non-linear equivalent analytical model of multi-region material balance, the material balance equation, at the late replenishment stage in lowpermeability region, can be described approximately by the pair of equations: sssHere Eqs. (A.11) and (A.12) describe the material balance in the high-permeability zone and the low-permeability zone separately.where ρ 1 is the fluid density in the highpermeability area, kg/m 3 ; ρ 2 is the fluid density in the lowpermeability area, kg/m 3 ; K * is the equivalent permeability between low-permeability area and high-permeability area, mD; G is the low-velocity non-Darcy seepage factor, MPa/ m; t is time, day; ρ sc is the fluid density under standard conditions, kg/m 3 ; q sc is the gas production rate under standard conditions, 10 4 m 3 /day.
Initial conditions: where ρ 0 is the gas density under the original conditions, kg/m 3 . Eq. (A.11)-Eq. (A.14) can be transformed into a secondorder homogeneous linear differential equation with constant coefficients: Among them, the parameters A, B, and γ are, respectively, expressed as: : ðA:16Þ where C g is the natural gas compression factor, MPa -1 ; A is the seepage resistance coefficient, decimal. From the above formula, we can get: Subscripts i: Gas reservoir initial conditions SC: The standard conditions 1: The high-permeability zone 2: The low-permeability zone g: Natural gas wf : Bottom hole P0: Starting to supply gas from the low-permeability area to the high-permeability area, 10 8 m 3 .
Superscripts * : the equivalent parameter between low-permeability area and high-permeability area.

Data Availability
The [Bio-reef sedimentary architecture map, pictures, pressure data] data used to support the findings of this study are included within the article.