Influence of Stress Sensitivity on Water-Gas Flow in Carbonate Rocks

Carbonate reservoirs signi ﬁ cantly contribute to exploitation. Due to their strong heterogeneity, it is of great signi ﬁ cance to study core seepage capacity and gas-water two-phase ﬂ ow of reservoirs with various pore structures under di ﬀ erent stresses for productivity prediction, gas reservoir development, and reservoir protection. We utilize micrometer-resolution X-ray tomography to obtain the digital rocks of porous, fractured-porous, and fractured-vuggy carbonate rocks during pressurized process and depressurization. The Lattice Boltzmann method and pore network model are used to simulate the permeability and gas-water two-phase ﬂ ow under di ﬀ erent con ﬁ ning pressures. We show that at the early stage of pressure increase, fractures, vugs, or large pores as the main ﬂ ow channels ﬁ rst undergo compaction deformation, and the permeability decreases obviously. Then, many isolated small pores are extruded and deformed; thus, the permeability reduction is relatively slow. As the con ﬁ ning pressure increases, the equal-permeability point of fractured-porous sample moves to right. At the same con ﬁ ning pressure, the water saturation corresponding to equal-permeability point during depressurization is greater than that of pressurized process. It is also proved that the pore size decreases irreversibly, and the capillary force increases, which is equivalent to the enhancement of water wettability. Therefore, the irreversible closure of pores leads to the decrease of permeability and the increase of gas-phase seepage resistance, especially in carbonate rocks with fractures, vugs, and large pores. The ﬁ ndings of this study are helpful to better understand the gas production law of depletion development of carbonate gas reservoirs and provide support for e ﬃ cient development.

complicates the internal flow as shown in Figure 1 [6]. To develop carbonate reservoirs more effectively, it is necessary to study their internal flow rule [7].
Affected by bottom and edge water, potential water-gas two-phase flows exist in carbonate reservoirs due to the possible water invasion [8,9]. More than 95% of the developed carbonate gas reservoirs have edge-bottom water, of which more than 75% exhibit strong water influx [10]. More recently, a few researchers focused on the study of gas-water two-phase flow in carbonate using different methods. The fundamental regulation for fluid flow is the NS equation. The three most commonly used methodologies to investigate the fluid-fluid interactions in porous media by numerical simulation are bundle of capillary tube modeling (BCTM) [11], direct pore scale modeling (DPSM), and pore network modeling (PNM) [10]. DPSM generally requires reconstruction of three-dimensional digital rock to describe porous structure using micro-X-ray computed tomography (micro-CT) and SEM, then performs the flow simulation based on more accurate numerical simulation models than capillary tube. To complete the simulation of water-gas two-phase flow, the Lattice Boltzmann method (LBM) [12], which employs particle distribution function, the smoothed particle hydrodynamics (SPH) methods [13][14][15], the level set (LS) [16,17], volume of fluid (VOF) method [18], and the phase-field method (PFM) [19][20][21] are used to describe the interaction between fluids. Besides, PNM is incorporated to describe the complicated pore structure to simulate the internal flow by solving flow and transport equations on the network, which represents the pores and throats [22][23][24].
However, previous researches usually focused on characterization at one scale to analyze the pore structure of carbonate. For large-scale vugs and fractures, several models have been developed to simulate the production of carbonate reservoirs [25,26]. The studies about the gas-water two-phase flow of carbonate rocks mainly relied on experimental results to obtain phase permeability curves at the core scale [8,27]. Most of the numerical models used in previous studies were mainly based on ideal pore structure or thin sections, basically in two-dimensional, which cannot restore real pore structure. Using advanced mico-imaging methodologies such as micro-CT and SEM, we can obtain the detailed microstructures of carbonate including micro-fractures, micro-vugs, and interparticles [28,29]. This paper is aimed at studying the effects of stress on two-phase flow characteristics during the pressurized process and depressurization based on digital rocks obtained by X-ray computed tomography.
The carbonate gas reservoirs with deep burial depth are characterized by high temperature, high pressure, and high stress. The fractures and vugs pose significant challenges to achieving efficient recovery rate mainly due to their stresssensitive nature [30,31]. In previous work, we have characterized the pore structure variation under the change of the confining pressure in vuggy and fractured carbonate rocks [32]. The results indicated that the porosity showed an exponential decreasing trend and then followed an increase that did not recover to its initial value during pressurized process and depressurization process, due to the stress effect on large structures. The changes of pore structures will inevitably affect fluid flow. Thus, the influence of pressure on flow characteristics cannot be neglected.
In this paper, we extend our previous study and investigate the effect of stress sensitivity on the flow characteristics of carbonate reservoirs [33]. We use the micrometer-resolution X-ray tomography to obtain the digital rocks of porous, fractured-porous and fractured-vuggy carbonate rocks during pressurized process and depressurization. Based on digital rocks and pore network models, we simulate single-phase and two-phase flow using the Lattice Boltzmann method and pore network modelling, respectively. We then obtain   2 Geofluids the microscopic flow characteristics mainly the permeability and relative permeability of carbonate reservoirs considering different types of pore space and confining pressure. This work reveals the flow behaviors and their influencing factors in the three-dimensional porous media of carbonate, which provides a theoretical basis for the development of carbonate reservoirs.

Experimental and Simulation Methods
We put a sample into the holder made of PEEK material, fixed the holder on the sample table, and scanned the sample under the initial state. Next, we gradually increased the confining pressure to 5 MPa, 10 MPa, 15 MPa, and 20 MPa, then gradually reduced the confining pressure to 15 MPa, 10 MPa, 5 MPa, and 0 MPa, and scanned the cores under each stress state. Finally, 9 groups of pore structure images under different confining pressures were obtained. On the basis of CT images, we carried out the following LBM and pore network model simulation to study the single-phase and two-phase flow.

Lattice Boltzmann Mathematical Model. The Lattice
Boltzmann method (LBM) is derived based on the molecule dynamics and the statistical mechanics [34][35][36]. A discrete velocity model is generated through the mass conservation, momentum conservation, and the energy conservation from particle scale to microscale [37,38]. Then, the particle distribution function is generated via the discretion of velocity model and is calculated by statistical mechanics, so that the macroscale velocity and pressure are obtained.
In LBM, according to the statistical mechanics, Boltzmann's equation is expressed by Equation (1): If the external force is ignored, the particle velocity of v is discretized as e i , and the Equation (1) can be rewritten as Equation (2): The uniform LBM model can be applied to porous media at microscale. This property gives LBM a great advantage in the study of nonequilibrium dynamic models, especially, when fluid flow involves interface dynamics and complex boundary geometry. To the uniform LBM model, the equilibrium distribution function is written as Equation (3).

Geofluids
In the uniform model, u ′ and u are defined as Equation (4) and Equation (5): When only the resistance tensor on the main diagonal is considered, the above equations (Equation (4) and Equation (5)) can be reduced to Equation (6): Therefore, the equilibrium distribution function of the uniform LBM model can be rewritten as Equation (7):  Figure 6: The images of streamlines in the porous sample.
2.2. Two-Phase Pore Network Simulation. Given that the image-based direct simulation results in ultrahigh computational cost, quasi-static pore network model is used to simulate the gas-water two-phase flow [39,40]. In this work, the drainage process in water-wet rock is simulated, in which water is the defended phase and is displaced by gas. At the initial state, the network model is saturated with water, then, gas is injected from the inlet pores. Combined with the percolation theory, the pores and throats are invaded in the order of capillary pressure [41,42]. After the fluid pressure exceeds the threshold pressure of a throat, the injected phase enters the throat totally or partially filled with defended phase. The threshold capillary pressure of a throat depends on fluid phase configuration inside it and the simulated displacement process. For this gas-driving-water process, only the piston-like displacement occurs. The capillary pressure of a throat is determined by the Young-Laplace equation: where σ gw is the gas-water interfacial tension; r 1 and r 2 are the principal radius of interface curvature. If the pore shape and contact angle are given, Equation (8) can be used to obtain the capillary pressure. The derivations of capillary pressures of throats with cross-sectional shape of circle, triangle, and square are displayed in supplementary materials.
As injecting gas phase invades the water saturated in the network model, water saturation decreases and flow capacity of gas increases. The water/gas saturation and relative permeability during the displacement process are calculated and plotted as two-phase relative permeability curve. The detailed calculation processes are shown in supplementary materials.

Results and Discussion
3.1. Influence of Stress Sensitivity on Single-Phase Flow Characteristics. The Lattice Boltzmann method (LBM) has the advantages of high simulation precision and easy implementation, but it has a high requirement on the running memory of the computer. We consider the requirements between computer operation efficiency and model size, and we extract a subvolume (200 × 200 × 200 voxel) from digital rocks (852 μm × 852 μm × 852 μm) of fractured-porous sample, porous sample, and fractured-vuggy sample (Figure 2). We simulate a single-phase flow based on the extracted digital rock sample at the same position under different confining pressures and analyze the influence of stress changes on the permeability.
In this study, we use Palabos, an open source LBM solver, to calculate the permeability. The first step is to prepare the image for the application. We take the digital rock of fractured-porous sample as an example, since Palabos can only read image files in BMP black and white binary format. Avizo is used to save 3D pore phase digital rock in BMP image format, as shown in Figure 3(a). Then we transform the BMP file to a DAT file with MATLAB as shown in Figre  3(b). Finally, the generated DAT file is imported into Palabos to calculate the permeability. The simulation results are imported into Paraview to draw the flow diagram, and the simulated streamlines of fractured-porous sample is shown in Figure 4.
By analyzing the above images, we note that the fractures are the main flowing channels for fractured-porous carbonate sample. With confining pressure increasing, the streamlines begin to decrease, indicating that the fracture conductivity begins to decline, while the streamlines in pores begin to disappear. When the confining pressure increases to 20 MPa, the streamlines except for fracture basically disappear, indicating that many pores have been closed. After the pressure is recovered (0 MPa), the number of streamlines in the fracture increases again, but the area with streamline is significantly reduced, showing that the rising and falling process of confining pressure makes the fracture close irreversibly, compressing fluid flow space. The result of LBM-simulated permeability data is presented in Table 1.
As the pressure increases, the permeability of fracturedporous sample declines rapidly at the low-stress stage ( Figure 5). This indicates that at the low-stress stage, the compressive deformation first occurs in the fractures with good connectivities, which leads to a poor connectivity of the main flow channels, and thus, the permeability drops significantly at the initial stage. However, in the pressure recovery process, the overall permeability recovers to about a half of the initial state after the complete pressure relief. The fracture recovery is slow. Even after the confining pressure is completely unloaded, the fracture morphology cannot return to the initial state, and an irreversible closure occurs. As a result, the 6 Geofluids permeability cannot return to the initial level after the stress cycle is completed and even drops significantly. Figure 6 shows that the streamline changes for the porous sample. The observation and analysis of the images show that, for porous carbonate rocks, the streamlines in pores are significantly reduced as the confining pressure increases. When the confining pressure increases to 20 MPa, no obviously connected streamlines are found, indicating that pores are constantly squeezed and deformed, and some flow channels are blocked. After the confining pressure disappears, it can be observed that the streamlines in the pores basically recover compared with the initial state, that indicates the rise and fall of confining pressure makes some pores close, but it has less impact on the main flow channels. Table 2 shows the permeability results of the carbonate rocks obtained from simulations.
As shown in Figure 7, the permeability of porous carbonate decreases rapidly at the low-stress stage. However, the permeability recovery is not obvious at the initial stage of the pressure drop. Until the pressure is 0, the permeability recovers greatly, which is basically consistent with the phenomenon of fractured-porous carbonate rocks. In addition, the permeability recovery ability of porous carbonate is greater than that of fractured-porous carbonate, which indicates that the stress sensitivity of porous sample is weak. Figure 8 shows the streamline of fractured-vuggy sample at the initial state confining pressure of 10 MPa. At the initial state, the streamline is concentrated around the vugs and discontinuous, which also leads to poor flow capacity of the fractured-vuggy sample. We find that some streamlines in pores disappear under a confining pressure of 10 MPa, indicating that the flow capacity of fractured-vuggy sample is significantly affected by pores, and the closure of pores may lead to the isolation of individual vugs and affect the overall permeability. Table 3 shows the simulated permeability of fractured-vuggy sample.
The permeability of the fractured-vuggy carbonate rock is lower than the other two samples due to its poor connectivity. With the increase of confining pressure, the decreasing trend of permeability is fast initially, and then slows down ( Figure 9). It is considered that at the initial stage of stress application, the large pores or vugs as the main flow channels are first compressed and deformed, which leads to a rapid increase of fluid flow resistance and obvious decrease of permeability. However, in the later stress application, the mainly squeezed pores are isolated small ones, accordingly, the permeability decreases slowly in the later period.
The stress sensitivity of reservoir rocks can be evaluated by irreversible permeability damage rate, and the calculation formula is as follows: where D p presents irreversible permeability loss rate, K i is the initial permeability, and K ' is the permeability at the initial static stress point. Based on Equation (9) Figure 8: The images of streamlines in the fractured-vuggy sample. sensitivity of fractured-porous sample is moderately strong, while that of porous sample is weak.

The Effect of Stress Sensitivity on Two-Phase Flow
Characteristics. Due to the poor connectivity of porous sample and fractured-vuggy sample, we select fractured-porous carbonate rock for two-phase flow simulation. Note that the volume of digital rock is in 400 × 400 × 400 voxels, and its actual physical size is 1704 μm × 1704 μm × 1704 μm. Then, we extract a corresponding pore network model and obtain the relative permeability curves for each confining pressure ( Figure 10). The specific parameters are as follows: (1) Air-water interfacial tension: 70 mN/m (2) Range of water-phase contact angle: 30°-60°( 3) Viscosity: 1 cp in water phase and 0.01806 cp in gas phase (4) Density: 1000 kg/m 3 for water phase and 0.72 kg/m 3 for gas phase The distribution range of equal-permeability point mainly concentrates at 65%-87%. As the confining pressure increases, the equal permeability points continue to move to right. Because the decrease of pore size leads to the increase of capillary force (capillary force curves are shown in Figure 11), and the equal permeability points move to right. Furthermore, during depressurization, the pore structure cannot be completely recovered. Thus, the equal perme-ability points of the depressurization process are closer to right than those of the pressurized process under the same confining pressure.

Summary and Conclusions
The Lattice Boltzmann method and the pore network model are used to simulate single-phase and two-phase flow under different confining pressures in fractured-porous, porous, and fractured-vuggy carbonate rocks. The variation rules of gaswater relative permeabilities of carbonate rocks with different pore structures and different confining pressures are summarized, which has a guiding significance for gas reservoir exploitation. However, there is one limitation of the current work. Due to the limited number of experimental samples, only three samples are selected to represent fractured-porous, porous, and fractured-vuggy carbonate rocks, respectively. If more representative cores are used for each pore structure to carry out pressure cycle test, the gas-water two-phase percolation law obtained by flow simulation on this basis will be more universal and practical.
The simulation results show that: (1) In the process of rising confining pressure, the decline trend of the permeabilities of the three types of carbonate rocks is fast initially, and then slows down. It indicates that in the low-stress stage, the deformation compaction first occurs in the main flow channels such as fractures, vugs, or large pores. As a result, the resistance of fluid flow increases rapidly and the permeability decreases significantly. However, at the high-stress stage, many isolated small pores are mainly squeezed, so the decrease in permeability in the later period of pressurization is relatively slow.
(2) For fractured-porous carbonate rocks, fractures are main flow channels. According to the permeability calculation, the conductivities of fractures are 2-3 orders of magnitude greater than those of pores. In the process of depressurization, the recovery ability of the fractures is poor. Even if the confining pressure is completely released, some  (3) With the confining pressure increasing, the permeability loss of the porous sample is lower than that of the fractured-porous sample, and its permeability recovery ability is stronger than that of fractured-porous carbonate rock during depressurization. It indicates that the stress sensitivity of the fractured-porous sample is stronger.
(4) Pore connectivity of the fractured-vuggy sample is poor, and its simulated permeability is relatively low. Streamlines in some pores and vugs disappear by increasing the con-fining pressure. It indicates that the flowing capacity of the fractured-vuggy sample is greatly affected by the closure of pores. The closure of pores may isolate some vugs being isolated and affect the overall permeability.
(5) Under different confining pressures, the distribution range of the equal permeability points of the fracturedporous sample mainly concentrates at 65%-87%. As the confining pressure increases, the equal permeability point moves to the right. At the same confining pressure, the equal permeability point of the gas-water relative permeability curve is closer to the right during depressurization. It indicates that  9 Geofluids the pore size decreases as the confining pressure increases; thus, the capillary force increases, and the saturation of the wetting phase (water) increases. Nomenclature f i ðx, tÞ: Particle distribution function of node x at time t in the direction of i f eq i : Equilibrium particle distribution function f : Single particle distribution functionr r: Coordination v: Particle velocity α: Acceleration caused by external force Qð f Þ: Collision operator f eq i : Equilibrium distribution function R: Ideal gas constant u: Velocity ρ: Fluid density T: Temperature e i : Discrete velocity ω: Correlated weight coefficient R: Resistance tensor k: Permeability tensor t: Time c s : Lattice sound speed I : Unit matrix τ: Relaxation time p: Pressure σ gw : Gas-water interfacial tension, N/m r 1 and r 2 : Principal radius of interface curvature, m θ gwr : Receding gas-water contact angle, rad r: The radius of pores and throats, m F d : Dimensionless correction factor R: The radius of curvature, m β i : The half angle of corner i, rad G: Shape factor of pores and throats, dimensionless A: The cross-sectional area of pores and throats, m 2 P: The perimeter of cross section of pores and throats, m A g : The effective area occupied by gas, m 2 A w : The effective area occupied by water, m 2 S w : Water saturation V iw : The volume occupied by water in pore (or throat) i, m 3 V i : The volume of pore (or throat) i, m 3 k ri : The relative permeability of phase i q tmi : Total flow rate of phase i for multiphase flow, m 3 /s q tsi : Total flow rate when the network model is saturated by single phase i, m 3 /s q i,jk : The flow rate of phase i from pore j to pore k, m 3 /s g i,jk : The hydraulic conductance of phase i between pore j and pore k, m 4 /(Pa•s) p i,j : Pressure of phase i in pore j, Pa p i,k : Pressure of phase i in pore k, Pa L jk : The distance between center of pore j and pore k, m L j : The length of half-pore j, m L k : The length of half-pore k, m L t : The length of the throat connecting pore j and pore k, m g i,t : Hydraulic conductance of the throat connecting pore j and pore k, m 4 /(Pa•s) g i,j : Hydraulic conductance of half-pore j, m 4 /(Pa•s) g i,k : Hydraulic conductance of half-pore k, m 4 /(Pa•s).

Conflicts of Interest
The authors declare that they have no conflicts of interest.