# Relative Permeability of Porous Media with Nonuniform Pores

Authors: He Chen; Minghui Yang; Ke Chen; Chao Zhang
Journal: Geofluids
 (2020)

---

## Abstract

Most classical predictive models of relative permeability conceptualize the pores in porous media as assemblies of uniform capillary tubes with different sizes. However, this simplification may overestimate the transport capacity of porous media due to overlooking the effects of the pore nonuniformity. This study presents a simple way to quantify the effect of the nonuniformity of pore cross section on the transport characteristic of unsaturated porous media. The way is based on the index relationship between the porosity of a newly defined reference cross section and that of porous media, which satisfies the intrinsic constraints for the nonuniform porosity of cross sections in porous media. Moreover, the index factor can be captured by a newly defined parameter, called the nonuniformity factor, which is used to establish an extended Darcy’s law. Based on these, a fractal-based continuous analytical model and a fractal-based Monte Carlo model of relative permeability as well as a permeability-porosity model are established. Experimental data of five wetting-nonwetting phase systems, including the water-air, water-steam, water-nitrogen, water-oil, and oil-gas systems, are selected to assess the performance of the proposed model. The results confirm the proposed model’s capacity in capturing the transport properties of various porous media. It is found that the nonuniformity of pores can significantly increase the resistance of fluid flow and thus reduce the transport capacity of porous media.


---

## Body

## 1. Introduction

The relative permeability of both wetting and nonwetting phases is a significant parameter to many engineering fields, e.g., chemical engineering, environmental engineering, and oil and gas reservoir [1, 2]. It is also essential for the numerical simulation of heat and mass transfer in unsaturated porous media [3–6].In the past few decades, numerous studies about the relative permeability for two-phase flow in porous media have been done, containing experimental measurements and model predictions. The steady-state method is a common method for measuring the relative permeability of porous media in the laboratory. It involves that the fixed ratio of the two phases is simultaneously driven in porous media at a given rate and pressure until the sample saturation and pressure differential become constant, e.g., [7–10].Predictive models describing relative permeability-saturation relationships can be divided into empirical models, statistical models, and fractal models [11–13]. The empirical models express the relative permeability of a given phase as a power function of the corresponding phase saturation, which involves the best fitting of existing mathematical formulas to the available experimental data, e.g., [14, 15]. Therefore, these models ignore the pore space characteristics (e.g., tortuosity, pore size distribution, and connectivity) of porous media. Typical representatives of statistical models are the Burdine model and the Mualem model which rely on the pore size distribution function obtained by the water retention curve [16–18]. However, these models usually involve complex integral operations, e.g., [9, 11, 18, 19]. The fractal feature of pores in porous media (e.g., rock, soil, sediment, chemical material, and biological tissue) has been widely recognized and successfully used to study the transport properties of porous media [20, 21]. Unlike the statistical model, since the pore size distribution of porous media is assumed to obey the fractal scaling law, the fractal permeability model can be obtained by Poisson’s law and Darcy’s law rather than the water retention curve, e.g., [22–25].The transport characteristic of multiphase flow in porous media is governed by pore space characteristics, e.g., pore size distribution, the roughness and wettability of pore inner surface, tortuosity, the nonuniformity of pore cross section, and connectivity [9, 26–29]. Compared with the empirical and statistical models, the fractal model is easier to simulate the effect of pore size distribution and pore tortuosity on the permeability of porous media. Moreover, it is also easier to obtain an analytical solution without complicated mathematical operations, which is beneficial to engineering applications. However, existing fractal models typically simplify the pores into parallel and uniform ideal capillary bundles, e.g., [12, 22, 23, 27, 30]. This is suspected of oversimplification, which ignores the effect of the roughness of pore inner surface and the nonuniformity of pore cross section on the transport of multiphase flow. Babadagli et al. [28, 29] have shown that the roughness of the pore inner surface plays an important role in the endpoint of relative permeability curves and in determining the shape of curves. More research on roughness can be found in the literature [31–33]. Nonuniformity of pore cross section can increase the resistance of fluid transport, so ignoring the effect of that may lead to established models that overestimate the transport capacity of porous media. Unfortunately, the effect is rarely discussed and considered in the study of transport properties of porous media.This paper is aimed at studying the effect of pore nonuniformity on the transport properties of unsaturated porous media. To do that, a simple pore model with different cross-sectional sizes is established to simulate the pore nonuniformity. And the nonuniformity can be further integrated into the extended Darcy’s law. Based on fractal theory and Monte Carlo technology, the relative permeability-saturation curves for wetting and nonwetting phases are derived. The proposed model is evaluated by the experimental data of five wetting-nonwetting phase systems (including water-air, water-steam, water-nitrogen, water-oil, and oil-gas systems), and results show good consistency. Additionally, the effect of model parameters on the relative permeability is addressed.


## 2. Theory

### 2.1. The Permeability of Nonuniform Cross-Sectional Pores

The complexity and irregularity of pore structure make it difficult to study the transport properties using natural shape of pores in porous media. As an alternative, the application of simplified pore models (e.g., capillary model and parallel rod model) is acceptable for studying the transport of porous media, e.g., [27, 34, 35]. Although these models are successful in estimating the approximate transport capacity of porous media, they cannot accurately calculate the permeability of wetting and nonwetting phases due to overlooking the nonuniformity of pores. Given this consideration, the capillary model containing different cross-sectional sizes is employed here to simulate the nonuniformity of pores, as shown in Figure 1. In this study, the reference cross section is processed as the research object and associate with other cross sections by parameters a and b. Note that to simplify the model, this study assumes that the capillary parameters of all sizes have the same parameters a and b.Figure 1
A representative elementary volume (REV). (a) The reference cross section of the REV filled by reference radii. (b) The pore space and flow paths in the REV. (c) The simplified capillary containing two different cross-sectional sizes.


#### 2.1.1. Saturated Flow

For capillary flow, the flow resistance of an individual tube can be approximated as follows [36]:
(1)Ω=∫01+aτL08μπr4ldl=8μχπτL0r4,where χ=1+a/b4is a resistance coefficient related to the nonuniformity of pore cross section, where aa>0 and b0<b≤1 are the length factor and size factor of the capillary-throat, respectively, and μ is the fluid viscosity. τ=Lt/L0τ≥1 is the pore tortuosity, where Lt and L0 are the tortuous length and representative length of capillary tubes, respectively. Note that under the assumption of laminar flow, the local energy loss caused by the sudden expansion or contraction of the pore size is negligible compared to the viscous energy loss (see Appendix A); therefore, only the viscous resistance is considered here.Since pores in the REV are simplified as capillary bundles, combining Equation (1), it can calculate the cumulative volumetric flow rate by the Hagen-Poiseuille equation as follows:
(2)Q=∑i=1nΔPΩi=ΔPπ8μχL0∑i=1nri4τi,where n is the total number of pores in the REV, ∆P is the pressure gradient, and subscript i indicates the ith capillary tube.As seen from Figure1, the cross-sectional porosity of the REV is not unique. In this figure, the reference porosity ϕr of the reference cross section is larger than that of the effective porosity ϕ of porous media determined by experiment. If the reference cross section is used as the research object, the relationship between the effective porosity ϕ and the reference porosity ϕr can be expressed as ϕ=ϕrc, where c is defined as the nonuniformity factor. However, in practice, the reference porosity is not necessarily greater than the effective porosity [37], so the range of the nonuniformity factor cshould be c>0, more specifically, when the reference porosity is greater than the effective porosity, c>1; when the reference porosity is less than the effective porosity, 0<c<1; and when the reference porosity equals to the effective porosity, c=1.Therefore, for a given REV with nonuniform pores, Darcy’s law [38] can be extended as follows:
(3)Q=ksΔPL0Apϕ=ksΔPL0Apϕrc,where ks is the permeability of the saturated REV, Ap is the average pore area of the REV cross section, and thus, Ap/ϕ=Ap/ϕrc=AREV, where AREV is the cross-sectional area of the REV. Furthermore, substituting Equation (2) into the rearranged Equation (3) yields the following:
(4)ks=πϕrc8μχAp∑i=1nri4τi.Equation (4) is the permeability model of the saturated REV on a statistical scale, which describes the relationship between the transport characteristics of porous media and pore structure parameters (e.g., porosity ϕ, pore size r, pore tortuosity τ, and the nonuniformity factor c).


#### 2.1.2. Unsaturated Flow

For the unsaturated flow in porous media, the liquid configuration of capillary water can be determined as filling the entire pore and being retained by the tension of meniscus, water pressure, and air pressure. According to the capillary law (also called Young-Laplace equation) [39], small pores are easier to be filled by the wetting phase and larger pores are occupied by the nonwetting phase. Here, the maximum radius of pores filled by the wetting phase is defined as the critical radius rc. Therefore, the volumetric flow rate of the wetting and nonwetting phases in unsaturated porous media can be determined as follows:
(5)Qw=ΔPwπ8μχL0∑i=1rminnwrcri4τi,(6)Qnw=ΔPnwπ8μχL0∑i=nw+1nrmaxri4τi,where ΔPw and ΔPnw are the pressure difference of wetting and nonwetting phases along the REV, respectively.Besides, Buckingham-Darcy’s law [40] describes the transport law of unsaturated flow. Similarly, considering the pore nonuniformity like Equation (3), the extended Buckingham-Darcy’s law can be written as follows [12, 16]:
(7)Qw=kwΔPwL0ApϕrSewc,where kw is the permeability of wetting phase in unsaturated porous media.Substituting Equation (5) into the rearranged Equation (7) yields the following:
(8)kw=πϕrSewc8μχApw∑i=1rminnwrcri4τi.Equation (8) is a generalized permeability function of wetting phase in unsaturated porous media.Additionally, based on the above analysis, the saturation of wetting phase can be determined asSew=∑i=1rminnwrcri2/∑i=1rminnrmaxri2. It should be noted that in this study, Sew is the saturation of capillary tubes and the effective saturation of porous media. Therefore, combining Equations (4) and (8), the relative permeability of wetting phase in porous media can be determined by kr=kw/ks.
(9)krw=Sewc⋅∑i=1rminnwrcri4/τi∑i=1rminnrmaxri4/τi=∑i=1rminnwrcri2∑i=1rminnrmaxri2c⋅∑i=1rminnwrcri4/τi∑i=1rminnrmaxri4/τi.The computation of the relative permeability of the nonwetting phase is identical to that of wetting phase except that the nonwetting phase is considered to occupy pores which size in the range ofrc−rmax.
(10)krn=1−∑i=1rminnwrcri2∑i=1rminnrmaxri2c⋅∑i=nw+1nrmaxri4/τi∑i=1rminnrmaxri4/τi.Equations (9) and (10) are the statistical scale models of relative permeability for wetting and nonwetting phases, respectively. These equations simulate the transport characteristics of porous media with nonuniform pores. However, Equations(9) and (10) cannot directly be used to predict the relative permeability. We need to define the pore size distribution (PSD) of porous media. Indeed, most existing PSDs can be embedded in these equations to develop a specific model, e.g., Gaussian and fractal distributions [12, 22, 38, 41]. In the next section, we will illustrate the application of the proposed statistical scale model by taking the fractal distribution of pore size as an example.


### 2.2. The Analytic Expressions of Relative Permeability Based on Fractal Theory

The pore size distribution of most naturally porous media has been proven to obey the fractal scaling law. In this study, we assumed that the pore size distribution of the reference cross section obeys the fractal scaling law. Therefore, the cumulative number of the pores whose size larger than a measured scaler can be expressed as follows [21, 22]:
(11)Nε≥r=rmaxrDf,where r is the radius of pores in the reference cross section (see Figure 1), N is the cumulative number of capillary tubes whose radius is greater than or equal to r, and Df is the area fractal dimension in the range of 1-2.Since Equation (11) is considered as a continuous and differentiable function, the number of pores in the infinitesimal range r to r+dr can be obtained by differentiating Equation (11) with respect to r as follows:
(12)−dN=DfrmaxDfr−Df+1dr,where dN>0 and the negative sign (-) implies that the number of capillary tubes decreases as the capillary radius increases.


#### 2.2.1. Permeability-Porosity Curve

In order to develop a specific permeability model, the number of pores with a sizer determined by Equation (12) can be substituted into Equation (4) as follows:
(13)ks=π8μχAREV1τa∫rminrmaxr4dNr=18μχτaπAREVDfrmaxDf4−Dfrmax4−Df−rmin4−Df,where τa is the average pore tortuosity of porous media.On the other hand, porosity of the REV can be calculated by substituting Equation (12) into Equation (B.1) (see Appendix B)
(14)ϕr=πAREV∫rminrmaxr2dNr=πAREVDfrmaxDf2−Dfrmax2−Df−rmin2−Df.Here,lnϕr=2−Dflnrmin/rmax [22] can be introduced into Equation (14) yields the following equation:
(15)ϕr=πAREVDf2−Df1−ϕrrmax2.Therefore, the permeability-porosity curve can be determined by substituting Equation (15) into Equation (13), that is,
(16)ks=Cϕ2/c1−ϕ1/c2,where
(17)C=18μχτaAREVπ2−Df24−DfDf.Note that in Equation (16), ϕr=ϕ1/c and rmin/rmax4−Df≈0 have been introduced to simplify the model [22]. Interestingly, Equation (16) is very similar to the Kozeny-Carman (KC) equation [42, 43], which reads as follows:
(18)k=CKCϕ31−ϕ2,where CKC is a fitting parameter related to pore-specific surface area and tortuosity. And the comparison of Equation (16) with the KC equation is discussed further in Section 3.


#### 2.2.2. Relative Permeability-Saturation Curve

Similarly, in order to develop a specific fractal model of relative permeability, the number of pores with a sizer computed by Equation (12) can be substituted into Equation (9).
(19)krw=∫rminrcr2dNr∫rminrmaxr2dNrc⋅∫rminrcr4/τadNr∫rminrmaxr4/τadNr=rc2−Df−rmin2−Dfrmax2−Df−rmin2−Dfc⋅rc4−Df−rmin4−Dfrmax4−Df−rmin4−Df.According to the capillary law,h=C/r, where h is the capillary head and C is a parameter related to the contact angle and liquid surface tension; Equation (19) can also be written as follows:
(20)krw=hcDf−2−hmaxDf−2hminDf−2−hmaxDf−2c⋅hcDf−4−hmaxDf−4hminDf−4−hmaxDf−4,where
(21)hcDf−2−hmaxDf−2hminDf−2−hmaxDf−2=Sew.And substituting Equation (21) into Equation (20), it can obtain the relative permeability krw of the wetting phase expressed by the saturation Sew.
(22)krw=SewcSew1−ϕ1/c+ϕ1/cη,where η=4−Df/2−Df and ϕ1/c=ϕr. Note that lnϕr=2−Dflnrmin/rmax and rmin/rmax4−Df≈0 have been introduced into Equation (22) to simplify the model. Similarly, the relative permeability of the nonwetting phase can be given as follows:
(23)krn=Senc1−1−Sen1−ϕ1/cη1−ϕη/c.Equations (22) and (23) describe the transport characteristics for two-phase flow in the fractal porous media. As seen from Equations (22) and (23) that the relative permeability for both wetting and nonwetting phases is governed by pore structure parameters, i.e., porosity ϕ, the nonuniformity factor c, and the fractal dimension Df.


#### 2.2.3. Fractal-Based Monte Carlo Solutions

As a powerful and unambiguous numerical technique, Monte Carlo simulation is used here to capture the transport characteristics of porous media with nonuniform pores.Based on the fractal scaling law, the total number of pores can be determined by Equation (11).
(24)Ntε≥r=rmax/rminDf.Dividing Equation (12) by Equation (24) gives
(25)−dN/Nt=frdr,where fr=DfrminDfr−Df+1. Yu et al. [23] argued that if the pore radius of the porous media satisfies rmin/rmaxDf≈0, fr can be regarded as a probability density function of pore size. Fortunately, for natural porous media (e.g., soil), rmin/rmax≪10−2 so that rmin/rmaxDf=0 holds approximately. Therefore, the cumulative probability (R), from the smallest radius rmin to the radius r, can be written as follows:
(26)Rr=∫rminrfrdr=1−rminrDf.Equation (20) describes the random probability of the pore size distribution in the range of rmin−rmax, and when r⟶rmin, R⟶0 as well as r⟶rmax, R⟶1.On the other hand, ifR is regarded as a random number from 0 to 1, the pore radius of porous media can be expressed as follows:
(27)ri=rmin1−Ri1/Df,where rmin<r<rmax. Equation (27) is a probability model for the pore size simulated by Monte Carlo method [12, 23, 44]. Therefore, the randomness of pore size distribution can be reproduced by computer-generated random numbers, as shown in Figure 2.Figure 2
Pore sizes simulated by the Monte Carlo method.Substituting Equation (27) into Equations (9) and (10), it can obtain the Monte Carlo solutions of the relative permeability for wetting and nonwetting phases.
(28)krw=∑i=1rmin≤r≤rcnw1/1−Ri2/Df∑i=1rmin≤r≤rmaxn1/1−Ri2/Dfc⋅∑i=1rmin≤r≤rcnw1/τi1−Ri4/Df∑i=1rmin≤r≤rmaxn1/τi1−Ri4/Df,(29)krn=∑i=nw+1rc≤r≤rmaxn1/1−Ri2/Df∑i=1rmin≤r≤rmaxn1/1−Ri2/Dfc⋅∑i=nw+1rc≤r≤rmaxn1/τi1−Ri4/Df∑i=1rmin≤r≤rmaxn1/τi1−Ri4/Df,where
(30)Sew=∑i=1rmin≤r≤rcnw1/1−Ri2/Df∑i=1rmin≤r≤rmaxn1/1−Ri2/Df,and R is a computer-generated random number in the range of 0-1, and the subscript i indicates the ith random number. n is the running total number of the Monte Carlo simulations and the subscript w indicates the number of pores occupied by the wetting phase. The maximum radius of pores occupied by the wetting phase is defined as the critical radius rc. rmin and rmax are the minimum and maximum pore sizes of the REV, respectively. τ is the tortuosity of pores in the REV and the subscript i indicates the ith capillary. Equations (28), (29), and (30) simulate the transport characteristics of porous media with nonuniform pores by the reconstructed pore space based on the Monte Carlo method. Therefore, these equations can capture the relative permeability of the wetting and nonwetting phases.


## 2.1. The Permeability of Nonuniform Cross-Sectional Pores

The complexity and irregularity of pore structure make it difficult to study the transport properties using natural shape of pores in porous media. As an alternative, the application of simplified pore models (e.g., capillary model and parallel rod model) is acceptable for studying the transport of porous media, e.g., [27, 34, 35]. Although these models are successful in estimating the approximate transport capacity of porous media, they cannot accurately calculate the permeability of wetting and nonwetting phases due to overlooking the nonuniformity of pores. Given this consideration, the capillary model containing different cross-sectional sizes is employed here to simulate the nonuniformity of pores, as shown in Figure 1. In this study, the reference cross section is processed as the research object and associate with other cross sections by parameters a and b. Note that to simplify the model, this study assumes that the capillary parameters of all sizes have the same parameters a and b.Figure 1
A representative elementary volume (REV). (a) The reference cross section of the REV filled by reference radii. (b) The pore space and flow paths in the REV. (c) The simplified capillary containing two different cross-sectional sizes.


### 2.1.1. Saturated Flow

For capillary flow, the flow resistance of an individual tube can be approximated as follows [36]:
(1)Ω=∫01+aτL08μπr4ldl=8μχπτL0r4,where χ=1+a/b4is a resistance coefficient related to the nonuniformity of pore cross section, where aa>0 and b0<b≤1 are the length factor and size factor of the capillary-throat, respectively, and μ is the fluid viscosity. τ=Lt/L0τ≥1 is the pore tortuosity, where Lt and L0 are the tortuous length and representative length of capillary tubes, respectively. Note that under the assumption of laminar flow, the local energy loss caused by the sudden expansion or contraction of the pore size is negligible compared to the viscous energy loss (see Appendix A); therefore, only the viscous resistance is considered here.Since pores in the REV are simplified as capillary bundles, combining Equation (1), it can calculate the cumulative volumetric flow rate by the Hagen-Poiseuille equation as follows:
(2)Q=∑i=1nΔPΩi=ΔPπ8μχL0∑i=1nri4τi,where n is the total number of pores in the REV, ∆P is the pressure gradient, and subscript i indicates the ith capillary tube.As seen from Figure1, the cross-sectional porosity of the REV is not unique. In this figure, the reference porosity ϕr of the reference cross section is larger than that of the effective porosity ϕ of porous media determined by experiment. If the reference cross section is used as the research object, the relationship between the effective porosity ϕ and the reference porosity ϕr can be expressed as ϕ=ϕrc, where c is defined as the nonuniformity factor. However, in practice, the reference porosity is not necessarily greater than the effective porosity [37], so the range of the nonuniformity factor cshould be c>0, more specifically, when the reference porosity is greater than the effective porosity, c>1; when the reference porosity is less than the effective porosity, 0<c<1; and when the reference porosity equals to the effective porosity, c=1.Therefore, for a given REV with nonuniform pores, Darcy’s law [38] can be extended as follows:
(3)Q=ksΔPL0Apϕ=ksΔPL0Apϕrc,where ks is the permeability of the saturated REV, Ap is the average pore area of the REV cross section, and thus, Ap/ϕ=Ap/ϕrc=AREV, where AREV is the cross-sectional area of the REV. Furthermore, substituting Equation (2) into the rearranged Equation (3) yields the following:
(4)ks=πϕrc8μχAp∑i=1nri4τi.Equation (4) is the permeability model of the saturated REV on a statistical scale, which describes the relationship between the transport characteristics of porous media and pore structure parameters (e.g., porosity ϕ, pore size r, pore tortuosity τ, and the nonuniformity factor c).


### 2.1.2. Unsaturated Flow

For the unsaturated flow in porous media, the liquid configuration of capillary water can be determined as filling the entire pore and being retained by the tension of meniscus, water pressure, and air pressure. According to the capillary law (also called Young-Laplace equation) [39], small pores are easier to be filled by the wetting phase and larger pores are occupied by the nonwetting phase. Here, the maximum radius of pores filled by the wetting phase is defined as the critical radius rc. Therefore, the volumetric flow rate of the wetting and nonwetting phases in unsaturated porous media can be determined as follows:
(5)Qw=ΔPwπ8μχL0∑i=1rminnwrcri4τi,(6)Qnw=ΔPnwπ8μχL0∑i=nw+1nrmaxri4τi,where ΔPw and ΔPnw are the pressure difference of wetting and nonwetting phases along the REV, respectively.Besides, Buckingham-Darcy’s law [40] describes the transport law of unsaturated flow. Similarly, considering the pore nonuniformity like Equation (3), the extended Buckingham-Darcy’s law can be written as follows [12, 16]:
(7)Qw=kwΔPwL0ApϕrSewc,where kw is the permeability of wetting phase in unsaturated porous media.Substituting Equation (5) into the rearranged Equation (7) yields the following:
(8)kw=πϕrSewc8μχApw∑i=1rminnwrcri4τi.Equation (8) is a generalized permeability function of wetting phase in unsaturated porous media.Additionally, based on the above analysis, the saturation of wetting phase can be determined asSew=∑i=1rminnwrcri2/∑i=1rminnrmaxri2. It should be noted that in this study, Sew is the saturation of capillary tubes and the effective saturation of porous media. Therefore, combining Equations (4) and (8), the relative permeability of wetting phase in porous media can be determined by kr=kw/ks.
(9)krw=Sewc⋅∑i=1rminnwrcri4/τi∑i=1rminnrmaxri4/τi=∑i=1rminnwrcri2∑i=1rminnrmaxri2c⋅∑i=1rminnwrcri4/τi∑i=1rminnrmaxri4/τi.The computation of the relative permeability of the nonwetting phase is identical to that of wetting phase except that the nonwetting phase is considered to occupy pores which size in the range ofrc−rmax.
(10)krn=1−∑i=1rminnwrcri2∑i=1rminnrmaxri2c⋅∑i=nw+1nrmaxri4/τi∑i=1rminnrmaxri4/τi.Equations (9) and (10) are the statistical scale models of relative permeability for wetting and nonwetting phases, respectively. These equations simulate the transport characteristics of porous media with nonuniform pores. However, Equations(9) and (10) cannot directly be used to predict the relative permeability. We need to define the pore size distribution (PSD) of porous media. Indeed, most existing PSDs can be embedded in these equations to develop a specific model, e.g., Gaussian and fractal distributions [12, 22, 38, 41]. In the next section, we will illustrate the application of the proposed statistical scale model by taking the fractal distribution of pore size as an example.


## 2.1.1. Saturated Flow

For capillary flow, the flow resistance of an individual tube can be approximated as follows [36]:
(1)Ω=∫01+aτL08μπr4ldl=8μχπτL0r4,where χ=1+a/b4is a resistance coefficient related to the nonuniformity of pore cross section, where aa>0 and b0<b≤1 are the length factor and size factor of the capillary-throat, respectively, and μ is the fluid viscosity. τ=Lt/L0τ≥1 is the pore tortuosity, where Lt and L0 are the tortuous length and representative length of capillary tubes, respectively. Note that under the assumption of laminar flow, the local energy loss caused by the sudden expansion or contraction of the pore size is negligible compared to the viscous energy loss (see Appendix A); therefore, only the viscous resistance is considered here.Since pores in the REV are simplified as capillary bundles, combining Equation (1), it can calculate the cumulative volumetric flow rate by the Hagen-Poiseuille equation as follows:
(2)Q=∑i=1nΔPΩi=ΔPπ8μχL0∑i=1nri4τi,where n is the total number of pores in the REV, ∆P is the pressure gradient, and subscript i indicates the ith capillary tube.As seen from Figure1, the cross-sectional porosity of the REV is not unique. In this figure, the reference porosity ϕr of the reference cross section is larger than that of the effective porosity ϕ of porous media determined by experiment. If the reference cross section is used as the research object, the relationship between the effective porosity ϕ and the reference porosity ϕr can be expressed as ϕ=ϕrc, where c is defined as the nonuniformity factor. However, in practice, the reference porosity is not necessarily greater than the effective porosity [37], so the range of the nonuniformity factor cshould be c>0, more specifically, when the reference porosity is greater than the effective porosity, c>1; when the reference porosity is less than the effective porosity, 0<c<1; and when the reference porosity equals to the effective porosity, c=1.Therefore, for a given REV with nonuniform pores, Darcy’s law [38] can be extended as follows:
(3)Q=ksΔPL0Apϕ=ksΔPL0Apϕrc,where ks is the permeability of the saturated REV, Ap is the average pore area of the REV cross section, and thus, Ap/ϕ=Ap/ϕrc=AREV, where AREV is the cross-sectional area of the REV. Furthermore, substituting Equation (2) into the rearranged Equation (3) yields the following:
(4)ks=πϕrc8μχAp∑i=1nri4τi.Equation (4) is the permeability model of the saturated REV on a statistical scale, which describes the relationship between the transport characteristics of porous media and pore structure parameters (e.g., porosity ϕ, pore size r, pore tortuosity τ, and the nonuniformity factor c).


## 2.1.2. Unsaturated Flow

For the unsaturated flow in porous media, the liquid configuration of capillary water can be determined as filling the entire pore and being retained by the tension of meniscus, water pressure, and air pressure. According to the capillary law (also called Young-Laplace equation) [39], small pores are easier to be filled by the wetting phase and larger pores are occupied by the nonwetting phase. Here, the maximum radius of pores filled by the wetting phase is defined as the critical radius rc. Therefore, the volumetric flow rate of the wetting and nonwetting phases in unsaturated porous media can be determined as follows:
(5)Qw=ΔPwπ8μχL0∑i=1rminnwrcri4τi,(6)Qnw=ΔPnwπ8μχL0∑i=nw+1nrmaxri4τi,where ΔPw and ΔPnw are the pressure difference of wetting and nonwetting phases along the REV, respectively.Besides, Buckingham-Darcy’s law [40] describes the transport law of unsaturated flow. Similarly, considering the pore nonuniformity like Equation (3), the extended Buckingham-Darcy’s law can be written as follows [12, 16]:
(7)Qw=kwΔPwL0ApϕrSewc,where kw is the permeability of wetting phase in unsaturated porous media.Substituting Equation (5) into the rearranged Equation (7) yields the following:
(8)kw=πϕrSewc8μχApw∑i=1rminnwrcri4τi.Equation (8) is a generalized permeability function of wetting phase in unsaturated porous media.Additionally, based on the above analysis, the saturation of wetting phase can be determined asSew=∑i=1rminnwrcri2/∑i=1rminnrmaxri2. It should be noted that in this study, Sew is the saturation of capillary tubes and the effective saturation of porous media. Therefore, combining Equations (4) and (8), the relative permeability of wetting phase in porous media can be determined by kr=kw/ks.
(9)krw=Sewc⋅∑i=1rminnwrcri4/τi∑i=1rminnrmaxri4/τi=∑i=1rminnwrcri2∑i=1rminnrmaxri2c⋅∑i=1rminnwrcri4/τi∑i=1rminnrmaxri4/τi.The computation of the relative permeability of the nonwetting phase is identical to that of wetting phase except that the nonwetting phase is considered to occupy pores which size in the range ofrc−rmax.
(10)krn=1−∑i=1rminnwrcri2∑i=1rminnrmaxri2c⋅∑i=nw+1nrmaxri4/τi∑i=1rminnrmaxri4/τi.Equations (9) and (10) are the statistical scale models of relative permeability for wetting and nonwetting phases, respectively. These equations simulate the transport characteristics of porous media with nonuniform pores. However, Equations(9) and (10) cannot directly be used to predict the relative permeability. We need to define the pore size distribution (PSD) of porous media. Indeed, most existing PSDs can be embedded in these equations to develop a specific model, e.g., Gaussian and fractal distributions [12, 22, 38, 41]. In the next section, we will illustrate the application of the proposed statistical scale model by taking the fractal distribution of pore size as an example.


## 2.2. The Analytic Expressions of Relative Permeability Based on Fractal Theory

The pore size distribution of most naturally porous media has been proven to obey the fractal scaling law. In this study, we assumed that the pore size distribution of the reference cross section obeys the fractal scaling law. Therefore, the cumulative number of the pores whose size larger than a measured scaler can be expressed as follows [21, 22]:
(11)Nε≥r=rmaxrDf,where r is the radius of pores in the reference cross section (see Figure 1), N is the cumulative number of capillary tubes whose radius is greater than or equal to r, and Df is the area fractal dimension in the range of 1-2.Since Equation (11) is considered as a continuous and differentiable function, the number of pores in the infinitesimal range r to r+dr can be obtained by differentiating Equation (11) with respect to r as follows:
(12)−dN=DfrmaxDfr−Df+1dr,where dN>0 and the negative sign (-) implies that the number of capillary tubes decreases as the capillary radius increases.


### 2.2.1. Permeability-Porosity Curve

In order to develop a specific permeability model, the number of pores with a sizer determined by Equation (12) can be substituted into Equation (4) as follows:
(13)ks=π8μχAREV1τa∫rminrmaxr4dNr=18μχτaπAREVDfrmaxDf4−Dfrmax4−Df−rmin4−Df,where τa is the average pore tortuosity of porous media.On the other hand, porosity of the REV can be calculated by substituting Equation (12) into Equation (B.1) (see Appendix B)
(14)ϕr=πAREV∫rminrmaxr2dNr=πAREVDfrmaxDf2−Dfrmax2−Df−rmin2−Df.Here,lnϕr=2−Dflnrmin/rmax [22] can be introduced into Equation (14) yields the following equation:
(15)ϕr=πAREVDf2−Df1−ϕrrmax2.Therefore, the permeability-porosity curve can be determined by substituting Equation (15) into Equation (13), that is,
(16)ks=Cϕ2/c1−ϕ1/c2,where
(17)C=18μχτaAREVπ2−Df24−DfDf.Note that in Equation (16), ϕr=ϕ1/c and rmin/rmax4−Df≈0 have been introduced to simplify the model [22]. Interestingly, Equation (16) is very similar to the Kozeny-Carman (KC) equation [42, 43], which reads as follows:
(18)k=CKCϕ31−ϕ2,where CKC is a fitting parameter related to pore-specific surface area and tortuosity. And the comparison of Equation (16) with the KC equation is discussed further in Section 3.


### 2.2.2. Relative Permeability-Saturation Curve

Similarly, in order to develop a specific fractal model of relative permeability, the number of pores with a sizer computed by Equation (12) can be substituted into Equation (9).
(19)krw=∫rminrcr2dNr∫rminrmaxr2dNrc⋅∫rminrcr4/τadNr∫rminrmaxr4/τadNr=rc2−Df−rmin2−Dfrmax2−Df−rmin2−Dfc⋅rc4−Df−rmin4−Dfrmax4−Df−rmin4−Df.According to the capillary law,h=C/r, where h is the capillary head and C is a parameter related to the contact angle and liquid surface tension; Equation (19) can also be written as follows:
(20)krw=hcDf−2−hmaxDf−2hminDf−2−hmaxDf−2c⋅hcDf−4−hmaxDf−4hminDf−4−hmaxDf−4,where
(21)hcDf−2−hmaxDf−2hminDf−2−hmaxDf−2=Sew.And substituting Equation (21) into Equation (20), it can obtain the relative permeability krw of the wetting phase expressed by the saturation Sew.
(22)krw=SewcSew1−ϕ1/c+ϕ1/cη,where η=4−Df/2−Df and ϕ1/c=ϕr. Note that lnϕr=2−Dflnrmin/rmax and rmin/rmax4−Df≈0 have been introduced into Equation (22) to simplify the model. Similarly, the relative permeability of the nonwetting phase can be given as follows:
(23)krn=Senc1−1−Sen1−ϕ1/cη1−ϕη/c.Equations (22) and (23) describe the transport characteristics for two-phase flow in the fractal porous media. As seen from Equations (22) and (23) that the relative permeability for both wetting and nonwetting phases is governed by pore structure parameters, i.e., porosity ϕ, the nonuniformity factor c, and the fractal dimension Df.


### 2.2.3. Fractal-Based Monte Carlo Solutions

As a powerful and unambiguous numerical technique, Monte Carlo simulation is used here to capture the transport characteristics of porous media with nonuniform pores.Based on the fractal scaling law, the total number of pores can be determined by Equation (11).
(24)Ntε≥r=rmax/rminDf.Dividing Equation (12) by Equation (24) gives
(25)−dN/Nt=frdr,where fr=DfrminDfr−Df+1. Yu et al. [23] argued that if the pore radius of the porous media satisfies rmin/rmaxDf≈0, fr can be regarded as a probability density function of pore size. Fortunately, for natural porous media (e.g., soil), rmin/rmax≪10−2 so that rmin/rmaxDf=0 holds approximately. Therefore, the cumulative probability (R), from the smallest radius rmin to the radius r, can be written as follows:
(26)Rr=∫rminrfrdr=1−rminrDf.Equation (20) describes the random probability of the pore size distribution in the range of rmin−rmax, and when r⟶rmin, R⟶0 as well as r⟶rmax, R⟶1.On the other hand, ifR is regarded as a random number from 0 to 1, the pore radius of porous media can be expressed as follows:
(27)ri=rmin1−Ri1/Df,where rmin<r<rmax. Equation (27) is a probability model for the pore size simulated by Monte Carlo method [12, 23, 44]. Therefore, the randomness of pore size distribution can be reproduced by computer-generated random numbers, as shown in Figure 2.Figure 2
Pore sizes simulated by the Monte Carlo method.Substituting Equation (27) into Equations (9) and (10), it can obtain the Monte Carlo solutions of the relative permeability for wetting and nonwetting phases.
(28)krw=∑i=1rmin≤r≤rcnw1/1−Ri2/Df∑i=1rmin≤r≤rmaxn1/1−Ri2/Dfc⋅∑i=1rmin≤r≤rcnw1/τi1−Ri4/Df∑i=1rmin≤r≤rmaxn1/τi1−Ri4/Df,(29)krn=∑i=nw+1rc≤r≤rmaxn1/1−Ri2/Df∑i=1rmin≤r≤rmaxn1/1−Ri2/Dfc⋅∑i=nw+1rc≤r≤rmaxn1/τi1−Ri4/Df∑i=1rmin≤r≤rmaxn1/τi1−Ri4/Df,where
(30)Sew=∑i=1rmin≤r≤rcnw1/1−Ri2/Df∑i=1rmin≤r≤rmaxn1/1−Ri2/Df,and R is a computer-generated random number in the range of 0-1, and the subscript i indicates the ith random number. n is the running total number of the Monte Carlo simulations and the subscript w indicates the number of pores occupied by the wetting phase. The maximum radius of pores occupied by the wetting phase is defined as the critical radius rc. rmin and rmax are the minimum and maximum pore sizes of the REV, respectively. τ is the tortuosity of pores in the REV and the subscript i indicates the ith capillary. Equations (28), (29), and (30) simulate the transport characteristics of porous media with nonuniform pores by the reconstructed pore space based on the Monte Carlo method. Therefore, these equations can capture the relative permeability of the wetting and nonwetting phases.


## 2.2.1. Permeability-Porosity Curve

In order to develop a specific permeability model, the number of pores with a sizer determined by Equation (12) can be substituted into Equation (4) as follows:
(13)ks=π8μχAREV1τa∫rminrmaxr4dNr=18μχτaπAREVDfrmaxDf4−Dfrmax4−Df−rmin4−Df,where τa is the average pore tortuosity of porous media.On the other hand, porosity of the REV can be calculated by substituting Equation (12) into Equation (B.1) (see Appendix B)
(14)ϕr=πAREV∫rminrmaxr2dNr=πAREVDfrmaxDf2−Dfrmax2−Df−rmin2−Df.Here,lnϕr=2−Dflnrmin/rmax [22] can be introduced into Equation (14) yields the following equation:
(15)ϕr=πAREVDf2−Df1−ϕrrmax2.Therefore, the permeability-porosity curve can be determined by substituting Equation (15) into Equation (13), that is,
(16)ks=Cϕ2/c1−ϕ1/c2,where
(17)C=18μχτaAREVπ2−Df24−DfDf.Note that in Equation (16), ϕr=ϕ1/c and rmin/rmax4−Df≈0 have been introduced to simplify the model [22]. Interestingly, Equation (16) is very similar to the Kozeny-Carman (KC) equation [42, 43], which reads as follows:
(18)k=CKCϕ31−ϕ2,where CKC is a fitting parameter related to pore-specific surface area and tortuosity. And the comparison of Equation (16) with the KC equation is discussed further in Section 3.


## 2.2.2. Relative Permeability-Saturation Curve

Similarly, in order to develop a specific fractal model of relative permeability, the number of pores with a sizer computed by Equation (12) can be substituted into Equation (9).
(19)krw=∫rminrcr2dNr∫rminrmaxr2dNrc⋅∫rminrcr4/τadNr∫rminrmaxr4/τadNr=rc2−Df−rmin2−Dfrmax2−Df−rmin2−Dfc⋅rc4−Df−rmin4−Dfrmax4−Df−rmin4−Df.According to the capillary law,h=C/r, where h is the capillary head and C is a parameter related to the contact angle and liquid surface tension; Equation (19) can also be written as follows:
(20)krw=hcDf−2−hmaxDf−2hminDf−2−hmaxDf−2c⋅hcDf−4−hmaxDf−4hminDf−4−hmaxDf−4,where
(21)hcDf−2−hmaxDf−2hminDf−2−hmaxDf−2=Sew.And substituting Equation (21) into Equation (20), it can obtain the relative permeability krw of the wetting phase expressed by the saturation Sew.
(22)krw=SewcSew1−ϕ1/c+ϕ1/cη,where η=4−Df/2−Df and ϕ1/c=ϕr. Note that lnϕr=2−Dflnrmin/rmax and rmin/rmax4−Df≈0 have been introduced into Equation (22) to simplify the model. Similarly, the relative permeability of the nonwetting phase can be given as follows:
(23)krn=Senc1−1−Sen1−ϕ1/cη1−ϕη/c.Equations (22) and (23) describe the transport characteristics for two-phase flow in the fractal porous media. As seen from Equations (22) and (23) that the relative permeability for both wetting and nonwetting phases is governed by pore structure parameters, i.e., porosity ϕ, the nonuniformity factor c, and the fractal dimension Df.


## 2.2.3. Fractal-Based Monte Carlo Solutions

As a powerful and unambiguous numerical technique, Monte Carlo simulation is used here to capture the transport characteristics of porous media with nonuniform pores.Based on the fractal scaling law, the total number of pores can be determined by Equation (11).
(24)Ntε≥r=rmax/rminDf.Dividing Equation (12) by Equation (24) gives
(25)−dN/Nt=frdr,where fr=DfrminDfr−Df+1. Yu et al. [23] argued that if the pore radius of the porous media satisfies rmin/rmaxDf≈0, fr can be regarded as a probability density function of pore size. Fortunately, for natural porous media (e.g., soil), rmin/rmax≪10−2 so that rmin/rmaxDf=0 holds approximately. Therefore, the cumulative probability (R), from the smallest radius rmin to the radius r, can be written as follows:
(26)Rr=∫rminrfrdr=1−rminrDf.Equation (20) describes the random probability of the pore size distribution in the range of rmin−rmax, and when r⟶rmin, R⟶0 as well as r⟶rmax, R⟶1.On the other hand, ifR is regarded as a random number from 0 to 1, the pore radius of porous media can be expressed as follows:
(27)ri=rmin1−Ri1/Df,where rmin<r<rmax. Equation (27) is a probability model for the pore size simulated by Monte Carlo method [12, 23, 44]. Therefore, the randomness of pore size distribution can be reproduced by computer-generated random numbers, as shown in Figure 2.Figure 2
Pore sizes simulated by the Monte Carlo method.Substituting Equation (27) into Equations (9) and (10), it can obtain the Monte Carlo solutions of the relative permeability for wetting and nonwetting phases.
(28)krw=∑i=1rmin≤r≤rcnw1/1−Ri2/Df∑i=1rmin≤r≤rmaxn1/1−Ri2/Dfc⋅∑i=1rmin≤r≤rcnw1/τi1−Ri4/Df∑i=1rmin≤r≤rmaxn1/τi1−Ri4/Df,(29)krn=∑i=nw+1rc≤r≤rmaxn1/1−Ri2/Df∑i=1rmin≤r≤rmaxn1/1−Ri2/Dfc⋅∑i=nw+1rc≤r≤rmaxn1/τi1−Ri4/Df∑i=1rmin≤r≤rmaxn1/τi1−Ri4/Df,where
(30)Sew=∑i=1rmin≤r≤rcnw1/1−Ri2/Df∑i=1rmin≤r≤rmaxn1/1−Ri2/Df,and R is a computer-generated random number in the range of 0-1, and the subscript i indicates the ith random number. n is the running total number of the Monte Carlo simulations and the subscript w indicates the number of pores occupied by the wetting phase. The maximum radius of pores occupied by the wetting phase is defined as the critical radius rc. rmin and rmax are the minimum and maximum pore sizes of the REV, respectively. τ is the tortuosity of pores in the REV and the subscript i indicates the ith capillary. Equations (28), (29), and (30) simulate the transport characteristics of porous media with nonuniform pores by the reconstructed pore space based on the Monte Carlo method. Therefore, these equations can capture the relative permeability of the wetting and nonwetting phases.


## 3. Applications

### 3.1. Permeability-Porosity Curve

Fractal-based permeability-porosity curve (Equation (16)) contains two parameters, i.e., C and the nonuniformity factor c. In order to evaluate the proposed permeability model, three data sets are selected, i.e., fine-grained sandstone, silty sandstones, and Timimoun Basin (“tight gas” sandstones) from Chilindar [45] and Hirst et al. [46], respectively. For each type of porous media, Equation (16) is fitted to the corresponding data set by calculating the minimum root-mean-square error (RMSE).
(31)RMSE=1M∑j=1Mlnkr,jmeas.−lnkr,jpred.2,where M is the number of the test data, kr,jmeas. is the jth measured data point, and kr,jpred. is the corresponding predicted relative permeability. Here, a logarithmic scale is applied to accommodate a wide range of variations in permeability. Besides, in this section, the KC equation (Equation (17)) is also compared with the proposed model. The fitted parameters and the corresponding RMSE of Equation (16) and the KC equation are listed in Table 1. It can be seen that the RMSE of Equation (16) is smaller than that of the KC equation for each material. Figure 3 shows that the proposed permeability-porosity model predicts good the measured data over a wide range of magnitude (about 4 to 11 orders of magnitude).Table 1
Values of the fitted parameters and the corresponding RMSE.
Porous media
c
C
RMSE
CKC
RMSE
Reference
Fine-grained sandstone
0.347
7.8×106
0.282
3.4×104
0.912
Chilindar [45]
Silty sandstones
0.382
5.3×105
0.402
7.7×103
0.648
Chilindar [45]
Timimoun Basin
0.395
3.1×105
1.794
1.6×102
2.811
Hirst et al. [46]Figure 3
The predicted permeability-porosity curve and the comparison to the KC equation and measured data.
(a)
(b)
(c)


### 3.2. Relative Permeability-Saturation Curve

#### 3.2.1. Methods

Fractal-based analytic expressions (Equations (22) and (23)) of the relative permeability contain three parameters, i.e., c, ϕ, and Df. ϕ is the porosity of porous media determined by experiments, c is the nonuniformity factor of pore cross section, and Df is the fractal dimension determined by the best-fitting or the box-counting method [47, 48]. In this section, model parameters are also determined by fitting measured data and calculating the minimum RMSE between prediction values and experimental data.
(32)RMSE=1M∑j=1Mkr,jmeas.−kr,jpred.2.Fractal-based Monte Carlo solutions (Equations (28) and (29)) of the relative permeability contain parameters c, Df, and τi. In this study, Monte Carlo method simulates the randomness of pore sizes based on fractal in porous media. The algorithm for the determination of the relative permeability for unsaturated porous media is summarized as follows:
(1)
Give a porosityϕ and find a nonuniformity factor c of pore cross section
(2)
Find fractal dimensionDf, rmin, rmax and calculate the running total number n of the Monte Carlo simulations according to Equation (24)
(3)
Generate a random numberRi between 0 and 1
(4)
Calculateri by Equation (27) and rmin≤ri≤rmax; otherwise, return to procedure (3)
(5)
Generate random tortuosityτi of 1.25-1.78 suggested by Xu et al. [12]
(6)
Calculate the saturationSew and relative permeability kr (krw and krn)
(7)
CalculateX times kr,x at a same saturation Sew, and then find an average value kr according to kr=1/X∑x=1Xkr,x
(8)
Calibrate model parameters by comparing prediction data with experimental dataAs seen from the above algorithm, based on fractal theory, procedures (1)-(5) reconstruct the pore size distribution in porous media. Evaluation of the model and comparison with other models will be further analyzed and discussed in the next sections.


#### 3.2.2. Results and Discussion

Eight data sets are selected for evaluating the proposed model. These data sets involve five fluid systems, i.e., water-air, water-nitrogen, water-steam, water-oil, and oil-gas systems. The pore structure parameters of these porous media are listed in Table2.Table 2
The pore structure parameters of porous media and the corresponding RMSE.
System (wetting-nonwetting)
Porous media
ϕ
Df
c
RMSE
Reference
krn
krw
Water-air
Hygiene sandstone
0.250
1.21
1.98
0.012
0.074
Brooks and Corey [7]
Water-air
Berea sandstone
0.206
1.14
3.73
0.041
0.059
Brooks and Corey [7]
Water-air
Grenoble sand
0.370
1.75
2.18
0.050
0.034
Touma and Vauclin [49]
Water-air
Oso Flaco sand
0.400
1.45
1.90
0.026
0.012
Tuli and Hopmans [9]
Water-nitrogen
Berea core
0.240
1.23
2.40
0.033
0.032
Li and Horne [10]
Water-steam
Core
0.250
1.77
1.35
0.108
0.020
O’Connor [8]
Water-oil
Berea core
0.225
1.18
1.96
0.021
0.110
Kleppe and Morse [50]
Oil-gas
Pyrex core
0.374
1.66
2.18
0.012
0.054
Gates and Leitz [51]Figure4 shows measured data of the water-air system and predicted relative permeability curves by the fractal analysis and Monte Carlo method. These data obtained by the instantaneous profile method in Grenoble sand, Berea sandstone, Hygiene sandstone, and Oso Flaco sand [7, 9, 49]. As shown, the predicted relative permeability (including wetting and nonwetting phases) is in good agreement with the experimental data.Figure 4
Predicted relative permeability for water-air system and the comparison to the experimental data.
(a)
(b)
(c)
(d)For water-nitrogen and water-steam systems, Figure5 shows that both the Monte Carlo solutions and the fractal analytical solutions make a good prediction of the relative permeability, except for the steam phase. This may attribute to the instability of water-steam system in which the immobile steam easily condenses into water so that there is a defect inaccuracy in determining the steam saturation [8].Figure 5
Predicted relative permeability for (a) water-nitrogen system and (b) water-steam system as well as the comparison to the experimental data.
(a)
(b)Studying the relative permeability of water-oil and oil-gas systems is significant for contamination in oil reservoirs as well as in soils. Figure6(a) shows the predicted relative permeability curves of the oil-water system against the corresponding experimental data [50]. The predicted relative permeability of nonwetting phase (oil) is in good agreement with the experimental data while that of wetting phase (water) is not such a good match. The prediction bias may be due to the fact that oil and water are coupled in terms of viscosity or tension, which makes the wetting phase permeability insensitive to changes in saturation. As seen in Figure 6(b), for the oil-gas system, the proposed model can also well determine the relative permeability of both wetting and nonwetting phase in the full range of saturation in Pyrex core [51].Figure 6
Predicted relative permeability for (a) water-oil system and (b) oil-gas system as well as the comparison to the experimental data.
(a)
(b)On the whole, these examples illustrate that the proposed physics-based model, which considers the nonuniformity of pore cross section, can generally capture the transport characteristics for both wetting and nonwetting phase in various porous media. It should be noted that due to the power-law feature of the fractal distribution of pore size, the proposed fractal models in this paper may not be applicable to the porous media of multimodal pore structure, which involves the multipeak distribution of pore size [41].


#### 3.2.3. Comparison with Other Models

To further verify the validity of the developed model, the proposed model is compared with the other two models. In the Brook and Corey-Mualem model (BC-M) and the van Genuchten-Mualem model (VG-M) [18, 52], the relative permeability of wetting and nonwetting phases is, respectively, calculated by BC-M:
(33)krw=Sew2.5+2/λ,krn=1−Sew0.51−Sew1+1/λ2,VG-M:(34)krw=Sew0.51−1−Sew1/mm2,krn=1−Sew0.51−Sew1/m2m,where λ is a model parameter related to the pore size distribution of porous media and m is a fitted parameter.Similarly, RMSE is used for the evaluation of these models. Table3 summarizes the RMSE of krw and krn predicted by each model (BC-M, VG-M, and proposed models) for various porous media. For the sake of interpretation, the smallest values of RMSE are highlighted in bold. As seen from Table 3, compared with the BC-M and VG-M models, the proposed model performs best, occupying eleven out of sixteen highlighted values. This overall performance can be reflected by the average values of RMSE (aveRMSE), which are only 0.049 and 0.038 for wetting and nonwetting phases, respectively.Table 3
RMSE values obtained with the BC-M, VG-M, and proposed models.
System (wetting-nonwetting)
Porous media
RMSE (BC-M)
RMSE (VG-M)
RMSE (model)
krn
krw
krn
krw
krn
krw
Water-air
Hygiene sandstone
0.022
0.070
0.020
0.072
0.012
0.074
Water-air
Berea sandstone
0.201
0.096
0.090
0.183
0.041
0.059
Water-air
Grenoble sand
0.112
0.116
0.110
0.079
0.050
0.034
Water-air
Oso Flaco sand
0.082
0.014
0.041
0.030
0.026
0.012
Water-nitrogen
Berea core
0.054
0.090
0.042
0.103
0.033
0.032
Water-steam
Core
0.103
0.122
0.123
0.052
0.108
0.020
Water-oil
Berea core
0.028
0.038
0.016
0.032
0.021
0.110
Oil-gas
Pyrex core
0.019
0.119
0.009
0.128
0.012
0.054
0.078
0.083
0.056
0.085
0.038
0.049Additionally, the differences between these models can be clearly visualized as scatter plots of measured and predicted relative permeability values shown in Figure7. As shown, the BC-M and VG-M models tend to generally overestimate the transport capacity of porous media over the full range of saturation, which may be attributed to the fact that the nonuniformity of pores was overlooked. In contrast, the predicted relative permeability values (both krw and krn) of the proposed model are close to the measured values, and these data sets show a much better linearship.Figure 7
Scatter plot of measured versus predicted values for BC-M, VG-M, and the proposed models.Note that thin-film flow is not considered in this paper, and the film flow is usually distributed on the particle surface or the inner surface of the capillary to form an annular fluid structure [53]. Therefore, the continuous capillary model developed in this study may overestimate the krn of the nonwetting phase and underestimate the krw of the wetting phase, especially in the low relative permeability stage, which can be seen in Figure 7.


#### 3.2.4. Parameter Analysis

In this section, the effect of model parametersDf and c on the relative permeability of porous media is also addressed. As seen from Figure 8, the relative permeability of wetting phases decreases as the fractal dimension Df increases, while enhances that of nonwetting phase. This is because a larger Df means a larger percent of small pores, which will decrease the critical radius with a given wetting saturation. Therefore, the transport of wetting phase requires a larger pressure gradient, which reduces its relative permeability and increases the relative permeability of nonwetting phase.Figure 8
The effect of the pore size distribution on relative permeability of porous media.From Figure9(a), the relative permeability of both wetting and nonwetting phases decreases as the nonuniformity factor c increases. Under the assumption of laminar flow, the pressure drop per unit length caused by the flow resistance can be obtained (see Appendix B).
(35)ΔPL′=4μτr21+ϕ2−2/c2−ϕ1/c−12vs.Figure 9
The effect of nonuniformity factorcon the relative permeability of porous media and pressure drop (μ=0.8937×10−3Pas, 25°C, r=1×10−3m, vs=1×10−3m/s, τ=1.5).
(a)
(b)And Figure9(b) shows that pressure drop increases significantly with the increase of c, which means that increasing in the nonuniformity factor c can enhance the fluid resistance and thus lower the transport capacity of porous media.


## 3.1. Permeability-Porosity Curve

Fractal-based permeability-porosity curve (Equation (16)) contains two parameters, i.e., C and the nonuniformity factor c. In order to evaluate the proposed permeability model, three data sets are selected, i.e., fine-grained sandstone, silty sandstones, and Timimoun Basin (“tight gas” sandstones) from Chilindar [45] and Hirst et al. [46], respectively. For each type of porous media, Equation (16) is fitted to the corresponding data set by calculating the minimum root-mean-square error (RMSE).
(31)RMSE=1M∑j=1Mlnkr,jmeas.−lnkr,jpred.2,where M is the number of the test data, kr,jmeas. is the jth measured data point, and kr,jpred. is the corresponding predicted relative permeability. Here, a logarithmic scale is applied to accommodate a wide range of variations in permeability. Besides, in this section, the KC equation (Equation (17)) is also compared with the proposed model. The fitted parameters and the corresponding RMSE of Equation (16) and the KC equation are listed in Table 1. It can be seen that the RMSE of Equation (16) is smaller than that of the KC equation for each material. Figure 3 shows that the proposed permeability-porosity model predicts good the measured data over a wide range of magnitude (about 4 to 11 orders of magnitude).Table 1
Values of the fitted parameters and the corresponding RMSE.
Porous media
c
C
RMSE
CKC
RMSE
Reference
Fine-grained sandstone
0.347
7.8×106
0.282
3.4×104
0.912
Chilindar [45]
Silty sandstones
0.382
5.3×105
0.402
7.7×103
0.648
Chilindar [45]
Timimoun Basin
0.395
3.1×105
1.794
1.6×102
2.811
Hirst et al. [46]Figure 3
The predicted permeability-porosity curve and the comparison to the KC equation and measured data.
(a)
(b)
(c)


## 3.2. Relative Permeability-Saturation Curve

### 3.2.1. Methods

Fractal-based analytic expressions (Equations (22) and (23)) of the relative permeability contain three parameters, i.e., c, ϕ, and Df. ϕ is the porosity of porous media determined by experiments, c is the nonuniformity factor of pore cross section, and Df is the fractal dimension determined by the best-fitting or the box-counting method [47, 48]. In this section, model parameters are also determined by fitting measured data and calculating the minimum RMSE between prediction values and experimental data.
(32)RMSE=1M∑j=1Mkr,jmeas.−kr,jpred.2.Fractal-based Monte Carlo solutions (Equations (28) and (29)) of the relative permeability contain parameters c, Df, and τi. In this study, Monte Carlo method simulates the randomness of pore sizes based on fractal in porous media. The algorithm for the determination of the relative permeability for unsaturated porous media is summarized as follows:
(1)
Give a porosityϕ and find a nonuniformity factor c of pore cross section
(2)
Find fractal dimensionDf, rmin, rmax and calculate the running total number n of the Monte Carlo simulations according to Equation (24)
(3)
Generate a random numberRi between 0 and 1
(4)
Calculateri by Equation (27) and rmin≤ri≤rmax; otherwise, return to procedure (3)
(5)
Generate random tortuosityτi of 1.25-1.78 suggested by Xu et al. [12]
(6)
Calculate the saturationSew and relative permeability kr (krw and krn)
(7)
CalculateX times kr,x at a same saturation Sew, and then find an average value kr according to kr=1/X∑x=1Xkr,x
(8)
Calibrate model parameters by comparing prediction data with experimental dataAs seen from the above algorithm, based on fractal theory, procedures (1)-(5) reconstruct the pore size distribution in porous media. Evaluation of the model and comparison with other models will be further analyzed and discussed in the next sections.


### 3.2.2. Results and Discussion

Eight data sets are selected for evaluating the proposed model. These data sets involve five fluid systems, i.e., water-air, water-nitrogen, water-steam, water-oil, and oil-gas systems. The pore structure parameters of these porous media are listed in Table2.Table 2
The pore structure parameters of porous media and the corresponding RMSE.
System (wetting-nonwetting)
Porous media
ϕ
Df
c
RMSE
Reference
krn
krw
Water-air
Hygiene sandstone
0.250
1.21
1.98
0.012
0.074
Brooks and Corey [7]
Water-air
Berea sandstone
0.206
1.14
3.73
0.041
0.059
Brooks and Corey [7]
Water-air
Grenoble sand
0.370
1.75
2.18
0.050
0.034
Touma and Vauclin [49]
Water-air
Oso Flaco sand
0.400
1.45
1.90
0.026
0.012
Tuli and Hopmans [9]
Water-nitrogen
Berea core
0.240
1.23
2.40
0.033
0.032
Li and Horne [10]
Water-steam
Core
0.250
1.77
1.35
0.108
0.020
O’Connor [8]
Water-oil
Berea core
0.225
1.18
1.96
0.021
0.110
Kleppe and Morse [50]
Oil-gas
Pyrex core
0.374
1.66
2.18
0.012
0.054
Gates and Leitz [51]Figure4 shows measured data of the water-air system and predicted relative permeability curves by the fractal analysis and Monte Carlo method. These data obtained by the instantaneous profile method in Grenoble sand, Berea sandstone, Hygiene sandstone, and Oso Flaco sand [7, 9, 49]. As shown, the predicted relative permeability (including wetting and nonwetting phases) is in good agreement with the experimental data.Figure 4
Predicted relative permeability for water-air system and the comparison to the experimental data.
(a)
(b)
(c)
(d)For water-nitrogen and water-steam systems, Figure5 shows that both the Monte Carlo solutions and the fractal analytical solutions make a good prediction of the relative permeability, except for the steam phase. This may attribute to the instability of water-steam system in which the immobile steam easily condenses into water so that there is a defect inaccuracy in determining the steam saturation [8].Figure 5
Predicted relative permeability for (a) water-nitrogen system and (b) water-steam system as well as the comparison to the experimental data.
(a)
(b)Studying the relative permeability of water-oil and oil-gas systems is significant for contamination in oil reservoirs as well as in soils. Figure6(a) shows the predicted relative permeability curves of the oil-water system against the corresponding experimental data [50]. The predicted relative permeability of nonwetting phase (oil) is in good agreement with the experimental data while that of wetting phase (water) is not such a good match. The prediction bias may be due to the fact that oil and water are coupled in terms of viscosity or tension, which makes the wetting phase permeability insensitive to changes in saturation. As seen in Figure 6(b), for the oil-gas system, the proposed model can also well determine the relative permeability of both wetting and nonwetting phase in the full range of saturation in Pyrex core [51].Figure 6
Predicted relative permeability for (a) water-oil system and (b) oil-gas system as well as the comparison to the experimental data.
(a)
(b)On the whole, these examples illustrate that the proposed physics-based model, which considers the nonuniformity of pore cross section, can generally capture the transport characteristics for both wetting and nonwetting phase in various porous media. It should be noted that due to the power-law feature of the fractal distribution of pore size, the proposed fractal models in this paper may not be applicable to the porous media of multimodal pore structure, which involves the multipeak distribution of pore size [41].


### 3.2.3. Comparison with Other Models

To further verify the validity of the developed model, the proposed model is compared with the other two models. In the Brook and Corey-Mualem model (BC-M) and the van Genuchten-Mualem model (VG-M) [18, 52], the relative permeability of wetting and nonwetting phases is, respectively, calculated by BC-M:
(33)krw=Sew2.5+2/λ,krn=1−Sew0.51−Sew1+1/λ2,VG-M:(34)krw=Sew0.51−1−Sew1/mm2,krn=1−Sew0.51−Sew1/m2m,where λ is a model parameter related to the pore size distribution of porous media and m is a fitted parameter.Similarly, RMSE is used for the evaluation of these models. Table3 summarizes the RMSE of krw and krn predicted by each model (BC-M, VG-M, and proposed models) for various porous media. For the sake of interpretation, the smallest values of RMSE are highlighted in bold. As seen from Table 3, compared with the BC-M and VG-M models, the proposed model performs best, occupying eleven out of sixteen highlighted values. This overall performance can be reflected by the average values of RMSE (aveRMSE), which are only 0.049 and 0.038 for wetting and nonwetting phases, respectively.Table 3
RMSE values obtained with the BC-M, VG-M, and proposed models.
System (wetting-nonwetting)
Porous media
RMSE (BC-M)
RMSE (VG-M)
RMSE (model)
krn
krw
krn
krw
krn
krw
Water-air
Hygiene sandstone
0.022
0.070
0.020
0.072
0.012
0.074
Water-air
Berea sandstone
0.201
0.096
0.090
0.183
0.041
0.059
Water-air
Grenoble sand
0.112
0.116
0.110
0.079
0.050
0.034
Water-air
Oso Flaco sand
0.082
0.014
0.041
0.030
0.026
0.012
Water-nitrogen
Berea core
0.054
0.090
0.042
0.103
0.033
0.032
Water-steam
Core
0.103
0.122
0.123
0.052
0.108
0.020
Water-oil
Berea core
0.028
0.038
0.016
0.032
0.021
0.110
Oil-gas
Pyrex core
0.019
0.119
0.009
0.128
0.012
0.054
0.078
0.083
0.056
0.085
0.038
0.049Additionally, the differences between these models can be clearly visualized as scatter plots of measured and predicted relative permeability values shown in Figure7. As shown, the BC-M and VG-M models tend to generally overestimate the transport capacity of porous media over the full range of saturation, which may be attributed to the fact that the nonuniformity of pores was overlooked. In contrast, the predicted relative permeability values (both krw and krn) of the proposed model are close to the measured values, and these data sets show a much better linearship.Figure 7
Scatter plot of measured versus predicted values for BC-M, VG-M, and the proposed models.Note that thin-film flow is not considered in this paper, and the film flow is usually distributed on the particle surface or the inner surface of the capillary to form an annular fluid structure [53]. Therefore, the continuous capillary model developed in this study may overestimate the krn of the nonwetting phase and underestimate the krw of the wetting phase, especially in the low relative permeability stage, which can be seen in Figure 7.


### 3.2.4. Parameter Analysis

In this section, the effect of model parametersDf and c on the relative permeability of porous media is also addressed. As seen from Figure 8, the relative permeability of wetting phases decreases as the fractal dimension Df increases, while enhances that of nonwetting phase. This is because a larger Df means a larger percent of small pores, which will decrease the critical radius with a given wetting saturation. Therefore, the transport of wetting phase requires a larger pressure gradient, which reduces its relative permeability and increases the relative permeability of nonwetting phase.Figure 8
The effect of the pore size distribution on relative permeability of porous media.From Figure9(a), the relative permeability of both wetting and nonwetting phases decreases as the nonuniformity factor c increases. Under the assumption of laminar flow, the pressure drop per unit length caused by the flow resistance can be obtained (see Appendix B).
(35)ΔPL′=4μτr21+ϕ2−2/c2−ϕ1/c−12vs.Figure 9
The effect of nonuniformity factorcon the relative permeability of porous media and pressure drop (μ=0.8937×10−3Pas, 25°C, r=1×10−3m, vs=1×10−3m/s, τ=1.5).
(a)
(b)And Figure9(b) shows that pressure drop increases significantly with the increase of c, which means that increasing in the nonuniformity factor c can enhance the fluid resistance and thus lower the transport capacity of porous media.


## 3.2.1. Methods

Fractal-based analytic expressions (Equations (22) and (23)) of the relative permeability contain three parameters, i.e., c, ϕ, and Df. ϕ is the porosity of porous media determined by experiments, c is the nonuniformity factor of pore cross section, and Df is the fractal dimension determined by the best-fitting or the box-counting method [47, 48]. In this section, model parameters are also determined by fitting measured data and calculating the minimum RMSE between prediction values and experimental data.
(32)RMSE=1M∑j=1Mkr,jmeas.−kr,jpred.2.Fractal-based Monte Carlo solutions (Equations (28) and (29)) of the relative permeability contain parameters c, Df, and τi. In this study, Monte Carlo method simulates the randomness of pore sizes based on fractal in porous media. The algorithm for the determination of the relative permeability for unsaturated porous media is summarized as follows:
(1)
Give a porosityϕ and find a nonuniformity factor c of pore cross section
(2)
Find fractal dimensionDf, rmin, rmax and calculate the running total number n of the Monte Carlo simulations according to Equation (24)
(3)
Generate a random numberRi between 0 and 1
(4)
Calculateri by Equation (27) and rmin≤ri≤rmax; otherwise, return to procedure (3)
(5)
Generate random tortuosityτi of 1.25-1.78 suggested by Xu et al. [12]
(6)
Calculate the saturationSew and relative permeability kr (krw and krn)
(7)
CalculateX times kr,x at a same saturation Sew, and then find an average value kr according to kr=1/X∑x=1Xkr,x
(8)
Calibrate model parameters by comparing prediction data with experimental dataAs seen from the above algorithm, based on fractal theory, procedures (1)-(5) reconstruct the pore size distribution in porous media. Evaluation of the model and comparison with other models will be further analyzed and discussed in the next sections.


## 3.2.2. Results and Discussion

Eight data sets are selected for evaluating the proposed model. These data sets involve five fluid systems, i.e., water-air, water-nitrogen, water-steam, water-oil, and oil-gas systems. The pore structure parameters of these porous media are listed in Table2.Table 2
The pore structure parameters of porous media and the corresponding RMSE.
System (wetting-nonwetting)
Porous media
ϕ
Df
c
RMSE
Reference
krn
krw
Water-air
Hygiene sandstone
0.250
1.21
1.98
0.012
0.074
Brooks and Corey [7]
Water-air
Berea sandstone
0.206
1.14
3.73
0.041
0.059
Brooks and Corey [7]
Water-air
Grenoble sand
0.370
1.75
2.18
0.050
0.034
Touma and Vauclin [49]
Water-air
Oso Flaco sand
0.400
1.45
1.90
0.026
0.012
Tuli and Hopmans [9]
Water-nitrogen
Berea core
0.240
1.23
2.40
0.033
0.032
Li and Horne [10]
Water-steam
Core
0.250
1.77
1.35
0.108
0.020
O’Connor [8]
Water-oil
Berea core
0.225
1.18
1.96
0.021
0.110
Kleppe and Morse [50]
Oil-gas
Pyrex core
0.374
1.66
2.18
0.012
0.054
Gates and Leitz [51]Figure4 shows measured data of the water-air system and predicted relative permeability curves by the fractal analysis and Monte Carlo method. These data obtained by the instantaneous profile method in Grenoble sand, Berea sandstone, Hygiene sandstone, and Oso Flaco sand [7, 9, 49]. As shown, the predicted relative permeability (including wetting and nonwetting phases) is in good agreement with the experimental data.Figure 4
Predicted relative permeability for water-air system and the comparison to the experimental data.
(a)
(b)
(c)
(d)For water-nitrogen and water-steam systems, Figure5 shows that both the Monte Carlo solutions and the fractal analytical solutions make a good prediction of the relative permeability, except for the steam phase. This may attribute to the instability of water-steam system in which the immobile steam easily condenses into water so that there is a defect inaccuracy in determining the steam saturation [8].Figure 5
Predicted relative permeability for (a) water-nitrogen system and (b) water-steam system as well as the comparison to the experimental data.
(a)
(b)Studying the relative permeability of water-oil and oil-gas systems is significant for contamination in oil reservoirs as well as in soils. Figure6(a) shows the predicted relative permeability curves of the oil-water system against the corresponding experimental data [50]. The predicted relative permeability of nonwetting phase (oil) is in good agreement with the experimental data while that of wetting phase (water) is not such a good match. The prediction bias may be due to the fact that oil and water are coupled in terms of viscosity or tension, which makes the wetting phase permeability insensitive to changes in saturation. As seen in Figure 6(b), for the oil-gas system, the proposed model can also well determine the relative permeability of both wetting and nonwetting phase in the full range of saturation in Pyrex core [51].Figure 6
Predicted relative permeability for (a) water-oil system and (b) oil-gas system as well as the comparison to the experimental data.
(a)
(b)On the whole, these examples illustrate that the proposed physics-based model, which considers the nonuniformity of pore cross section, can generally capture the transport characteristics for both wetting and nonwetting phase in various porous media. It should be noted that due to the power-law feature of the fractal distribution of pore size, the proposed fractal models in this paper may not be applicable to the porous media of multimodal pore structure, which involves the multipeak distribution of pore size [41].


## 3.2.3. Comparison with Other Models

To further verify the validity of the developed model, the proposed model is compared with the other two models. In the Brook and Corey-Mualem model (BC-M) and the van Genuchten-Mualem model (VG-M) [18, 52], the relative permeability of wetting and nonwetting phases is, respectively, calculated by BC-M:
(33)krw=Sew2.5+2/λ,krn=1−Sew0.51−Sew1+1/λ2,VG-M:(34)krw=Sew0.51−1−Sew1/mm2,krn=1−Sew0.51−Sew1/m2m,where λ is a model parameter related to the pore size distribution of porous media and m is a fitted parameter.Similarly, RMSE is used for the evaluation of these models. Table3 summarizes the RMSE of krw and krn predicted by each model (BC-M, VG-M, and proposed models) for various porous media. For the sake of interpretation, the smallest values of RMSE are highlighted in bold. As seen from Table 3, compared with the BC-M and VG-M models, the proposed model performs best, occupying eleven out of sixteen highlighted values. This overall performance can be reflected by the average values of RMSE (aveRMSE), which are only 0.049 and 0.038 for wetting and nonwetting phases, respectively.Table 3
RMSE values obtained with the BC-M, VG-M, and proposed models.
System (wetting-nonwetting)
Porous media
RMSE (BC-M)
RMSE (VG-M)
RMSE (model)
krn
krw
krn
krw
krn
krw
Water-air
Hygiene sandstone
0.022
0.070
0.020
0.072
0.012
0.074
Water-air
Berea sandstone
0.201
0.096
0.090
0.183
0.041
0.059
Water-air
Grenoble sand
0.112
0.116
0.110
0.079
0.050
0.034
Water-air
Oso Flaco sand
0.082
0.014
0.041
0.030
0.026
0.012
Water-nitrogen
Berea core
0.054
0.090
0.042
0.103
0.033
0.032
Water-steam
Core
0.103
0.122
0.123
0.052
0.108
0.020
Water-oil
Berea core
0.028
0.038
0.016
0.032
0.021
0.110
Oil-gas
Pyrex core
0.019
0.119
0.009
0.128
0.012
0.054
0.078
0.083
0.056
0.085
0.038
0.049Additionally, the differences between these models can be clearly visualized as scatter plots of measured and predicted relative permeability values shown in Figure7. As shown, the BC-M and VG-M models tend to generally overestimate the transport capacity of porous media over the full range of saturation, which may be attributed to the fact that the nonuniformity of pores was overlooked. In contrast, the predicted relative permeability values (both krw and krn) of the proposed model are close to the measured values, and these data sets show a much better linearship.Figure 7
Scatter plot of measured versus predicted values for BC-M, VG-M, and the proposed models.Note that thin-film flow is not considered in this paper, and the film flow is usually distributed on the particle surface or the inner surface of the capillary to form an annular fluid structure [53]. Therefore, the continuous capillary model developed in this study may overestimate the krn of the nonwetting phase and underestimate the krw of the wetting phase, especially in the low relative permeability stage, which can be seen in Figure 7.


## 3.2.4. Parameter Analysis

In this section, the effect of model parametersDf and c on the relative permeability of porous media is also addressed. As seen from Figure 8, the relative permeability of wetting phases decreases as the fractal dimension Df increases, while enhances that of nonwetting phase. This is because a larger Df means a larger percent of small pores, which will decrease the critical radius with a given wetting saturation. Therefore, the transport of wetting phase requires a larger pressure gradient, which reduces its relative permeability and increases the relative permeability of nonwetting phase.Figure 8
The effect of the pore size distribution on relative permeability of porous media.From Figure9(a), the relative permeability of both wetting and nonwetting phases decreases as the nonuniformity factor c increases. Under the assumption of laminar flow, the pressure drop per unit length caused by the flow resistance can be obtained (see Appendix B).
(35)ΔPL′=4μτr21+ϕ2−2/c2−ϕ1/c−12vs.Figure 9
The effect of nonuniformity factorcon the relative permeability of porous media and pressure drop (μ=0.8937×10−3Pas, 25°C, r=1×10−3m, vs=1×10−3m/s, τ=1.5).
(a)
(b)And Figure9(b) shows that pressure drop increases significantly with the increase of c, which means that increasing in the nonuniformity factor c can enhance the fluid resistance and thus lower the transport capacity of porous media.


## 4. Conclusions

In this paper, a new permeability model is established to describe the transport characteristics of porous media with nonuniform pores. The nonuniformity of pores is characterized by capillary tubes with different cross-sectional sizes which cause cross-sectional porosity to be not unique. Furthermore, the nonuniformity factor is introduced to capture the nonunique of porosity and extend Darcy’s law. Based on these and combined with the fractal theory and Monte Carlo technique, the fractal analytical solutions and Monte Carlo solutions of relative permeability are derived, respectively.Sample calculations for eight data sets involving five wetting-nonwetting systems, i.e., water-air, water-steam, water-nitrogen, water-oil, and oil-gas systems, show the rationality of the proposed model in the entire range of wetting saturation. Besides, compared with the BC-M model and the VG-M model, the proposed model has great advantages in predicting the accuracy of relative permeability, especially for nonwetting phase. It has also been found that the nonuniformity of pores can significantly increase the resistance of fluid flow and thus lower the transport capacity of porous media. Besides, the permeability-porosity curve is also derived and it is found that the proposed permeability-porosity model is better than the KC equation in predicting measured data sets.The study in this paper better reflects the general law of fluid transport in porous media with nonuniform pores, which helps to reveal the inherent physical mechanism of heat and mass transfer in porous media.



---

*Source: 5705424-2020-01-31.xml*