Relative Permeability of Porous Media with Nonuniform Pores

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, waternitrogen, 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.


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][4][5][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][8][9][10].
Predictive models describing relative permeabilitysaturation relationships can be divided into empirical models, statistical models, and fractal models [11][12][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][17][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][23][24][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][27][28][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][32][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 inte-grated into the extended Darcy's law. Based on fractal theory and Monte Carlo technology, the relative permeabilitysaturation 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.

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.  II-II   II-II   II-II   I-I  I-I [36]: where χ = ð1 + a/b 4 Þis a resistance coefficient related to the nonuniformity of pore cross section, where a ða > 0Þ and b ð0 < b ≤ 1Þ are the length factor and size factor of the capillary-throat, respectively, and μ is the fluid viscosity. τ = L t / L 0 ðτ ≥ 1Þ is the pore tortuosity, where L t and L 0 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: where n is the total number of pores in the REV, ΔP is the pressure gradient, and subscript i indicates the i th capillary tube. As seen from Figure 1, 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 ϕ = ϕ r c , 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 c should 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: where k s is the permeability of the saturated REV, A p is the average pore area of the REV cross section, and thus, 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).

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 r c . Therefore, the volumetric flow rate of the wetting and nonwetting phases in unsaturated porous media can be determined as follows: where ΔP w and ΔP nw 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]: where k w is the permeability of wetting phase in unsaturated porous media. Substituting Equation (5) into the rearranged Equation (7) yields the following: 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 as It should be noted that in this study, S ew 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 k r = k w /k s .

Geofluids
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 of r c − r max .
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.

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 scale r can be expressed as follows [21,22]: 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 D f 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: where dN > 0 and the negative sign (-) implies that the number of capillary tubes decreases as the capillary radius increases.

Permeability-Porosity Curve.
In order to develop a specific permeability model, the number of pores with a size r determined by Equation (12) can be substituted into Equation (4) as follows: 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) Here, ln ϕ r = ð2 − D f Þ ln ðr min /r max Þ [22] can be introduced into Equation (14) yields the following equation: Therefore, the permeability-porosity curve can be determined by substituting Equation (15) into Equation (13), that is, where Note that in Equation (16), ϕ r = ϕ 1/c and ðr min /r max Þ 4−D f ≈ 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: where C KC 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.

Relative Permeability-Saturation Curve.
Similarly, in order to develop a specific fractal model of relative permeability, the number of pores with a size r computed by Equation (12) can be substituted into Equation (9).
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: where And substituting Equation (21) into Equation (20), it can obtain the relative permeability k rw of the wetting phase expressed by the saturation S ew .
where (22) to simplify the model. Similarly, the relative permeability of the nonwetting phase can be given as follows: 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 D f .

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).
Dividing Equation (12) by Equation (24) gives where f ðrÞ = D f r min D f r −ðD f +1Þ . Yu et al. [23] argued that if the pore radius of the porous media satisfies ðr min /r max Þ D f ≈ 0, f ðrÞ can be regarded as a probability density function of pore size. Fortunately, for natural porous media (e.g., soil), r min / r max ≪10 −2 so that ðr min /r max Þ Df = 0 holds approximately. Therefore, the cumulative probability (R), from the smallest radius r min to the radius r, can be written as follows: Equation (20) describes the random probability of the pore size distribution in the range of r min − r max , and when r ⟶ r min , R ⟶ 0 as well as r ⟶ r max , R ⟶ 1.
On the other hand, if R is regarded as a random number from 0 to 1, the pore radius of porous media can be expressed as follows: where r min < r < r max . 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. Substituting Equation (27) into Equations (9) and (10), it can obtain the Monte Carlo solutions of the relative permeability for wetting and nonwetting phases. where and R is a computer-generated random number in the range of 0-1, and the subscript i indicates the i th 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 r c . r min and r max are the minimum and 5 Geofluids maximum pore sizes of the REV, respectively. τ is the tortuosity of pores in the REV and the subscript i indicates the i th 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.

Applications
3.1. Permeability-Porosity Curve. Fractal-based permeabilityporosity 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).
ln k meas: where M is the number of the test data, k r,j meas: is the j th measured data point, and k r,j pred: 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 pre-dicts good the measured data over a wide range of magnitude (about 4 to 11 orders of magnitude). (22) and (23)) of the relative permeability contain three parameters, i.e., c, ϕ, and D f . ϕ is the porosity of porous media determined by experiments, c is the nonuniformity factor of pore cross section, and D f 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.

Methods. Fractal-based analytic expressions (Equations
Fractal-based Monte Carlo solutions (Equations (28) and (29)) of the relative permeability contain parameters c, D f , 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 dimension D f , r min , r max and calculate the running total number n of the Monte Carlo simulations according to Equation (24) (3) Generate a random number R i between 0 and 1 (4) Calculate r i by Equation (27) and r min ≤ r i ≤ r max ; otherwise, return to procedure (3) (5) Generate random tortuosity τ i of 1.25-1.78 suggested by Xu et al. [12] (6) Calculate the saturation S ew and relative permeability k r (k rw and k rn ) (7) Calculate X times k r,x at a same saturation S ew , and then find an average value k r according to k r = ð1/XÞ ∑ X x=1 k r,x (8) Calibrate model parameters by comparing prediction data with experimental data As 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.

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 Table 2.  6 Geofluids Figure 4 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.
For water-nitrogen and water-steam systems, Figure 5 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].
Studying the relative permeability of water-oil and oil-gas systems is significant for contamination in oil reservoirs as well as in soils. Figure 6(a) shows the predicted relative permeability curves of the oil-water system against the corresponding experimental data [50]. The predicted relative    7 Geofluids 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].
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].  22) and (23)) Fractal analytical (Eqs. (16) and (17) 22) and (23)) Fractal analytical (Eqs. (16) and (17) Meas. Monte Carlo solutions (Eqs. (22) and (23)) Fractal analytical (Eqs. (16) and (17)) Hygiene sandstone  22) and (23)) Fractal analytical (Eqs. (16) and (17)) (d) Figure 4: Predicted relative permeability for water-air system and the comparison to the experimental data.
8 Geofluids VG-M: 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. Table 3 summarizes the RMSE of k rw and k rn 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.
Additionally, the differences between these models can be clearly visualized as scatter plots of measured and predicted relative permeability values shown in Figure 7. 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 Meas. Monte Carlo solutions (Eqs. (22) and (23)) Fractal analytical (Eqs. (16) and (17) Meas. Monte Carlo solutions (Eqs. (22) and (23)) Fractal analytical (Eqs. (16) and (17)) (b) Figure 6: Predicted relative permeability for (a) water-oil system and (b) oil-gas system as well as the comparison to the experimental data. Monte Carlo solutions (Eqs. (22) and (23)) Fractal analytical (Eqs. (16) and (17) Monte Carlo solutions (Eqs. (22) and (23)) Fractal analytical (Eqs. (16) and (17)) (b) Figure 5: Predicted relative permeability for (a) water-nitrogen system and (b) water-steam system as well as the comparison to the experimental data.
9 Geofluids nonuniformity of pores was overlooked. In contrast, the predicted relative permeability values (both k rw and k rn ) of the proposed model are close to the measured values, and these data sets show a much better linearship.
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 k rn of the nonwetting phase and underestimate the k rw of the wetting phase, especially in the low relative permeability stage, which can be seen in Figure 7.

Parameter Analysis.
In this section, the effect of model parameters D f 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 D f increases, while enhances that of nonwetting phase. This is because a larger D f 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.
From Figure 9(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).  10 Geofluids And Figure 9(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.

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 permeabilityporosity 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. The reference radius  Therefore, the pressure drop caused by viscous energy loss can be written as follows: where L′ = ð1 + aÞL 0 and χ′ = χ/ð1 + aÞ = ð1 + a/b 4 Þ/ð1 + aÞ, which is defined as the resistance coefficient.
A.2 Local (Kinetic) Energy Loss. In order to determine the local energy loss during the fluid through the cross section of sudden expansion or contraction, we consider the capillary as shown in Figure 10 and assume that the fluid flows from the A 1 section to the A 3 section. The loss of pressure head caused by the change of cross section can be expressed as follows [54]: where v is the real average velocity of the cross sections A 1 or A 3 . ξ is the coefficient of local energy loss. For a sudden contracting pore of cross section A 2 , the coefficient of local energy loss can be written as follows: In this case, the loss of head can be determined by On the other hand, for a sudden expending pore of cross section A 3 , the coefficient of local energy loss can be written as follows: Hence, the loss of head during the fluid flows through the expending pore can be calculated as follows: : ðA:7Þ Combining Equations (A.5) and (A.7) as well as ΔP 2 = ρgh l , the pressure drop caused by local energy loss can be written as follows: ðA:8Þ Note that v s = v/ϕ has been introduced into Equation (A.8). Therefore, we can obtain the total pressure drop per unit length during the fluid flows through the capillary tube (see Figure 10) by Equations (A.2) and (A.8).

ΔP
ðA:9Þ As seen from Equation (A.9), the energy loss increases as the fluid flow rate increases. However, for soils and rocks, the superficial velocity v s of fluids is typically much less than 10 -2 m/s, in which case the local energy loss caused by the change in pore cross-sectional size can be neglected compared to the viscous energy loss.

B. Pressure Drop-Nonuniformity Factor Curve
As shown in Figure 1, the porosity of the reference cross section can be calculated as follows: Similarly, the porosity of the cross section with a radius br can be calculated as follows: π br i ð Þ 2 = b 2 ϕ r : ðB:2Þ Therefore, the porosity of the REV in Figure 1 can be written as follows: It can be found from Equation (B.4), when b ⟶ 1, c equals 1, which means that the size of pore section is unique. If the parameter a is set to a constant (e.g., let a = 1), combined with Equations (A.9) and (B.4) and neglecting the local

Data Availability
Previously reported relative permeability data were used to support this study. These prior data sets are cited at relevant places within the text as References [7-10, 45, 46, 49-51].