Two-phase Flow Model and Productivity Evaluation of Gas and Water for Dual-Medium Carbonate Gas Reservoirs

Fluid flow in the dual-medium carbonate gas reservoir is characterized by stress sensitivity and non-Darcy flow effect. In order to accurately describe the unsteady flow of gas and water in the dual-medium gas reservoir, a two-phase flow model of gas and water is built. First, reservoir space and fluid flow characteristics of carbonate gas reservoirs are investigated, and the flow model that considers both the stress sensitivity and non-Darcy flow is built based on the fundamental flow theory, after fully investigating the reservoir space and fluid flow characteristics of carbonate gas reservoirs. Then, the perturbation theory is introduced, and the model is solved in the Laplace space, after which the obtained Laplace space analytical solution is converted into the real-space solution. Finally, the productivity evaluation model for the dual-medium gas reservoir with the gas-water two-phase flow is built, based on the flowing material balance method and Newton iteration. The presented productivity evaluation model is applied to analyze the effects of stress sensitivity and non-Darcy flow on the two-phase flow model of gas and water for the dual-medium gas reservoir and the reservoir productivity. The results indicate that a higher stress sensitivity coefficient is demonstrated to indicate higher stress sensitivity and accelerated production decline of the reservoir, while a lower non-Darcy flow effect coefficient represents a stronger non-Darcy effect and boosted drop of initial production of the reservoir. Hence, it is not reasonable to neglect the effects of stress sensitivity and non-Darcy flow during the evaluation of the productivity of a dual-medium carbonate gas reservoir. The model presented in this research provides important references for improving the recovery performance of dual-medium gas reservoirs.


INTRODUCTION
Carbonate reservoirs are important targets for hydrocarbon exploration and exploitation, because they contain abundant oil and gas resources (Xu et al., 2007). Compared with conventional gas reservoirs, the dual-medium carbonate gas reservoir is associated with greatly different reservoir characteristics that may result from the relatively developed pores and fractures. Its flow is typically characterized by stress sensitivity, non-Darcy effect, and unique relative gas-water permeability relationship.
In terms of stress sensitivity, Fan et al. (2011) develop a flow model incorporating the stress-sensitive permeability of fracture systems via analyzing the medium deformation of the dualmedium gas reservoir. Zhao et al. (2013) carry out stress sensitivity experiments using four types of cores, namely the matrix type, the fully-filled-fracture type, the semi-filled-fracture type, and the non-filled-fracture type. They rank the stress sensitivity of these four types of cores in the increasing order of matrix, fully-filled-fracture, semi-filled-fracture, and nonfilled-fracture types. Moreover, they identify a power-law relationship between permeability and stress. Peng et al. (2015) experimentally analyzed the stress sensitivity of cores under the constant confining pressure and varied pressure and its effects on the reservoir productivity. Zhang R et al. (2016) compare five frequently-used formulas for stress sensitivity via an experimental-theoretical-integrated approach and conclude that the stress sensitivity of fractured reservoirs can be characterized using the power-law formula. Luo et al. (2021) considered the stress-sensitive effects of natural fractures based on the dualmedium flow theory model.
Also, studies on the characteristics of non-Darcy flow have been reported. Cheng and Chen (1998) probe the characteristics of low-speed non-Darcy flow in the case of two-phase flow of oil and gas and point out that non-Darcy flow to some extent impacts (reduces) oil production. Moreover, with the same water content, the productivity index of non-Darcy flow is smaller than that of Darcy flow. Lu (2010) experimentally investigates the fluid flow characteristics under the varied fracture and fracture-vug conditions and identifies the critical parameter values for non-Dary flow, and builds the flow model by introducing Forchheimer number to characterize non-Darcy flow. Zhang F et al. (2016) study the effects of high-speed non-Darcy flow on the productivity of gas wells with certain water cut, via theoretical derivation. Furthermore, the analysis based on the derived productivity equation and actual data of production wells demonstrates that the productivity of gas wells with a certain water cut is considerably lower than that of gas wells with no water production; the non-Darcy effect of productivity of gas wells declines after water production starts. Javadpour et al. (2021) considered the non-Darcy flow of shale gas, and reviewed the dominant gas-flow processes in a single nanopore based on theoretical models and molecular dynamics simulations, and Lattice Boltzmann modeling.
As for the relative permeability of gas and water, Esmaeili et al. (2020) use a newmethodology and conduct limited steady-state relative permeability measurements at different temperatures to confirm the validity of displacement-based relative permeability. Fang et al. (2015) study the pattern of the relative permeability curves of gas and water in high-temperature, high-pressure tight gas reservoirs. Their research identifies that gas permeability under high temperature and high pressure is found to be higher than that under conventional conditions, which means high temperature and high pressure are in favor of gas flow. Liao (2016) investigates the characteristics of the two-phase flow of gas and water in the fractured gas reservoir using numerical models. In his research, he discussed the effects of numerous factors on the production characteristics of the gas reservoir based on the mechanism model, including geological parameters, aquifer parameters, production schemes, water intrusion patterns, and heterogeneity.
At last, for the flow model and productivity evaluation model of gas reservoirs, Zhang et al. (2017a) combine the material balance equation and stress-dependent production equation in the case of over-pressurization to analyze the stress sensitivity of formations and derive the productivity prediction model. The result indicate the critical role of stress sensitivity in reducing the production of gas wells and shortening the duration of stable production. Brown et al. (2009) analyze the pressure and productivity characteristics of stage-fractured horizontal wells of shale gas reservoirs using an analytical tri-linear flow model. Deng et al. (2011), based on their analysis of the gas-water relative permeability regularity, modify the productivity prediction model of water-producing gas wells and the production decline analysis method for water-involved gas reservoirs with the non-Darcy flow. Chen (2016) probes the fluid flow patterns of lowpermeability gas reservoirs and develops the productivity calculation formulas for fractured horizontal wells and those with inclined fractures and gas-water two-phase flow, respectively. The presented models are validated by the good application performance of the models to the production data of actual production wells. In 2017, a productivity evaluation method for abnormally over-pressurized gas reservoirs, incorporating both the reservoir rock deformation and gaswater relative permeability, is developed by Zhang et al. (2017b), which is then applied to investigating the factors affecting the production performance of gas wells via actual production data. Based on the material balance theory, Zhang W et al. (2017) derive the pressure calculation formula and production prediction equation for gas wells featuring twophase flow, which provide theoretical support to production forecasting of gas wells with coexisting gas and water.
The above review demonstrates substantial progress in studying the flow regularity and productivity prediction of carbonate gas reservoirs. Nonetheless, most research only considers the effects of a single factor, and rare studies have been reported to investigate the joint impacts of stress sensitivity, non-Darcy flow, and two-phase relative permeability of gas and water on carbonate gas reservoirs. Moreover, complete, systematic productivity evaluation methods for carbonate gas reservoirs haven't been presented yet. Given these, this study aims and manages to build a gas-water two-phase flow model of carbonate gas reservoirs, incorporating both stress sensitivity and non-Darcy flow. Results of this study are expected to provide references for improving the recovery performance of dual-medium carbonate gas reservoirs.

FLUID FLOW PATTERNS FOR CARBONATE GAS RESERVOIRS
Dual-medium carbonate gas reservoirs are different from conventional reservoirs, and they are composed of matrix and fracture systems -the low-porosity, high-permeability fracture system is the main flow channel for fluids, while the high-porosity, low-permeability matrix system is the main storage space for gas. Such dual-medium carbonate reservoirs have greatly different pore-throat structures and more complicated fluid flow than conventional pore-dominated gas reservoirs (Jia et al., 2013).

Stress Sensitivity of Reservoirs
During the recovery process of carbonate gas reservoirs, the reservoir rock framework deforms, as the effective stress grows (e.g., the overburden stress) owing to the sustained production of natural gas and consequent pore pressure decline . The ultimate resultant variation of the reservoir physical parameters such as permeability and porosity is referred to as stress sensitivity (Peng et al., 2015). It is known that with the progress in hydrocarbon recovery, the reservoir permeability and productivity are both degraded.
The core porosity and permeability measurements under the overburden pressure of the L carbonate gas reservoir of the A gas field ( Figure 1) shows that the tendencies of permeability to decline with increasing overburden pressure are generally consistent; as the overburden pressure reaches 30 MPa, the core permeability declines by 60%; as the overburden pressure restores to the initial value, the permeability recovers to about 60% of the initial permeability (in other words, a permanent permeability loss of 40% is caused). The dimensionless permeability in Figure 1 is the ratio of permeability under certain pressure conditions to initial permeability (k i /k 0 ).
The stress sensitivity can be characterized using the power-law function (Petrosa, 1986): where k f is the permeability of the fracture system, mD; k fi is the fracture system permeability corresponding to the original formation conditions, mD; γ f is the permeability modulus, MPa −1 ; ψ i is the initial pseudo-pressure of formations, MPa 2 / (mPa·s); ψ f is the pseudo-pressure of the fracture system, MPa 2 / (mPa·s).
The pseudo-permeability modulus is defined as:

Effects of High-Speed Non-Darcy Flow
Due to the reservoir rock property, the dual-medium carbonate reservoir is subjected to considerable gas turbulence and thus the intensive non-Darcy effect (Chen et al., 2015). The non-Darcy flow regularity of the fracture-pore type reservoir rock of the L gas reservoir in the A gas field is illustrated in Figure 2. When the flow rate grows, the non-linearity between the pressure square difference and the gas flow rate gradually climbs up. Gas wells, subjected to the non-Darcy effect, are expected to present a certain degree of production decline. The non-Darcy effect can be characterized by introducing a correction factor dependent on fluid flow rates (Huang et al., 2004).
FIGURE 1 | Relationship between core permeability and overburden pressure.
FIGURE 2 | Non-Darcy feature of the fracture-pore type reservoir rock.
FIGURE 3 | Two-phase relative permeability of gas and water for the dual-medium.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 823764 where μ g is the fluid viscosity, mPa·s; u g is the flow rate, cm/s; δ is the non-Darcy coefficient, cm −1 ; ρ is the fluid density, kg/m 3 .

Two-phase Relative Permeability Pattern of Gas and Water
Because of the presence of dual media composed of pores and fractures and the substantial difference of fluid flow space between them, the fracture-pore type carbonate reservoir presents highly complex flow characteristics. Accordingly, the two-phase relative permeability pattern of gas and water is also greatly different from that in the conventional gas reservoir ( Figure 3). Cores in the study area are proven to be highly hydrophilic, with the relative permeability intersection point occurring at the water saturation of above 70%. This is caused by the micropore structure of the gas reservoir-some tiny throats in the pore structure are too small to allow water to form a continuous flowing phase after penetrating such space.

TWO-PHASE PRODUCTIVITY EVALUATION MODEL FOR GAS AND WATER IN DUAL-MEDIUM GAS RESERVOIRS Flow Model
The fracture-pore type dual-medium has been demonstrated to primarily consist of two types of pore structures, namely primary inter-granular pores and secondary pores. Primary pores are largely affected by deposition and diagenesis processes, dimensions, distributions, and shapes of reservoir rock grains, thereby presenting heterogeneous features of permeability and porosity. As for secondary pores, they are composed of fractures, joints, and dissolved pores. In other words, the dual-medium consists of two systems, namely pore and fracture systems. In most cases, the fracture system serves as the main fluid flow channel, while the rock matrix provides the main storage space. The fractured reservoir has the following characteristics: 1) Dual structural characteristics. The reservoir rock is associated with the extensive development of both primary and secondary pores (in other words, it consists of dual media of pores and fractures, respectively). Both the pore and fracture media have their own porosity and permeability features and the two media are connected.
2) The compressibility of the fracture system is considerably higher than that of the pore system. Therefore, as the fluid (pore) pressure decreases, the fracture may close under the overburden pressure; moreover, such fracture closure is irreversible. Ultimately, the fracture permeability is reduced. 3) Anisotropy. The reservoir permeability is highly anisotropic.
For productivity evaluation, if the goal is only to investigate the formation characteristics in a vertical well, the medium can be considered isotropic.
The fracture-pore medium model assumes that the reservoir is a dual-medium for fluid storage and flow composed of the matrix and fracture systems. The assumptions of the physical model for production decline of wells in the fracture-pore type reservoir are summarized below: 1) The vertical well is producing at a constant production rate from a reservoir with the closed outer boundary. 2) The payzone is producing across its whole thickness and the fluid flows into the wellbore via the radial flow regime.
3) The rock compressibility is neglected, while the gas compressibility and stress-dependency of the fracture system permeability are considered. 4) Fluids have two components, namely gas and water. 5) The effects of capillary pressure and gravity are ignored. 6) Fractures are the only channel for fluids to flow into the wellbore. Fluids in the matrix system flow into fractures via the pseudo-steady cross flow and then enter the wellbore through fractures.
The mathematic model of the fracture system is described below: The continuity equation: where ρ g is the gas density, kg/m 3 ; v g is the gas flow rate, cm/s; ρ w is the water density, kg/m 3 ; v w is the water flow rate cm/s; S g is the gas saturation, decimals; S w is the water saturation, decimals. Equations of motion: where k rg is the relative permeability of gas, decimals; μ g is the gas viscosity, mPas; B g is the formation volume coefficient of gas, decimals.
v w −δ kk rw μ w B w zp zr (6) where k rw is the relative permeability of water, decimals; μ w is the water viscosity, mPa·s; B w is the formation volume coefficient of water, decimals. The cross-flow rate between the matrix and fracture systems can be determined by: By combing the above equations, basic control equations of flow for the fracture system can be obtained: where k rgm is the relative permeability of gas for the matrix system, decimals; k rgf is that for the fracture system, decimals.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 823764 where k rwm is the relative permeability of water for the matrix system, decimals; k rwf is that for the fracture system, decimals.
Similarly, the basic control equations of flow for the matrix system can be obtained: The pseudo-pressure is defined as: where λ g 1 μ g B g . Thus, the basic control equation for gas flow in the fracture system is 1 r The basic control equation for gas flow in the matrix system is The initial condition is: The inner boundary condition is: The outer boundary condition is:

Solving the Model
The duration for productivity prediction is divided into numerous time intervals. The parameters dependent on saturations are calculated explicitly, which means that within a time step, the saturation remains the same. In contrast, parameters dependent on pressure are determined implicitly. The gas production is determined by solving our model, and meanwhile, the water production is predicted using the gas/water ratio. In other words, the model solving involves only the gas equation. After the gas production is determined, the average pressure and saturation across the production-affected (pressure drawdown) zone of the gas reservoir are computed via the flowing material balance method, according to which the model parameters are updated. Finally, the semi-analytical solutions of the model are obtained after numerous iterations. The solving workflow is summarized below:

1) Parameter normalization and equation linearization
The parameter normalization is shown below: The dimensionless stress sensitivity coefficient is defined as: The dimensionless time: The dimensionless radius: The dimensionless pseudo-pressure: Accordingly, the fracture control equation can be rewritten as: The matrix control equation can be rewritten as: The rewritten expression of the initial condition: The inner boundary condition: Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 823764 e −γ mD ψ fD δ zψ fD zr D rD 1 −1 (28) The outer boundary condition: The perturbation theory is introduced: The zeroth-order solution for perturbation is adopted and simplified as: (33) Correspondingly, the initial condition: The inner boundary condition: The outer boundary condition: zη 0fD zr D rD reD 0 (36)

2) Laplace transform and solutions of the mathematic model
For the convenience of derivation, we perform the Laplace transform to solve the model in the Laplace space. Subsequently, the Stehfest numerical inversion method is applied to obtain the solutions of the real space of the model.
The equation group is re-arranged using the Laplace transform of t D : The initial condition: The inner boundary condition: The outer boundary condition: Combining Eqs 37, 38 yields: Substituting the inner and outer boundary conditions into Eq. 42, the solution is obtained: The production is expressed as below: 3) Calculating the average formation pressure using the material balance method We manage to linearize the two-phase equation by introducing the pseudo-time and pressure. However, the obtained general solution is a function of the pseudo-time, while the pseudo-time is a function of pressure. Hence, the iterative calculation is required to compute the productivity corresponding to the real time and its composition (gas and water production). In iteration, the average gas reservoir pressure needs to be computed at each time step.
The overall material balance equation of the gas system is presented below: where IGIP is the initial gas in place, m 3 ; RGIP is the remaining gas in place, m 3 ; G p is the cumulative gas production, m 3 . Moreover, the IGIP can be calculated using the following equation:  (48) where r ginv is the drainage radius of gas, m; S gi is the initial gas saturation, decimals; B gi is the initial formation volume coefficient of gas, decimals. The RGIP is determined as below: where S g is the average gas saturation, decimals; B g 为is the average formation volume coefficient of gas, decimals.
The cumulative gas production: where q g is the daily gas production, m 3 /d. By substituting Eqs 48-50 into Eq. 47, we have: The overall material balance equation of water is: where IWIP is the initial water in place, m 3 ; RWIP is the remaining water in place, m 3 ; N w is the cumulative water production, m 3 . The IWIP can be calculated using the following equation: where r winv is the drainage radius of water, m; S wi is the initial water saturation, decimals; B wi is the initial formation volume coefficient of water, decimals. The RWIP can be calculated as below: where S w is the average water saturation, decimals; B w is the average formation volume coefficient of water, decimals.
The cumulative water production: where q w is the daily water production, m 3 /d. Substituting Eqs 53-55 into Eq. 50 yields: FIGURE 4 | Productivity prediction workflow based on the semi-analytic model for two-phase flow of gas and water.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 823764 S w B w The correlation between the water and gas saturation is shown below: By combining Eqs 51, 56, 57, we have: The derivative of Eq. 58 is: Then the Newton iteration equation can be written as Eq. 62: The two-phase flow model of gas and water is solved in a semianalytical approach. The productivity prediction workflow for the two-phase flow of gas and water is illustrated in Figure 4. The duration for productivity prediction is divided into numerous time steps. The saturation-dependent parameters are calculated explicitly, while the pressure-dependent ones are computed implicitly. Specifically, within one time step, the saturation is defined as a constant, and the pressure-dependent parameter is calculated using the pseudo-pressure. After obtaining the gas production, the average reservoir pressure and saturation across the pressure drawdownaffected zone are determined using the flowing material balance method. Subsequently, the model parameters are updated according to the new average reservoir pressure and saturation and the production is calculated again. The semi-analytical solutions of the model are finally obtained after multiple iterations.

ANALYSIS OF FACTORS AFFECTING PRODUCTIVITY OF CARBONATE GAS RESERVOIRS
The value of the permeability modulus γ represents how much the production performance of gas wells is affected by stress sensitivity ( Figure 5). A higher permeability modulus indicates higher stress sensitivity of the reservoir, which is associated with lowering down of the type curve and faster production decline. The dual-medium carbonate reservoir is found with extensive development of fractures. During the recovery process, the increase of the effective overburden pressure results in deformation of the rock framework and decline of      permeability, accelerating production decline. Moreover, with the elapsed time of recovery, the effect of stress sensitivity grows. Hence, the recovery of carbonate gas reservoirs shall adopt a proper recovery rate and drawdown pressure to minimize the impact of stress sensitivity and improve recovery performance.
The high-speed non-Darcy coefficient δ represents the effect of the high-speed non-Darcy flow on the type curve ( Figure 6). With a smaller value of the high-speed non-Darcy coefficient, the non-Darcy effect grows, the type curve is lowered, and the initial production is reduced. Therefore, the productivity evaluation of the dual-medium carbonate gas reservoir cannot neglect the non-Darcy effect.

CASE STUDY
Here the L carbonate gas reservoir of the A gas field is taken as an example and the developed two-phase flow model of gas and water based on the dual-medium is applied to productivity prediction of the production well, Well M2. This well presents an initial gas production of 17 × 104 m 3 /d during the production testing, and after 1 month of production, the production is raised to 28 × 104 m 3 /d. After 3 years of production, the daily gas production of this well is 9.6 × 104 m 3 /d, the daily water is 9.8 m 3 /d, and the water/gas ratio is 0.98 m 3 /10 4 m 3 . The basic parameters of this well are summarized in Table 1.
The developed productivity evaluation model is used for fitting the interpretation of the reservoir drilled by this well, based on its production performance. The fitting interpretation results are concluded in Figure 7; Table 2.
The future production of Well M2 at a constant flowing bottomhole pressure is forecasted, based on the model parameters interpreted from its actual production performance. The results are summarized in Figure 8; Table 3.
The error analysis of the predicted production shows that on Day 800, the predicted daily production of Well M2 is 9.47 × 104 m 3 /d, while the actual production rate is 9.61 × 104 m 3 /d. The error of the calculation is 1.5%, which is within the allowed error range and thus means that the productivity prediction results meet the engineering criterion.

CONCLUSION
1) The dual-medium carbonate gas reservoir is grealtly different from the conventional gas reservoir, in terms of the storage space and fluid flow characteristics. Therefore, the flow model shall consider both the stress sensitivity and high-speed non-Darcy flow effect. 2) The stress sensitivity and non-Darcy effect coefficients are introduced into the flow model of the dual-medium gas reservoir, which results in the nonlinearity of the flow equations. Then, the perturbation theory is applied to eliminate the nonlinearity and the flow equations are transformed into the Laplace space for solutions. The solution of fluid production is obtained via the Stehfest numerical inversion, and the semianalytical solutions of the model are calculated via the flowing material balance method and Newton iteration.
3) The analysis of the production decline type curve of the dualmedium carbonate gas reservoir shows that the dual-medium gas reservoir with higher stress sensitivity is associated with accelerated production decline; that with a more intensive non-Darcy flow effect is found with lower initial production. 4) The developed model is used for the case study of an actual production well. The analysis results show that the error of the predicted production of the analyzed well is within the engineering-permissible range, which validates the high applicability and practicability of our model.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
LK, WG, XZ, and YL contributed to design of the study and method establishment, wrote the first draft of the manuscript. JG, ZS, ML, and YS wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
The work was funded by R&D Department of Petrochina (No. 2021DJ2 005). The APC was funded by R&D Department of Petrochina (No. 2021DJ2 005). The funders had no role in the design, execution, interpretation, or writing of the study.