An Improved Relative Permeability Model for Gas-Water Displacement in Fractal Porous Media

Abstract

Many researchers have revealed that relative permeability depends on the gas-water-rock interactions and ultimately affects the fluid flow regime. However, the way that relative permeability changes with fractal porous media has been unclear so far. In this paper, an improved gas-water relative permeability model was proposed to investigate the mechanism of gas-water displacement in fractal porous media. First, this model took the complexity of pore structure, geometric correction factor, water film, and the real gas effect into account. Then, this model was compared with two classical models and verified against available experimental data. Finally, the effects of structural parameters (pore-size distribution fractal dimension and tortuosity fractal dimension) on gas-water relative permeability were investigated. It was found that the sticking water film on the surface of fracture has a negative effect on water relative permeability. The increase of geometric correction factor and the ignorance of real gas effect cause a decrease of gas relative permeability.

1. Introduction

Unconventional natural gas, especially shale gas, has become increasingly important in the global energy supply in the past decade [1,2]. Hydraulic fracturing technology has been successfully applied to the exploitation of unconventional natural gas for the permeability enhancement of shale gas reservoirs. The fracturing fluid retains in the target formation and changes the flow mechanism in a longer production time [3,4]. In this flow regime, the relative permeability, as an essential parameter for characterizing the behavior of two-phase flow in a complex porous medium, plays a critical role in shale gas production and CO₂ geological storage [5,6]. Previous studies focused on the flow characteristics of gas in porous media and ignored the significant effect of water in complex micropores on gas effective permeability. This behavior is harmful to understanding the mechanism of gas-water displacement in complex pore structures. For instance, it may lead to an inaccurate prediction of fracturing fluid flowback in shale gas production. Therefore, it is meaningful to study the gas-water relative permeability in complex fractal porous media for gas-water displacement.

The gas-water relative permeability has been investigated through experiments, theoretical modeling, and numerical simulations. Figure 1 summarizes some classical relative permeability models and provides a brief introduction to the development of relative permeability model in recent years. The Purcell model calculated the permeability through capillary pressure data and formulated a relatively primitive permeability model [7]. This model assumes that the sum of the gas and water relative permeability is equal to 1, which does not agree with the actual fluid flow in porous media. Subsequently, this model was further modified by many scholars. Gates and Leitz [8] integrated the relative permeability from the capillary pressure. Burdine [9] thought that the assumption of capillary parallelism was not accurate and introduced the tortuosity factor into the Purcell model. Mualem [10] proposed an integration method based on the Purcell model. For simplicity, Romm [11] developed a relative permeability model which ignored the phase interference between the gas and water phases. This model is also called the X model. Brooks and Corey [12] introduced a pore-size distribution index to modify the capillary pressure function and proposed a more generalized gas-water relative permeability model. The Brooks–Corey model has been widely used in modeling two-phase flow in complex porous media. Fourar and Lenormand [13] derived a viscous coupling model after integrating Stokes’ equation and the effect of viscosity. Combined with momentum balance, viscosity, and cubic law, Chima and Geiger [14] developed a new relative permeability model. After ignoring the capillary pressure, Chima’s model became a unique form of the viscous coupling model. Li et al. [15] improved Chima’s model through the consideration of the influences of capillary tortuosity and irreducible water saturation on relative permeability. The above relative permeability models have made significant contributions to the prediction of two-phase flow in porous media. They also found that the relative permeability of the gas phase is more susceptible to the geometry of flow channels. However, these relative permeability models did not thoroughly consider the effects of complex pore microstructures on fluid flow.

Complex nanoscale pore structures of shale gas reservoirs can be directly observed by the advanced image technology such as focused ion beam-scanning electron microscopes (FIB-SEM) and broad ion beam-scanning electron microscopes (BIB-SEM), but the inclusion of these complex microstructures into the relative permeability model has become a challenging problem. There are several mathematical and physical issues to be solved. The first issue is the description of complex microstructures. Mandelbrot [16] first developed the fractal theory for self-similarity in complex porous media. The fractal theory has been regarded as a useful tool to describe nanoscale pore structures [17,18,19]. The second issue is the interaction between complex microstructures and fluid flow in porous media. Many scholars have discussed the influence of the complex pore structure on the flow mechanism in the relative permeability model [20,21,22]. For example, Li et al. [23] considered some non-interconnected tortuous capillaries in fractal-like porous media and proposed a novel gas-water displacing model. Lei et al. [17] developed a relative permeability model in multiscale porous media, which takes lognormal distribution function and residual water saturation into consideration. Ghanbarian et al. [21] developed a relative permeability model to consider the tortuosity-connectivity factor and effective saturation. The influences of drawdown pressure on the relative permeability of water and gas phase were experimentally investigated by Mo et al. [24]. Their experimental results indicated that gas relative permeability increases with drawdown pressure. Additionally, the increase in drawdown pressure can reduce irreducible water saturation. Huo and Benson [25] experimentally investigated the dependence of relative permeability on stress. They found that both irreducible water saturation and endpoint gas phase relative permeability decrease with the increase of effective stress at the same flow rate. Li et al. [1] stated that the interface effect plays a significant role in gas-water relative permeability. When the contact angle is higher than 90 °, the surface of hydrophobic material decreases the resistance of fluid flow, and the relative permeability correspondingly increases with the contact angle. Lei et al. [26] proposed a fractal model for the prediction of gas-water relative permeability in fractures. They found that increasing the complexity of pore structure (fractal dimension) would result in a decrease of water relative permeability, while the gas relative permeability has a corresponding increase with the ratio of fluid viscosity. Geng et al. [18] deduced a fractal apparent permeability model for shale gas flow in the heterogeneous matrix (organic and inorganic matters). If the organic matter in the matrix is ignored, the apparent permeability will be overestimated up to 24.1%. Therefore, pore shape, maximum pore size, and porosity are primary factors that have significant effects on apparent permeability [27]. Cai et al. [28] proposed an improved permeability model to consider viscous flow, molecular diffusion, and surface diffusion. However, this model uses a uniform diffusion coefficient to deduce gas permeability, ignoring the effect of water on the transport mechanism. Wang et al. [29] considered surface effects (hydrogen bond, double layer repulsive force, Vander Waals force, and structure repulsive force) into the gas-water relative permeability model. They found that the effect of hydrogen bond can significantly enhance fluid flow in porous media. However, the interaction of microstructures and fluid flow regimes has not been well-described in relative permeability models.

Because of the extremely low permeability of shale and complex pore structure, the relative permeability is difficult to measure in experiments. Thus, an analytical method for deriving gas-water relative permeability is critical. In the above-mentioned relative permeability models, the complexity of the pore structure, geometric correction factor, water film, and real gas effect are not comprehensively considered. It is still not clear how the irreducible water in the shale gas reservoir affects the gas flow in the complex pore structure. In this paper, an improved gas-water relative permeability model was proposed based on fractal theory. This improved fractal model takes the complexity of pore structure, geometric correction factor, water film, and real gas effect into consideration. Then, the improved fractal model was compared with two classical models and verified against available experimental data. Finally, the effects of structural parameters, geometric correction factor, water film, and real gas effect on gas-water relative permeability were investigated.

2. Proposed Gas-Water Relative Permeability Model

In this section, the two-phase flow mechanism in complex pore structures is investigated. The threshold pressure gradient and the critical capillary radius are depicted in a conceptual model to emphasize the influence of capillary size on the gas-water displacing process. Then, the fractal dimensions are introduced to describe the complexity of porous media. Finally, an improved gas-water relative permeability model is derived from accounting for the influences of the complex pore structure, geometric correction factor, water film and the real gas effect.

2.1. Two Classical Relative Permeability Models

2.1.1. Brooks–Corey Model

In the Corey model, relative permeability is a function of capillary pressure as

$$K_{rw}={\left(\frac{h}{p_c}\right)}^m$$

(1)

where $h$ is a constant, $p_c$ is the capillary pressure in porous media, $m$ is a correction coefficient, and $K_{rw}$ is the water relative permeability.

For complex pore structures, this equation cannot predict the relative permeability with sufficient accuracy. Researchers have found that relative permeability is sensitive to pore-size distribution. The classical Brooks–Corey model introduces a macroscopic parameter of pore-size distribution index $\lambda$ to describe the heterogeneity of pore structure. Their relative permeability model is

$$K_{rw}={\left(S_w^{*}\right)}^{\left(2+3\lambda \right)/\lambda }$$

(2)

$$K_{rg}={\left(1-S_w^{*}\right)}^2\left[1-{\left(S_w^{*}\right)}^{\left(2+\lambda \right)/\lambda }\right]$$

(3)

where $S_w^{*}$ represents the effective water saturation and is a function of capillary pressure $p_c$.

2.1.2. Viscous Coupling Model

Considering the effect of viscosity on fluid flow, Fourar and Lenormand [13] developed a gas-water relative permeability model as

$$K_{rw}=\frac{S_w^2}{2}(3-S_w)$$

(4)

$$K_{rg}=(1-S_w)+\frac{3}{2}{\mu }_r{S}_w(1-S_w)(2-S_w)$$

(5)

where $S_w$ is the water saturation. ${\mu }_r$ is the viscosity ratio of gas to water and is ${\mu }_r={\mu }_g/{\mu }_w$.

2.2. Gas-Water Displacement and the Variation of Capillary Pressure with Water Saturation

Figure 2 presents a simplified conceptual diagram to express the variation of capillary pressure with water saturation. The solid red line represents the process of gas displacing water. The initial state of the rock formation is fully water-saturated, and the pressure increases with gas injection. When the capillary entry pressure is overcome, the front of two-phase flow starts to move. The accumulated pressure will cause the gas to break through the upper boundary of the rock formation. This pressure is called breakthrough pressure. Finally, the increase of pressure does not further reduce water saturation in the fracture. This minimum water saturation is irreducible water saturation (IRWS). The blue dotted line represents the process of water imbibition. As the capillary pressure decreases, a portion of the gas is trapped in tinny pores or fractures under the wetting phase pressure. This minimum capillary pressure is often termed as snap-off pressure. The corresponding saturation is irreducible gas saturation (IRGS).

Both the pore-size distribution index and the capillary entry pressure are macroscopic parameters. They can express the complexity of rock pore structure. Such a model can macroscopically match the experimental data of permeability but does not describe the interaction of gas and water from the microscopic mechanism. In fact, the effects of pore-size distribution and tortuosity on permeability should be discussed from the viewpoint of fractal dimension. From a capillary theoretical model to the integration of capillaries in different radius sizes, fractal theory can better explore the effects of complex pore structure on permeability. Besides, the capillary entry pressure is usually regarded as a fixed value in the classical Brooks–Corey model. However, this value is difficult to determine using experimental tests. In this paper, we defined the critical capillary radius, which corresponds to the capillary entry pressure. When the capillary tube size is smaller than the critical radius of the capillary, the gas cannot displace the water, and the water flow fills the entire capillary flow channel. On this sense, this part of water is the irreducible water and its saturation is called the irreducible water saturation in the two-phase flow.

2.3. Irreducible Water Saturation and Threshold Pressure Gradient

This concept of irreducible water saturation is very similar to the threshold pressure gradient (TPG) that was mentioned in our previous study [6]. From the macroscopic view, the residual water in the target reservoir blocks the pore-throat. Only when the capillary pressure is higher than the capillary entry pressure, the two-phase flow can be observed, as shown in Figure 3a. TPG is a vital parameter in describing the gas starting flow against viscous forces. This TPG makes the Darcy flow gradually transform to a low-velocity, non-Darcy flow when water saturation is gradually increasing. Many researchers have indicated that permeability, water saturation, and pore pressure are controlling parameters for TPG in tight gas reservoirs [23,25,30]. On this sense, we proposed the following formula for TPG:

$$\begin{array}{l}TPG=a{e}^{s_w}\cdot e^{b\ast p_f+c\ast p_f^2}\\ a=e^{-(5.4491\times {10}^{15}\times k+1.13293)}\\ b=8.81931\times {10}^{13}\times k-0.06892\\ c=-5.83713\times {10}^{11}\times k+0.0129\end{array}$$

(6)

where $k$ is the intrinsic permeability, $p_f$ is the pore pressure, $a$ is a correction coefficient for water saturation, and $b$ and $c$ are correction coefficients for pore pressure. These correction coefficients are also linked with shale permeability.

The accuracy of Equation (6) was checked using the experimental data from Ding et al. [30]. Their comparison is shown in Figure 4. Figure 4a presents the rise of TPG with the increase of water saturation. Figure 4b shows the effect of different permeabilities on TPG. Although TPG describes the gas-water displacement in capillary tubes, it does not consider the complex pore structure of porous media.

2.4. Calculation of Water Saturation and Irreducible Water Saturation

Based on the Laplace formula, capillary pressure can be expressed as

$$p_c=\frac{2\sigma \cos \theta }{r_c}$$

(7)

where $p_c$ is the capillary pressure, $\sigma$ represents the surface tension, and $\theta$ is the contact angle of gas-water-rock interface.

If the surface tension and contact angle are constant, the capillary pressure has an inverse proportional relationship with the pore radius. There is a critical capillary radius $r_c$ in the complex distribution of capillary pore size. When the pore size is larger than the critical capillary radius, the gas flow can occur in the capillary tube, and the water film is attached to the capillary wall. If the capillary tube size is smaller than the critical radius $r_c$, the capillary tube is fully saturated with water and gas flow does not occur. Therefore, from the microscopic view, the gas can displace water only when the capillary tube radius is greater than the critical capillary tube radius, as described in Figure 3b.

The water sources come from two parts: One is a single water flow $\left(r<r_c\right)$, and the other is the gas-water two-phase flow $\left(r>r_c\right)$. According to the area ratio of pores, the water saturation can be calculated by pore-size distribution and water film as [20]

$$S_w=\frac{{\displaystyle {\int }_{r_{\min }}^{r_c}\pi r^2(-dN)+{\displaystyle {\int }_{r_c}^{r_{\max }}\pi \left(r^2-r_{wg}^2\right)(-dN)}}}{{\displaystyle {\int }_{r_{\min }}^{r_{\max }}\pi r^2(-dN)}}$$

(8)

When the pore size is less than the critical capillary tube radius, there is not enough capillary pressure to further displace the water in pores. At this time, the water saturation is the irreducible water saturation and is calculated by

$$S_{wi}=\frac{{\displaystyle {\int }_{r_{\min }}^{r_c}\pi r^2(-dN)}}{{\displaystyle {\int }_{r_{\min }}^{r_{\max }}\pi r^2(-dN)}}=\frac{r_c^{2-D_f}-r_{\min }^{2-D_f}}{r_{\max }^{2-D_f}-r_{\min }^{2-D_f}}$$

(9)

Figure 5 shows the variation of irreducible water saturation with critical capillary tube radius. No linear relationship was observed between the critical capillary radius and irreducible water saturation. When the critical capillary radius is large enough, the irreducible water saturation grows slowly and tends to be stable. In the case of the same capillary critical tube radius, the increase in the maximum capillary radius leads to a decrease in irreducible water saturation. The water flow is enhanced by the increase of capillary tube radius rather than being sealed in the capillary tube. When the critical capillary tube radius is small relative to the capillary tube size, the difference is not obvious. However, the irreducible water saturation changes from 0.2 to 0.4 when the critical capillary tube radius is 10 nm.

2.5. Fractal Distribution of Pores

Porous media have numerous capillary tubes with different pore sizes and tortuous lengths. These inter-pore structures are of self-similarity and can be described by fractal theory. The number of pore sizes greater than $r$ can be expressed as

$$N\left(\delta \ge r\right)={\left(\frac{r_{\max }}{r}\right)}^{D_f}$$

(10)

where $\delta$ represents the scale of length, $r_{\max }$ and $r$ are the maximum pore radius and pore radius, respectively, and $D_f$ is the pore-size distribution (PSD) fractal dimension and in the range of $0<D_f<2$ in two dimensions. Yu and Li [31] found that this fractal dimension is related to the porosity $\varphi$ and the radius ratio $r_{\min }/r_{\max }$.

$$D_f=d-\frac{\ln \varphi }{\ln \left(r_{\min }/r_{\max }\right)}$$

(11)

where $d$ represents the Euclidean dimension.

Differentiating Equation (10) with respect to $r$ obtains the following expression:

$$-dN=D_f{r}_{\max }^{D_f}{r}^{-(D_f+1)}dr$$

(12)

Thus, the probability density function of pore-size distribution is expressed as

$$f(r)=D_f{r}_{\max }^{D_f}{r}^{-(D_f+1)}$$

(13)

Due to the tortuosity, the length of the capillary tube is higher than the straight capillary tube in complex porous media. Capillary length is related to the radius size with following fractal law [31]:

$$L(r)={(2r)}^{1-D_t}{L}_0^{D_t}$$

(14)

where $L$ and $L_0$ represent the actual capillary length and representative capillary length, respectively. $D_t$ is the fractal dimension of the tortuosity, and its range is between 1 and 2 in two-dimensional space. Therefore, the cross-sectional area of all pores could be written as [20]

$$A=\frac{\left(1-\varphi \right)}{\varphi }\frac{D_f\pi r_{\max }^2}{2-D_f}$$

(15)

Thus, the representative capillary length can be expressed as

$$L_0=r_{\max }^{}\sqrt{\frac{\left(1-\varphi \right)}{\varphi }\frac{D_f\pi }{2-D_f}}$$

(16)

2.6. Two-Phase Flow Model in Nanopores

2.6.1. Flow Rate of Gas

The flow channels in Figure 6 were assumed to be capillary tubes with tortuosity. SEM images showed that the pore shape was not only circular, but also slit. Ignoring the influence of pore shape will affect the understanding of the fluid flow mechanism in complex porous media. Mortensen et al. [32] stated that the capillary channel with a circular cross-section does not reflect the complexity of various cross-sectional shapes, such as elliptic, rectangular, triangular, and harmonic-perturbed circles. We introduced a geometrical correction factor to describe different cross-sectional shapes in a unified form.

When a steady-state laminar flow occurs within capillary tubes, Figure 7 gives the flow distribution of the gas and water phases. This conceptual model has three different regions: The gas flow region, the typical viscosity water flow region, and a high-viscosity water film region. The fluid flow in these three regions has different flow mechanisms. They are mathematically described below.

Based on the Hagen–Poiseuille equation, the gas flow rate in a capillary tube is obtained below if slip flow and free molecular flow are considered in the continuum flow regime.

$$q_g=\frac{\pi r_{}^4}{128{\mu }_g}(1+\beta K_n^{*})\left(1+\frac{4K_n^{*}}{1-\alpha K_n^{*}}\right)\frac{\Delta p}{L}$$

(17)

where $r$ and $L$ are the capillary tube radius and capillary tube length, respectively, $\Delta p$ is the pressure difference applied to the external capillary tube, ${\mu }_g$ is the gas viscosity, and $\beta$ represents a dimensionless rarefaction coefficient. When $\beta$ changes from 0 to $\frac{64}{3\pi \left(1-4/\alpha \right)}$, a slip flow regime transforms into a free molecular flow regime. $\alpha$ represents the slip coefficient and can be obtained from the linearized Boltzmann solution or velocity distribution given by the Poiseuille flow in a cylindrical tube [35]. $K_n^{*}$ is the modified Knudsen number, which is the ratio of the mean free path of gas flow to the characteristic size of gas flow channel. Details are provided in the next section.

2.6.2. Geometrical Correction Factor and Real Gas Effect

Shale formation has various shapes of pores. This study introduced a geometric correction factor to express the deviation of pore shape from a circular pore channel [36,37].

$$K_n^{*}=\frac{\eta }{2\gamma r}$$

(18)

where $\gamma$ is the geometric correction factor. The $\gamma$ is 1 for the circular pore, 1.094 for the square pore, and is 1.1186 for the equilateral triangle pore [18]. $\eta$ is the gas mean free path.

$$\eta =\frac{{\mu }_g}{p_g}\sqrt{\frac{\pi ZRT}{2M}}$$

(19)

where $p_g$ is the gas pressure, $R$ is the universal gas constant, $T$ is the reservoir temperature, $M$ is the gas molecule mass, and $Z$ is the gas compressible factor, which expresses the real gas effect on gas density and viscosity. It was noted that previous models usually neglected this real gas effect (that is Z = 1), but this effect becomes more notable at high pressure and temperature. $Z$ can be expressed as [38]:

$$\begin{array}{l}Z=1+A_1{p}_r+A_2{p}_r^2+\frac{A_3{p}_r^{A_4}}{T_r^{A_5}}+\frac{A_6{p}_r^{\left(A_4+1\right)}}{T_r^{A_7}}+\frac{A_8{p}_r^{\left(A_4+2\right)}}{T_r^{\left(A_7+1\right)}}\\ p_r=p/p_{gc},T_r=T/T_{gc}\end{array}$$

(20)

where $p_{gc}$ and $T_{gc}$ are the gas critical pressure and temperature, respectively, and $A_1$–$A_8$ are tuned coefficients.

Jarrahian and Heidaryan [39] developed a generalized gas viscosity model with a range of pressures (0.1–137.8 MPa) and temperatures (241–473 K). This study described their expression of gas viscosity as

$${\mu }_g={\mu }_a\left[1+\frac{A_{11}}{T_r^5}\left(\frac{p_r^4}{T_r^{20}+P_{\mathrm{r}}^4}\right)+A_{12}{\left(\frac{P_r}{T_r}\right)}^2+A_{13}\left(\frac{P_r}{T_r}\right)\right]$$

(21)

where ${\mu }_a$ is the gas viscosity under the standard condition, and $A_{11}$, $A_{12}$, and $A_{13}$ are temperature, pressure, and molar weight independently tuned coefficients, respectively.

Beskok and Karniadakis [40] found that gas flow is in the slip flow regime when $\beta =0,\alpha =-1$. Combining Equations (17)–(21), we extended the gas flow rate to consider the geometric correction factor and real gas effect in a capillary tube and obtain the flow rate as

$$q_g=\frac{\pi r_{wg}^4}{8{\mu }_g}\left(1+\frac{4K_n^{*}}{1+K_n^{*}}\right)\frac{\Delta p}{L}$$

(22)

Integrating Equation (22), the total gas volumetric flow rate is expressed as

$$\begin{array}{ll}{Q}_g& =D_f{r}_{\max }^{D_f}{\displaystyle {\int }_{r_c}^{r_{\max }}\frac{\pi r_{wg}^4}{8{\mu }_g}}\left(1+\frac{4K_n^{*}}{1+K_n^{*}}\right)\frac{\Delta p}{L}\cdot r^{-(D_f+1)}dr\\ & ={\epsilon }_g{\displaystyle {\int }_{r_c}^{r_{\max }}{r}_{wg}^4\cdot r^{D_t-D_f-2}}(1+4K_n^{*})dr\end{array}$$

(23)

where $r_{wg}$ is the distance from the center of the circular pore to the gas-water interface and $r_c$ is the critical capillary tube radius. ${\epsilon }_g=D_f{r}_{\max }^{D_f}\frac{\pi \Delta p}{2^{4-D_t}{\mu }_g{L}_o^{D_t}}$.

2.6.3. Water Film Effect

The water with normal viscosity occupies the most of capillary channels. Due to influences of the hydrophilicity of the matrix, the high-viscosity water forms a water film on the surface of capillary tubes. For the flow region of normal viscosity water, the governing equation can be obtained based on the Hagen–Poiseuille equation.

$${\mu }_w\frac{1}{r}\frac{d}{dr}\left(r\frac{d{v}_w}{dr}\right)=\frac{\Delta p}{L}r\in \left[r_{wg},r_{wf}\right]$$

(24)

The governing equation of water film region is

$${\mu }_{wf}\frac{1}{r}\frac{d}{dr}\left(r\frac{d\left(v_{wf}+\zeta v_v\right)}{dr}\right)=\frac{\Delta p}{L}r\in \left[r_{wf},r_0\right]$$

(25)

where ${\mu }_w$ and ${\mu }_{wf}$ are the water viscosity for the typical viscosity water region and water film region, respectively, $\zeta$ is the correction factor of the vapor flow rate, and $v_v$ represents the vapor flow rate.

The ratio of these two viscosities is related to the contact angle $\theta$ [41]:

$$\frac{{\mu }_w}{{\mu }_{wf}}=-0.018\theta +3.25$$

(26)

With the combined effects of temperature, pressure, and gas flow rate, the water film adhering to the surface of the channel will diffuse into the gas under evaporation reaction and be extracted with the gas. The influence of water vapor on the pressure of the gas is small and was thus ignored. However, the evaporation of the water film directly affects the flow rate of the water film. If the effect of water vapor flow rate on the water film flow rate is not considered, the water film flow rate will be overestimated. This vapor flow rate could be calculated as [42]

$$v_v=\frac{\pi M_v{D}_{va}}{RT}\frac{\Delta p_v}{L}$$

(27)

where $M_v$ is the molecule mass of vapor, $D_{va}$ is the diffusion coefficient, and $\Delta p_v$ is the pressure difference of vapor.

The boundary conditions for water flow and water film regions are as follows:

$$\{\begin{array}{l}{v}_w|{}_{r=r_{wf}}=v_{wf}|{}_{r=r_{wf}}\\ v_w\frac{\partial v_w}{\partial r}|{}_{r=r_{wf}}=v_{wf}\frac{\partial v_{wf}}{\partial r}|{}_{r=r_{wf}}\\ v_{wf}|{}_{r=r_0}=0\end{array}$$

(28)

By substituting Equations (26)–(28) into Equations (24) and (25), the final flow velocities are obtained as

$$v_w=\frac{\left(r_{wf}^2-r^2\right)}{4{\mu }_w}\frac{\Delta p}{L}+\frac{\left(r_0^2-r_{wf}{}^2\right)}{4{\mu }_{wf}}\frac{\Delta p}{L}r\in \left[r_{wg},r_{wf}\right]$$

(29)

$$v_{wf}=\frac{\left(r_0^2-r^2\right)}{4{\mu }_{wf}}\frac{\Delta p}{L}r\in \left[r_{wf},r_0\right]$$

(30)

By integrating Equations (29) and (30) on the cross-section of the circular tube, the volumetric flow rate of the water are obtained as

$$\begin{array}{ll}{q}_w& ={\displaystyle {\int }_{r_{wg}}^{r_{wf}}{v}_w}dA+{\displaystyle {\int }_{r_{wf}}^{r_0}{v}_{wf}}dA\\ & =\frac{\pi }{8{\mu }_w}\frac{\Delta p}{L}{\left(r_{wf}^2-r_{wg}^2\right)}^2+\frac{\pi }{4{\mu }_{wf}}\frac{\Delta p}{L}\left(r_{wf}^2-r_{wg}^2\right)\left(r_0^2-r_{wf}^2\right)+\frac{\pi }{8{\mu }_{wf}}\frac{\Delta p}{L}{\left(r_0^2-r_{wf}^2\right)}^2\end{array}$$

(31)

When only water flows in the capillary tube ($r_{wg}=0$), the volumetric flow rate of water can be simplified as

$$q_{w-\sin gle}=\frac{\pi \Delta p}{8{\mu }_{w}L}{r}_{wf}^4+\frac{\pi \Delta p}{8{\mu }_{wf}L}\left(r_0^2-r_{wf}^2\right)\left(r_0^2+r_{wf}^2\right)$$

(32)

Similarly, combining with fractal theory obtains the volumetric flow rate of water as

$$\begin{array}{c}{Q}_w=D_f{r}_{\max }^{D_f}{\displaystyle {\int }_{r_{\min }}^{r_c}\left[\frac{\pi }{8{\mu }_w}\frac{\Delta p}{L}{\left(r_{wf}^2-r_{wg}^2\right)}^2+\frac{\pi }{4{\mu }_{wf}}\frac{\Delta p}{L}\left(r_{wf}^2-r_{wg}^2\right)\left(r_0^2-r_{wf}^2\right)\right]}\\ +D_f{r}_{\max }^{D_f}{\displaystyle {\int }_{r_c}^{r_{\max }}\left[\frac{\pi }{8{\mu }_w}\frac{\Delta p}{L}{\left(r_{wf}^2-r_{wg}^2\right)}^2+\frac{\pi }{4{\mu }_{wf}}\frac{\Delta p}{L}\left(r_{wf}^2-r_{wg}^2\right)\left(r_0^2-r_{wf}^2\right)+\frac{\pi }{8{\mu }_{wf}}\frac{\Delta p}{L}{\left(r_0^2-r_{wf}^2\right)}^2\right]}\\ ={\epsilon }_w\left[{\displaystyle {\int }_{r_{\min }}^{r_{\max }}\left(r_{wf}^2-r_{wg}^2\right)}\left(r_{}^2-r_{wg}^2\right)\cdot r^{D_t-D_f-2}dr+{\displaystyle {\int }_{r_c}^{r_{\max }}{\left(r_{}^2-r_{wg}^2\right)}^2}\cdot r^{D_t-D_f-2}dr\right]\end{array}$$

(33)

where ${\epsilon }_w=D_f{r}_{\max }^{D_f}\frac{\pi \Delta p}{2^{4-D_t}{\mu }_w{L}_0^{D_t}}$.

2.7. Intrinsic Permeability

Intrinsic permeability is an inherent property of rock. It is only related to the rock microstructures and is thus independent of fluid type. For two-phase flow, the effective permeability of each phase depends on its saturation and intrinsic permeability. Experimental results have suggested that the effective permeability of non-wetting phase tends to be higher than the effective permeability of wetting phase [43]. Relevant literature has summarized following reasons [44]: (1) The rock is rich in swelling clay minerals. When the effective permeability of the wetting phase is measured, the clay mineral in the matrix absorbs water and swells. This leads to the reduction of the flow channel of fractures; (2) The existence of a water film on the structure surface reduces the effective radius of the flow channel; (3) The high viscosity of the movable water film zone increases the viscous force of normal flowing water. Thus, higher shear stress is required to start the water flow. Figure 8 compares several sets of experimental data on the relationship between water and gas permeability. These experimental data were obtained from Busch and Hildenbrand [43]. These data showed that the gas effective permeability is slightly higher than the water effective permeability. This is related to the complex micropore structure and the flow mechanism of two-phase flow.

Two mainstream methods are usually used to calculate the intrinsic permeability of rock. The first is the Kozeny–Carman equation or its modified form. The intrinsic permeability of rocks, especially with high permeability, is expressed as

$$k=\frac{B{r}_{mean}^2}{K}\frac{{\varphi }^3}{1-{\varphi }^2}$$

(34)

where $B$ is a constant, $r_{mean}$ is the average radius of grains. $K$ represents the Kozeny–Carman constant and has

$$K=c_K\tau$$

(35)

where $c_K$ is the dimensionless Kozeny constant and $\tau$ is the tortuosity of rock.

Another method was proposed by Yang and Aplin. This method does not only consider the pore size and porosity, but also includes the effect of mineral composition (clay content) in the intrinsic permeability:

$$k=2.2446\times {10}^{-17}{J}_v^{1.18}{r}_{mean}^{1.074}$$

(36)

where

$$\begin{array}{l}{J}_v=\frac{9}{8}\varphi {\left(\sin o\right)}^2\frac{J_1^3}{{\left(1+J_1+J_1^2\right)}^2}\\ J_1=2.371-1.626{\chi }^2+153.8{\varphi }^4\\ o={45}^{°}-{10.24}^{°}\left(e_{100}-e\right)\\ e_{100}=0.3024+1.687\chi +1.951{\chi }^2\\ e=\varphi /\left(1-\varphi \right)\end{array}$$

(37)

$e$ is the volume ratio of pores to matrix and $e_{100}$ is the volume ratio of pores to matrix at the effective stress of 0.1 MPa, $\chi$ is the clay mineral content, and $o$ represents the angle between the average pore alignment angle and bedding direction.

2.8. Gas-Water Relative Permeability Model

The water flow velocity is expressed in Darcy’s law:

$$Q_w=\frac{k_w\varphi }{{\mu }_w{S}_w}\frac{\Delta p}{L}$$

(38)

By substituting Equation (33) into Equation (38), the effective permeability of water is obtained as

$$k_w=-\frac{{\mu }_w{S}_w{L}_0{\epsilon }_w}{A\Delta p}\left[{\displaystyle {\int }_{r_{\min }}^{r_{\max }}\left(r_{wf}^2-r_{wg}^2\right)}\left(r_{}^2-r_{wg}^2\right)\cdot r^{D_t-D_f-2}dr+{\displaystyle {\int }_{r_c}^{r_{\max }}{\left(r_{}^2-r_{wg}^2\right)}^2}\cdot r^{D_t-D_f-2}dr\right]$$

(39)

Similarly, the effective permeability of gas is

$$k_g=-\frac{{\mu }_g\left(1-S_w\right)L_0}{A\Delta p}{\epsilon }_g{\displaystyle {\int }_{r_c}^{r_{\max }}{r}_{wg}^4\cdot r^{D_t-D_f-2}}(1+4K_n^{*})dr$$

(40)

Thus, the relative permeability for water is

$$k_{rw}=\frac{-\frac{{\mu }_w{S}_w{L}_0{\epsilon }_w}{A\Delta p}\left[{\displaystyle {\int }_{r_{\min }}^{r_{\max }}\left(r_{wf}^2-r_{wg}^2\right)}\left(r_{}^2-r_{wg}^2\right)\cdot r^{D_t-D_f-2}dr+{\displaystyle {\int }_{r_c}^{r_{\max }}{\left(r_{}^2-r_{wg}^2\right)}^2}\cdot r^{D_t-D_f-2}dr\right]}{k}$$

(41)

and the relative permeability for gas is

$$k_{rg}=\frac{-\frac{{\mu }_g\left(1-S_w\right)L_0}{A\Delta p}{\epsilon }_g{\displaystyle {\int }_{r_c}^{r_{\max }}{r}_{wg}^4\cdot r^{D_t-D_f-2}}(1+4K_n^{*})dr}{k}$$

(42)

In this study, our intrinsic permeability was calculated by Yang and Aplin’s method. The clay content was 0.1 and the porosity of rock was 0.2.

3. Model Validation

The above-improved fractal relative permeability model for gas and water was compared with the classic Brooks–Corey model and viscous coupling model. Our fractal model was further verified against experimental data. The parameters used in this validation are listed in

Table 1. Parameters for model validation.

Parameters	Unit	Value	Physical Meanings	Reference
$r_{\min }$	m	${10}^{-9}$	Maximum pore radius	Zhang et al. [49]
$r_{\max }$	m	${10}^{-7}$	Minimum pore radius	Zhang et al. [49]
$R$	J/(mol·K)	8.314	Universal gas constant	Given
$M$	Kg/mol	0.016	Methane molecular weight	Given
$T$	K	298.15	Temperature	Given
${\mu }_{wf}$	Pa·s	$2.6811\times {10}^{-3}$	Water film viscosity	Wu et al. [41]
${\mu }_w$	Pa·s	$8.937\times {10}^{-4}$	Water viscosity	Wu et al. [41]
$\gamma$		1	Geometric correction factor	Geng et al. [18]
$p_g$	MPa	1	Gas pressure	Given
$\varphi$		0.2	Porosity	Zhang et al. [49]
$A_1$		0.007698	Tuned coefficient	Sanjari and Lay [38]
$A_2$		0.003839	Tuned coefficient	Sanjari and Lay [38]
$A_3$		−0.467212	Tuned coefficient	Sanjari and Lay [38]
$A_4$		1.018801	Tuned coefficient	Sanjari and Lay [38]
$A_5$		3.805723	Tuned coefficient	Sanjari and Lay [38]
$A_6$		−0.087361	Tuned coefficient	Sanjari and Lay [38]
$A_7$		7.138305	Tuned coefficient	Sanjari and Lay [38]
$A_8$		0.083440	Tuned coefficient	Sanjari and Lay [38]
$A_{11}$		7.9	Temperature tuned coefficient	Jarrahian and Heidaryan [39]
$A_{12}$		$9\times {10}^{-6}$	Pressure tuned coefficient	Jarrahian and Heidaryan [39]
$A_{13}$		0.28	Molar weight tuned coefficient	Jarrahian and Heidaryan [39]

.

The comparison of our fractal model with the classic Brooks–Corey model and the viscous coupling model is presented in Figure 9. The residual water saturation of the Brooks–Corey model was 0.2, and the pore-size distribution index was infinite.

The viscous coupling model takes the viscosity ratio of the gas and water phase into account. When the viscosity difference between the gas and liquid phase is significant, the viscous coupling model is similar to the Brooks–Corey model. Our fractal model had similar trends to the two classic models. However, our fractal model had some features that the classic models missed: (1) The pore-size distribution index in the classical model is constant and difficult to measure. However, pore-size distribution fractal dimension and tortuosity fractal dimension can be calculated by the capillary pore distribution and pore size. Through fractal dimensions, complex porous media can be better expressed; and (2) The real gas effect and water film flow are not considered in the classical models. These factors are particularly important in high-temperature, pressure, and humidity formations. Further, our fractal model was compared with experimental data, and the results are presented in Figure 10. These experimental data were measured by Mahiya [45] for the relative permeability of water and vapor at 120 °C. The permeability and porosity of Berea core samples used in this experiment were 1200 md and 0.24. Subsequently, the relative permeability of nitrogen and water were measured by Li and Horne [46]. Diomampo et al. [47] conducted drainage experiments for the rough-walled fractures, and relative permeability data of nitrogen and water were obtained. In their experiments, the fracture was initially saturated with water, and the non-wetting phase was injected with a constant flow rate. Wang et al. [48] measured the relative permeability of oil and water in carbonate core samples. The permeability and porosity of this sample were 33.9 md and 0.18. Our fractal relative permeability model is in a good agreement with these experimental data. This indicates that this improved fractal model can describe the evolution of gas-water relative permeability in complex porous media.

4. Discussions

4.1. Effects of Pore-Size Distribution and Tortuosity Fractal Dimension

The structural parameters, pore-size distribution (PSD) fractal dimension, and tortuosity fractal dimension can be used to describe the structural complexity of porous media. Porous media have multiscale capillary structures from micron to nanoscale. The distribution and heterogeneity of the pore structure can be described by the PSD fractal dimension. Figure 11 shows that the effective permeability of gas and water decreased significantly as the PSD fractal dimension increased from 1.2 to 1.8. When the tortuosity fractal dimension was constant at 1.48, the increase of the PSD fractal dimension increased the complexity of the pore-size distribution and resulted in a decrease of effective permeability for both gas and water.

The tortuosity fractal dimension was first introduced into this improved fractal model to reflect the variation of capillary tortuosity with pore size. Figure 12 provides the variation of effective permeability of gas and water with the increase of PSD fractal dimension. This trend was similar to the variation of the PSD fractal dimension, but the degree of decline was not significant. This indicates that the influence of the PSD fractal dimension on both gas and water permeability can be greater than the tortuosity fractal dimension.

4.2. Effect of Nanotube Size

The size of nanotubes is directly related to the gas-water effective permeability. When the maximum radius of the nanotubes was changed from 100 nm to 50 nm, the gas-water effective permeability significantly declined, as shown in Figure 13. Because the capillary tube radius was still very large relative to the size of critical capillary tube radius and water film thickness, residual water cannot block the capillary tube and single-phase flow did not occur. The effective permeability of the gas and water decreased at almost the same speed.

4.3. Effect of Water Film Thickness

Our fractal model has some parameters in the gas or water effective permeability. For example, the water film adhered to the surface of the capillary tube is an important parameter that can affect the effective permeability of water. The viscosity of the water film is greater than that of water in the normal flowing zone. An increase in the water film thickness directly reduces the water flow flux. Moreover, the water film produces a drag force on the normally flowing water. This can reduce the flow rate of the wetting phase. Some studies have considered the water film as the non-flowable film, which absorbs on the surface of the matrix to form residual water in the porous medium. This assumption exaggerates the effect of water film on the relative permeability of the wetting phase. Equations (29) and (30) provide the expressions of water flow rate and water film flow rate, respectively. If a capillary radius is 100 nm, water saturation is 0.5, and a pressure difference is 1 MPa, the effect of water film thickness on the water and water film flow rates are presented in Figure 14, where the evolution of water flow rate is plotted under different water film thicknesses. When the water film thickness increases from 0.7 nm to 7 nm, a significant decrease in water flow rate is observed. In addition, a slight decrease of the normal water flow zone is observed due to the increase of water film thickness.

Figure 15 shows the effect of water film thickness on the water film flow rate. It is worth noting that the flow rate of the water film is not affected by the change in water film thickness. There are two reasons: (1) The thickness of the water film is much smaller than the radius of the capillary tube; and (2) The water film is a high-viscosity water flow zone, which mainly causes dragging force on the normally flowing water. However, the increase of the water film thickness significantly increases the zone where the water film flows, and the flow flux of the water film obtained by this integration increases correspondingly. Finally, Figure 16 is the evolution of gas-water effective permeability with different water film thickness. The increase of water film thickness has little influence on the effective permeability of gas, but has a negative influence on the effective permeability of the wetting phase.

4.4. Effect of Geometric Correction Factor

The pores in porous media have various shapes. In this paper, the relative permeability was derived based on the circular tube, and the geometric correction factor was introduced to redefine the Knudsen number (see Equation (18)). The Knudsen number reflects the motion state of gas molecules. Therefore, the change of the Knudsen number directly affects the gas flow rate, resulting in a change of the gas effective permeability. Figure 17 provides the variation of Knudsen number with a geometric correction factor. It states that the shape irregularity of pores reduced the Knudsen number. As the capillary radius increases, this effect will gradually decrease. When the flow channel is large enough, the free path of gas molecules is much smaller than the characteristic length, and the shape irregularity of the pores is no longer significant.

Figure 18 shows the variation of gas-water effective permeability with the increase of geometric correction factor. The increase of the geometric correction factor directly reflects the increase of capillary irregularity. Therefore, the effective gas permeability decreases when the geometric correction factor increases. The effective permeability of the wetting phase is not affected by the shape correction factor. Furthermore, the effect of the geometric correction factor is less than the effect of structural parameters (PSD fractal dimension and tortuosity fractal dimension). This may be why the capillary is generally assumed to be a circular tube in many numerical simulations, regardless of its shape variation.

4.5. Real Gas Effect

The real gas effect refers to the variation of gas compressibility and viscosity relative to ideal gas. Both temperature and pressure affect the real gas effect and ultimately lead to a change of gas flow. Figure 19 presents the importance of real gas effect on the effective permeability of water and gas. The effective permeability of gas is slightly higher after considering the real gas effect. Although the viscosity of the gas increases, the gas is more compressible when the real gas effect is considered. In addition, the increase of the Knudsen number means an enhancement of gas slip effect, which is beneficial to the effective gas permeability. These indicate that the real gas effect cannot be ignored in the evolution of gas effective permeability. This effect will become more significant in high-pressure and temperature conditions.

5. Conclusions

In this paper, an improved fractal model was developed for the gas-water relative permeability in porous media to incorporate the structural parameters (PSD fractal dimension and tortuosity fractal dimension), water film, geometric correction factor, and real gas effect. First, the threshold pressure gradient was introduced to describe the gas starting flow against viscous forces in two-phase flow. This threshold pressure gradient was expressed by some macroscopic parameters, such as water saturation, pore pressure, and permeability. Their linkage with complex pore structures (critical capillary radius) was explored through a microscopic theoretical model. Then, the gas-water relative permeability model was derived based on the circular cross-section flow channels and complex fractal pore structures. The improved fractal model was compared with two classical models and verified against available experimental data. Finally, the sensitivities of structural parameters, water film, geometric correction factor, and real gas effect to gas-water effective permeability were investigated. Based on these investigations, the following conclusions can be drawn:

• Pore-size distribution fractal dimension and tortuosity fractal dimension can well-describe the complexity of pore-size distribution and tortuosity distribution of porous media, respectively. Their higher dimensions cause a decrease of the water-gas effective permeability. It was found that the pore-size distribution determines the flow pattern: Single-phase flow or two-phase flow and the pore-size distribution fractal dimension have a more significant influence on gas-water effective permeability.
• The water film thickness has a negative effect on the water effective permeability. It increases the water film zone and loses the space for normal water flows, thus causing a reduction in water flow rate. Moreover, the high-viscosity water film also applies a viscous force to normal water, reducing its flow rate. If the water film thickness is much smaller than the pore size, the effect of water film thickness on the gas effective permeability is negligible.
• The pore geometry directly affects the mechanism of gas flow. When the irregularity of the pore geometry increases, the Knudsen number decreases. The collision between gas molecules is intensified and the gas flow transitions to the continuous medium flow. The geometric correction factor can well-describe the gas effective permeability for any shape of pores.
• The gas effective permeability can be slightly higher if the real gas effect is considered. Notably, this real gas effect on the gas effective permeability becomes obvious at high-temperature and pressure conditions, and is thus not ignorable.