Permeability Prediction of Saturated Geomaterials with Revised Pore–Solid Fractal Model and Critical Path Analysis

: The revised pore–solid fractal (PSF) model is presented by using the largest inscribed circle-based geometries of squares or cubes to replace the original pore or solid subregions as the new pore or solid phase in porous media. The revised PSF model changes the discrete lacunar pore and solid phases in the original PSF model to integrated. Permeability is an intrinsic property of geomaterials and has broad applications in exploring ﬂuid ﬂow and species transport. Based on the revised PSF model and critical path analysis, a fractal model for predicting the permeability of saturated geomaterials is proposed. The permeability prediction model is veriﬁed by comparison with the existing predicted model and the laboratory testing. The results show that the predicted permeabilities match the measured values very well. This work provides a theoretical framework for the revised PSF model and its application in predicting the permeability of geomaterials.


Introduction
Accurate prediction of the permeability for porous media such as geomaterials is an interesting topic in hydrology, geophysics, petroleum engineering, geotechnical engineering, etc. Permeability is the reflection and measurement of the complex relationship between the geometry and the topology of the pore phase in porous media, and it has been commonly used in exploring fluid flow and species transportation in porous media. The literature on predicting permeability is extensive. In terms of experimental testing, mercury intrusion porosimetry [1,2] and nuclear magnetic resonance [3,4] have been widely used in quantifying the pore properties of a porous medium and have become the effective verification methods for the theoretical models of estimating permeability. In terms of calculation methods, empirical models [5][6][7][8] and theoretical models considering porosity and pore or solid particle size distribution have been proposed; examples include percolation theory [9][10][11][12], effective medium approximation [11][12][13], critical path analysis [14,15] and fractal geometry [16]. Critical path analysis (CPA) can explore both throat connectivity and tortuosity of pore channels in porous media, which has been used to estimate the permeability of heterogeneous geomaterials [17][18][19][20].
The micro geometric structures of porous media have self-similar properties that are better interpreted and described with fractal geometry. A porous medium can be regarded as a fractal body, and its micro geometric structures can be analyzed by fractal models such as a solid or pore fractal model (Sierpinski carpet or Menger sponge) [21], pore-solid fractal model [22], bunched capillaries fractal model [16,23,24], multifractal model [25] and intermingled fractal units model [26]. The bunched capillaries fractal model [16,23,24] is mostly applied to estimate the permeability of a porous medium. This model considers the porous medium to consist of a bunch of cylindrical pore capillaries of different diameters; the cumulative pore size distribution and the length of capillaries follow the fractal scaling law. However, the permeability prediction models based on bunched capillaries fractal Fractal Fract. 2022, 6, 351 2 of 11 model lack a characterization of the pore interconnections in real existence in a porous medium, which is an idealization of a natural porous medium.
The pore or solid fractal model based on the standard two-dimensional Sierpinski carpet or three-dimensional Menger sponge has been mainly applied to focus on the cumulative pore or solid particle size distribution, the interface between the pore and solid, the mass (volume) and the water retention curve. However, based on the monofractal structure, these simplified fractal models cannot reflect the real porous media, especially geomaterials.
A geometrical, multiscale model of soil structure, which is named a pore-solid fractal (PSF) model [22], characterizes the structure of soil exhibits self-similarity to a degree; that is, where local structure appears, it is regarded as similar to the whole structure. It is achieved through the incorporation of a pore phase, solid phase and fractal phase that can be decomposed into smaller pore, solid and fractal subregions when iteration increases. The pore phase and the solid phase both have volume fractions greater than 0 and less than 1 even for infinite iterations, and they coexist at each level in the hierarchy, which is a main difference from the classical pore or solid fractal model with only the pore and solid phases based on the standard Sierpinski carpet or Menger sponge. In addition, when the fractal phase is considered as a solid phase or pore phase, the PSF model is simplified to a solid fractal model or a pore fractal model, respectively. At present, the PSF model is mostly used to estimate the soil water retention function [27,28] and determine the unsaturated hydraulic conductivity [29].
In this article, a fractal model for predicting the permeability of saturated geomaterials is proposed. Section 2 establishes a PSF model with circular-based pores to reflect the pore phase of different geometric pores and a PSF model with circular-based solid particles to describe the solid phase of different geometric solid particles. The permeability for the fragmented pore or solid phase is derived using the Kozeny-Carman equation, and the permeability for continuous geomaterials is predicted by combining the PSF model with circular-based pores and critical path analysis in Section 3. The revised PSF models degenerate to the general solid or pore fractal model, and the predicted model is verified by comparison with the existing predicted model and the laboratory testing in Section 4. The conclusion is expressed in Section 5.

Revised PSF Model
The original PSF model uses a square (or cube) as a basic unit that is the same as the pore or solid fractal model, which makes the contact between the pore and the solid complicated and makes theoretical analysis and numerical simulation difficult. In order to improve the applicability of the PSF model, it should be revised.

PSF Model with Circular-Based Pores
The model construction is consistent with the original PSF model. Firstly, divide a region of size L into N equal subregions in the Euclidean space of dimension d as the initiator. N is calculated by where r is the linear length of each subregion and n is the number of divisions in each direction. Define Nz subregions as the self-similar replicated phase, Ns subregions as the solid phase of porous media and Np subregions as the pore phase.
where p, s and z are the proportions of the pore, solid and self-similar replicated phases in the initiator, respectively. Take the largest inscribed circles (cylinders) in the Np subregions as the new pore phase and consider the remaining parts of the Np subregions as the solid phase. The above operations are performed in 1st generation as a generator. In the next ith (i ≥ 2) generation, a recursive process replaces (Nz) i−1 subregions with the reduced replicate of the generator with the ratio (1/n) i−1 . An example of the construction process of the PSF model with circular-based pores is illustrated in Figure 1. where p , s and z are the proportions of the pore, solid and self-similar replicated phases in the initiator, respectively.
Take the largest inscribed circles (cylinders) in the Np subregions as the new pore phase and consider the remaining parts of the Np subregions as the solid phase. The above operations are performed in 1st generation as a generator. In the next ith (i 2 generation, a recursive process replaces ( ) where D is the fractal dimension of the original PSF model; ( ) is the cumulative number of pores after ith generation; p i r is the pore diameter after ith generation, An inscribed circle (cylinder) of maximum diameter in a square (cube) subregion of size i r instead of the square (cube) pore subregion leads to a decrease in porosity as where α is geometric parameter, circle (cylinder) = / 4 α π , semicircle (semicylinder) At the end of the ith generation, the decreased porosity is Therefore, the porosity i φ of the PSF model with circular-based pores is shown as The cumulative pore size distribution P ≥ r p i of the PSF model with circular-based pores is still calculated using the corresponding equation of the original PSF model [22]: where D is the fractal dimension of the original PSF model; P ≥ r p i is the cumulative number of pores after ith generation; r p i is the pore diameter after ith generation, r p i = L n i . An inscribed circle (cylinder) of maximum diameter in a square (cube) subregion of size r i instead of the square (cube) pore subregion leads to a decrease in porosity as where α is geometric parameter, circle (cylinder) α = π/4, semicircle (semicylinder) α = π/8 and square (cube) α = 1. At the end of the ith generation, the decreased porosity is Therefore, the porosity φ i of the PSF model with circular-based pores is shown as When i → ∞ , r p i → 0 , and then φ = pα p+s . The surface area of circular-based pores S p (i) is obtained by Noting N = n d , Equation (6) changes to Since the construction process of the PSF model with circular-based pores is consistent with the original PSF, refer to the original PSF model (Perrier et al., 1999 [22]); the pore-solid interface area S sp (i) is achieved by So, the surface area of the solid phase S s (i) is written as Assume M s ≤ r p i is the cumulative mass of solid phase when the pore diameter is less than or equal to r p i , which is obtained by multiplying the solid volume of pore diameter smaller than or equal to r p i and solid density ρ s . After T − 1 iterations (Tth generation) of the PSF model with circular-based pores, M s ≤ r p i is derived as follows: When T → ∞ , r p T → 0 , and then

PSF Model with Circular-Based Solid Particles
The construction process is similar to the PSF model of circular-based pores except that the largest inscribed circles (spheres) in the Ns(Nz) i−1 subregions of size r i are used to represent the new solid phase and the remaining parts of the Ns(Nz) i−1 subregions of size r i are regarded as the pore phase.
The cumulative solid particle size distribution S ≥ r s i of the PSF model with circularbased solid particles is still calculated as where S ≥ r s i is the cumulative number of solid particles after ith generation; r s i is the solid particle diameter after ith generation, r s i = L n i . An inscribed circle (sphere) of maximum diameter in a square (cube) subregion of size r i replaces the square (cube) solid subregion, which causes an increase in porosity as where β is a geometric parameter. When the subregion is a circle, β = π/4; sphere, β = π/6; square (cube), β = 1. At the end of the ith generation, the increased porosity is Therefore, the porosity φ i of the PSF model with circular-based solid particles is obtained as when i → ∞ , r s i → 0 , and then φ = p+sυ p+s . According to the PSF model with circular-based pores, the surface area of circularbased solid particles S s (i) is written as The surface area of the pore phase S p (i) is obtained as The cumulative mass of the solid phase M s ≤ r s i when the solid particle diameter is less than or equal to r s i is derived as When T → ∞ , then When the subregions of solid phase are squares (cubes), β = 1. Equation (20) coincides with the cumulative mass of the solid phase given by Bird et al. [30].

Fragmented Pore or Solid Phase
The Kozeny-Carman equation is widely used to estimate the permeability of porous media [5,6], which is expressed as where k is the permeability, φ is the porosity, C is the empirical parameter and S is the specific surface area per unit volume of the solid phase.
From the revised PSF model construction, a porous medium can be fragmented into circular-based pores of different sizes embedded in an aggregated solid phase or circularbased solid particles surrounded by pore space. For circular-based pores, according to Equation (11), the specific surface area per unit volume of solid phase in ith iteration S i is If nz = 1, If nz = 1, For circular-based solid particles, based on Equation (18), S i is obtained by If nz = 1, If nz = 1,

Continuous Geomaterials
Equations (4) and (13) show that the cumulative discrete pore size distribution and solid particle size distribution of the revised PSF model obey P ∼ r −D and S ∼ r −D , respectively. Hunt argues that pore size is commonly used to characterize continuous porous media and can be translated into the case of using solid particles [31,32]. Learning from this, the continuous probability density function of pore size f (r) is obtained as where r max is the maximum pore diameter. When α = 1, Equation (29) degenerates to the formula proposed by Ghanbarian et al. [11,29].
The cumulative probability distribution function for pore diameter smaller than r is achieved by integrating Equation (29) as If r = r max , It is obvious that Equation (31) is the equivalent of Equation (7). The water retention function θ is obtained with the process of the original PSF model [30] as where h min and h max are the water tension heads at which the maximum and minimum pore diameters occur in the porous media, respectively. According to critical path analysis [17,19,20], the percolation threshold p c of pore volume for saturated condition (θ = φ) is defined as where r c is the critical pore diameter for saturated porous media. Transform Equation (33) to obtain The percolation threshold p c and the critical pore diameter r c for an unsaturated condition (θ = φ) are ( and (36) The generalized equation for permeability k derived by CPA is expressed as [17,18,20,29,33] where F is the formation factor and C is a constant according to the specific porous materials. The formation factor F of saturated porous media is approximated from the water retention curve for the saturated condition (θ = φ) and is written as [20,33] (38) Equation (38) has been successfully used for predicting the permeability of geomaterials [20,33].
Substitute the critical pore diameter Equation (34) and the formation factor Equation (38) into Equation (37) to yield the permeability as (39)

General Pore or Solid Fractal Model
When s = 0, p + z = 1, and the PSF model with circular-based pores degenerates to the general solid fractal model. P and z are the proportions of pore and solid. An example of a general solid fractal model with circular-based pores is illustrated in Figure 2.

General Pore or Solid Fractal Model
When 0 s = , 1 p z + = , and the PSF model with circular-based pores degenerates to the general solid fractal model. P and z are the proportions of pore and solid. An example of a general solid fractal model with circular-based pores is illustrated in Figure  2. The porosity of the solid fractal model with circular-based pores is where ' D is the fractal dimension of the general solid fractal model.
For the solid fractal model based on a standard two-dimensional Sierpinski carpet or three-dimensional Menger sponge, similarly, take the largest inscribed circles (cylinders) in the square (cube) pore subregions of size i r as the new pore phase; the porosity is ( ) ( ) It can clearly be seen that the solid fractal model with circular-based pores based on standard Sierpinski carpet or Menger sponge is a special case of the general solid fractal model changed from PSF. When 0 p = , 1 s z + = , and the PSF model with circular-based solid particles changes to the general pore fractal model. s and z are the proportions of solid and pore. The porosity of the pore fractal model with circular-based solid particles is where '' D is the fractal dimension of the general pore fractal model. The porosity of the solid fractal model with circular-based pores is where D is the fractal dimension of the general solid fractal model.
For the solid fractal model based on a standard two-dimensional Sierpinski carpet or three-dimensional Menger sponge, similarly, take the largest inscribed circles (cylinders) in the square (cube) pore subregions of size r i as the new pore phase; the porosity is It can clearly be seen that the solid fractal model with circular-based pores based on standard Sierpinski carpet or Menger sponge is a special case of the general solid fractal model changed from PSF.
When p = 0, s + z = 1, and the PSF model with circular-based solid particles changes to the general pore fractal model. s and z are the proportions of solid and pore. The porosity of the pore fractal model with circular-based solid particles is where D is the fractal dimension of the general pore fractal model.
For the pore fractal model based on a Sierpinski carpet or Menger sponge, the largest inscribed circles (spheres) in the square (cube) solid subregions of size r i are used to represent the new solid phase; the porosity is Equation (43) coincides with the Khabbazi formula [34]. Obviously, the pore fractal model with circular-based solid particles based on a Sierpinski carpet or Menger sponge is also a special case of the general pore fractal model transformed from PSF.

Permeability
Equation (39) is applied to predict the permeabilities of the geomaterials samples used in Daigle's prediction [17]. The comparison of predicted permeabilities and the measured values is shown in Figure 3. .
Equation (43) coincides with the Khabbazi formula [34]. Obviously, the pore fractal model with circular-based solid particles based on a Sierpinski carpet or Menger sponge is also a special case of the general pore fractal model transformed from PSF.

Permeability
Equation (39) is applied to predict the permeabilities of the geomaterials samples used in Daigle's prediction [17]. The comparison of predicted permeabilities and the measured values is shown in Figure 3.  It can be seen from Figure 3 that the predicted permeabilities with Equation (36) are almost consistent with the measured permeabilities and k mea − k pre = ±0.4 × 10 −11 is the upper and lower boundary. Meanwhile, the ratio k pre k mea is between 0.9 and 1.1 when log r max r min ≥ 2.7; for log r max r min < 2.7, the error between the predicted permeabilities and the measured values is 20% to 10%.
The difference between the predicted values in this article and Daigle's predictions is mainly due to the difference in the calculation method of formation factor F. Therefore, the accurate expression of formation factor F determines the accuracy of permeability prediction. Figure 4 shows the predicted permeability k pre is highly correlated with the critical pore diameter r c . The permeability is not only influenced by micro geometric characteristics such as pore diameter, pore size distribution and interface area of pores and solid particles, but also influenced by topological properties such as connectivity between pore and solid phases in the porous medium, which were also obtained by Ghanbarian et al. [20]. Figure 4 shows the predicted permeability pre k is highly correlated with the critical pore diameter c r . The permeability is not only influenced by micro geometric characteristics such as pore diameter, pore size distribution and interface area of pores and solid particles, but also influenced by topological properties such as connectivity between pore and solid phases in the porous medium, which were also obtained by Ghanbarian et al. [20].

Conclusions
The revised PSF models use the largest inscribed circle-based geometries in pore or solid subregions instead of original squares or cubes to represent pores or solid particles in porous media, which increases the integrated contact between pores and the solid phase.
A fractal model for predicting the permeability of geomaterials based on the PSF model with circular-based pores and critical path analysis is proposed; the permeability is highly correlated with the critical pore diameter.

Conclusions
The revised PSF models use the largest inscribed circle-based geometries in pore or solid subregions instead of original squares or cubes to represent pores or solid particles in porous media, which increases the integrated contact between pores and the solid phase.
A fractal model for predicting the permeability of geomaterials based on the PSF model with circular-based pores and critical path analysis is proposed; the permeability is highly correlated with the critical pore diameter.