# A New Fractal Permeability Model for Porous Media Based on Rough Capillary Channels

Authors: Yibo Zhao
Journal: Advances in Materials Science and Engineering
 (2022)

---

## Abstract

Porous media are assumed as a bunch of curved capillaries of rough pore-solid surface, with capillary size distributions and surface roughness following the fractal scaling laws, for which a permeability model is derived to capture both tortuosity and roughness of the pore space in this work. First, the fractal geometry theory and the regarding methods are used to simulate tortuosity and roughness, and then, the permeability of porous media is linked to pore area fractal dimension, tortuosity fractal dimension, relative mean roughness, and other structural parameters (e.g., characteristic length and maximum and minimum pore size). Each parameter in the proposed model has specific physical meaning, which is able to reveal certain mechanisms that affect permeability comprehensively. For several porous media samples, the predicted permeability data based on the current fractal model are compared with the experimental measurement data and the permeability model predictions for other porous media with smooth capillary channels. A good agreement was found between the predicted values of the new permeability fractal model and the experimental data.


---

## Body

## 1. Introduction

Solute transport in porous media is one of the most common phenomena, which is often represented by electrical conductivity, hydraulic conductivity, or diffusivity. The prediction of the ability for any fluid to flow through porous media, that is, (absolute) permeability, is of great importance to many research fields such as petroleum engineering [1], nuclear waste disposal [2], and geothermal exploration [3, 4]. The well-known Kozeny-Carman (KC) equation [5, 6], a conventional means of modelling permeability, is an analytical nonlinear function of porosity and surficial area, which is widely applied to many areas and frequently modified by coupling with additional features for different interests, aiming at improving the accuracy in the permeability estimation [7–13]. Yu and Cheng pioneered a new method for correcting the K-C equation based on the fractal geometry theory [14]. The biggest advantage of this original fractal model is that there is no empirical parameter, which needs to be obtained by fitting out through regional data. However, the sizes of particles and clusters are included in the permeability model, which is inconvenient to measure. To solve this problem, Xu and Yu [15] derived the analytical expression of the permeability of homogeneous porous media based on the capillary channel model of porous media, which is expressed as a function of fractal dimension, porosity, and maximum porosity. However, the KC constant in this model is derived from the assumption of a square geometric model, so it is not representative. Du et al. [16] developed a fractal permeability prediction model for water flow through hydrate-bearing porous media, which incorporated the pore fractal dimension, tortuosity fractal dimension, porosity, hydrate saturation, and maximum pore size. However, some parameters in the model (such as tortuosity fractal dimension, the water phase fractal dimension, and tortuosity) are calculated using approximate models, and multiple parameters are related to each other, which will lead to large errors in the results for specific rock samples. Xiao et al. [17] proposed a new fractal permeability model to study the permeabilities and KC constants of porous media consisting of solid particles and porous fibers. The model incorporated the pore fractal dimension, tortuosity fractal dimension, porosity, hydrate saturation, and maximum pore size.The above results on the permeability were all obtained under the assumption of smooth inner wall surfaces of capillaries or smooth surfaces of the pores within porous media. The influence of the roughness of the inner walls of capillaries has not been considered. In fact, the internal surfaces of natural porous media are rough, which significantly influences the permeability. Many researchers have also used the fractal geometry theory to characterize the roughness of pore-solid surfaces. Li et al. [18] established a fractal dimension model of rough surface by considering rough elements as cone shapes, which provided a model basis for studying the fractal characteristics of rough elements in porous media. Yang et al. [19] studied the fractal model of roughness using conical rough elements too and deduced the pressure gradient, friction coefficient, and Poiseuille number of laminar flow in the rough surfaces of microchannels. However, their fractal models of rough surfaces contain too many structural parameters that are extremely difficult to determine in practice, so their models are commonly considered useful in the theory rather than applicable in applications like predicting material properties or analyzing flow behaviors.In this work, a permeability model with “relative mean roughness” is established based on the fractal theory, and the tortuosity of the capillary channels in the model is also reconstructed, improving the accuracy of the results. Finally, the prediction results of the newly established permeability model are compared with the experimental values and the prediction results of other models without roughness factors to demonstrate the effectiveness of the new permeability model.


## 2. Permeability Model of Porous Media with Rough Capillary Walls

### 2.1. Fractal Geometry Theory of Porous Media

The cumulative distributions of pore sizes follow the power law [20], and the number of pores that have diameters greater than or equal to λ can be expressed as [21](1)NL≥λ=λmaxλDf,where N is the number of pores, L is the length scale, and λ is pore diameter; λmin≤λ≤λmax, where λmin and λmax are the minimum and maximum diameters in a porous medium, respectively, and Df is the pore area fractal dimension. Apparently, 1<Df<3 in two dimensions and 2<Df<3 in three dimensions. With (1), the total number of pores between λmin and λmax in a porous medium can be obtained by(2)NtL≥λmin=λmaxλminDf.From the derivative of (1) with respect to λ, we have the number of pores that have diameters ranging between λ and λ+dλ:(3)−dN=DfλmaxDfλ−Df+1dλ.Combining (2) and (3), we can obtain the following expression:(4)−dNNt=DfλminDfλ−Df+1dλ=fλdλ,where fλ=DfλminDfλ−Df+1 is the probability density function associated with pore size distribution. According to the property of the probability density function ∫λminλmaxfλdλ=1, it is clear that λmin/λmaxDf=0 for fractal porous media, and the porous media in nature are proved to satisfy λmin/λmax<10−2 [14], in principle. Furthermore, Yu and Li [21] established a relationship between porosity and pore area fractal dimension based on fractal geometry theory, that is,(5)ϕ=λminλmaxdE−Df,where dE is the Euclidean dimension, dE = 2, 3 in two dimensions and three dimensions, respectively.From (5), Df can be rewritten as(6)Df=dE−lnϕlnλmin−lnλmax.Considering pore tortuosity, the average valueτ¯ in a porous medium is defined as(7)τ¯=LtL0,where Lt and L0 are the actual length and the straight length of a flow path, respectively. The fractal dimension DT of tortuosity indicates the bending extent of streamline paths. DT = 1 means that the flow path is straight; otherwise, it is tortuous; it is clear that the greater the tortuosity, the greater the DT is.The tortuosity fractal dimension can then be described by [22](8)DT=1+lnτ¯lnL0−lnλ¯,where λ¯ is the expectation of pore sizes, that is,(9)λ¯=∫λminλmaxλfλdλ=DfDf−1λmin1−λminλmaxDf−1≈DfλminDf−1.Equations (1)–(9) are from the theoretical base of the present work.


### 2.2. Fractal Characterization of Surface Roughness

Conical roughness elements are studied herein, as shown in Figure1 [23]. Apparently, the fluid flow in rough capillaries is affected by the roughness that is modelled by randomly distributed elements, causing the flow path to bend to a certain extent, like the case of fluid flowing through natural porous media.Figure 1
(a) Capillary bundle model for simulating porous media, whereLt denotes the actual length of capillaries and L0 denotes the linear length of the streamline, which is also the side length of the elementary volume. (b) The roughness elements on the pore wall surface.Due to the tortuous feature of capillary tubes, the actual capillary lengthLtλ is greater than or equal to the straight length L0, and Ltλ fits the power law [14], that is,(10)Ltλ=L0DTλ1−DT.For a rough capillary with a lengthLtλ, according to (3) and (10), the total surface area in an elementary can be obtained as(11)S0=∫λminλmaxπλLtλ−dN=πL0DTDfλmax2−DT2−DT−Df1−λminλmax2−DT−Df.Because of the self-similarity of pore channels, in an elementary volume of a porous medium, the volumeVr of all conical roughness elements is(12)Vr=L031−ϕ.Furthermore, according to (11) and (12), the mean height of the roughness elements can be effectively expressed as(13)h¯=VrS0=L03−DT1−ϕ2−DT−DfπDfλmax2−DT1−λmin/λmax2−DT−Df.The relative mean roughnessξ¯ in a rough cylindrical capillary can then be defined as(14)ξ¯=2h¯λ¯=2L03−DT2−DT−DfDf−1πDf2λmax2−DTλmin1−λmin/λmax2−DT−Df1−ϕ,where λ¯ is given by (9).According to (14), the relative mean roughness of a porous medium can be expressed as a function of porosity, pore area fractal dimension, tortuosity fractal dimension, characteristic length, and pore size (e.g., maximum and minimum pore size). Equation (14) reveals some of the mechanisms related to roughness in rough channels/capillaries.


### 2.3. Permeability Model for Porous Media

Modifying the well-known Hagen–Poiseuille equation, Yu et al. [14] derived an expression for the flow rate through a single curved smooth cylindrical capillary:(15)qλ=π128ΔPLtλλ4μ,where μ is fluid viscosity, ΔP is pressure drop, and λ is pore diameter.Due to the roughness, the effective diameters for fluid flow in capillaries certainly are smaller than those in smooth capillaries. Hence, (15) can be refined as(16)qrλ=π128ΔPLtλλ−2h¯4μ,where h¯ is the effective height of conical roughness elements on average, as described by (13). (16) suggests that due to surface roughness in capillary, the flow rate decreases if the pressure difference maintains constant, or this pressure difference increases if the flow rate is kept constant.Therefore, the total flow in an elementary volume can be expressed as(17)Qr=∫λminλmaxqrλ−dN=π⋅ΔP⋅Df⋅λmaxDf128L0DTμ∫λminλmaxλ−2h¯4λDT−Df−2dλ,where(18)∫λminλmaxλ−2h¯4λDT−Df−2dλ=λmax3+DT−Df3+DT−Df1−λminλmax3+DT−Df−4ξ¯λ¯λmax2+DT−Df2+DT−Df1−λminλmax2+DT−Df+6ξ¯2λ¯2λmax1+DT−Df1+DT−Df1−λminλmax1+DT−Df−4ξ¯3λ¯3λmaxDT−DfDT−Df1−λminλmaxDT−Df+ξ¯4λ¯4λmaxDT−Df−1DT−Df−11−λminλmaxDT−Df−1.In fractal porous media model,1<DT<3,1<Df<3 so −4<2−DT<Df<0. λmin/λmax<0.01, and thus, according to (14) and considering the pore size of the actual porous medium, we get ξ¯2≅0, ξ¯3≅0, ξ¯4≅0.Furthermore, as3+DT−Df>1,2+DT−Df>0, we have λmin/λmax3+DT−Df≅0, λmin/λmax2+DT−Df≅0. Considering the roughness of the capillary wall, the total flow in an elementary volume can be simplified as(19)Qr=πΔPDfλmax2+DT128L0DTμλmax3+DT−Df−4ξ¯λ¯2+DT−Df.According to Darcy’s law, the absolute permeabilityK of a porous medium can be given by [14](20)K=μL0QrΔPA,where A is the cross-sectional area of an elementary volume (cube) with a side length of L0.By substituting (19) into (20) with A=L02, the permeability K∗ of the porous medium with rough pore channels can be obtained as(21)K∗=πDfλmax3+DT128L01+DT3+DT−Df1−Cξ¯,where C is a constant, which is given as(22)C=3+DT−Df2+DT−Df⋅4DfDf−1⋅λminλmax.(21) shows that the permeability of a porous medium with rough capillaries is jointly determined by pore fractal dimension Df, tortuosity fractal dimension DT, structural parameters (L0, λmin, and λmax), and the relative mean roughness ξ¯. The structural parameters are apparently not independent of each other; instead, there is a certain correlation between them. (21) also demonstrates that the larger the relative mean roughness is, the smaller the permeability becomes. This is because the flow resistance increases with the increase in roughness, which is consistent with our understanding about multiphase flow in porous media. The permeability is proportional to the fourth power (3+DT>4) of the maximum pore diameter λmax, so the maximum pore diameter has a big impact on permeability.If the inner walls of capillaries in porous media are smooth, that is,ξ¯=0, the permeability model can be expressed as(23)K=πDfλmax3+DT128L01+DT3+DT−Df.Equation (23) is identical to the fractal permeability model by Yu and Cheng [14] if only smooth capillaries are taken into account.


## 2.1. Fractal Geometry Theory of Porous Media

The cumulative distributions of pore sizes follow the power law [20], and the number of pores that have diameters greater than or equal to λ can be expressed as [21](1)NL≥λ=λmaxλDf,where N is the number of pores, L is the length scale, and λ is pore diameter; λmin≤λ≤λmax, where λmin and λmax are the minimum and maximum diameters in a porous medium, respectively, and Df is the pore area fractal dimension. Apparently, 1<Df<3 in two dimensions and 2<Df<3 in three dimensions. With (1), the total number of pores between λmin and λmax in a porous medium can be obtained by(2)NtL≥λmin=λmaxλminDf.From the derivative of (1) with respect to λ, we have the number of pores that have diameters ranging between λ and λ+dλ:(3)−dN=DfλmaxDfλ−Df+1dλ.Combining (2) and (3), we can obtain the following expression:(4)−dNNt=DfλminDfλ−Df+1dλ=fλdλ,where fλ=DfλminDfλ−Df+1 is the probability density function associated with pore size distribution. According to the property of the probability density function ∫λminλmaxfλdλ=1, it is clear that λmin/λmaxDf=0 for fractal porous media, and the porous media in nature are proved to satisfy λmin/λmax<10−2 [14], in principle. Furthermore, Yu and Li [21] established a relationship between porosity and pore area fractal dimension based on fractal geometry theory, that is,(5)ϕ=λminλmaxdE−Df,where dE is the Euclidean dimension, dE = 2, 3 in two dimensions and three dimensions, respectively.From (5), Df can be rewritten as(6)Df=dE−lnϕlnλmin−lnλmax.Considering pore tortuosity, the average valueτ¯ in a porous medium is defined as(7)τ¯=LtL0,where Lt and L0 are the actual length and the straight length of a flow path, respectively. The fractal dimension DT of tortuosity indicates the bending extent of streamline paths. DT = 1 means that the flow path is straight; otherwise, it is tortuous; it is clear that the greater the tortuosity, the greater the DT is.The tortuosity fractal dimension can then be described by [22](8)DT=1+lnτ¯lnL0−lnλ¯,where λ¯ is the expectation of pore sizes, that is,(9)λ¯=∫λminλmaxλfλdλ=DfDf−1λmin1−λminλmaxDf−1≈DfλminDf−1.Equations (1)–(9) are from the theoretical base of the present work.


## 2.2. Fractal Characterization of Surface Roughness

Conical roughness elements are studied herein, as shown in Figure1 [23]. Apparently, the fluid flow in rough capillaries is affected by the roughness that is modelled by randomly distributed elements, causing the flow path to bend to a certain extent, like the case of fluid flowing through natural porous media.Figure 1
(a) Capillary bundle model for simulating porous media, whereLt denotes the actual length of capillaries and L0 denotes the linear length of the streamline, which is also the side length of the elementary volume. (b) The roughness elements on the pore wall surface.Due to the tortuous feature of capillary tubes, the actual capillary lengthLtλ is greater than or equal to the straight length L0, and Ltλ fits the power law [14], that is,(10)Ltλ=L0DTλ1−DT.For a rough capillary with a lengthLtλ, according to (3) and (10), the total surface area in an elementary can be obtained as(11)S0=∫λminλmaxπλLtλ−dN=πL0DTDfλmax2−DT2−DT−Df1−λminλmax2−DT−Df.Because of the self-similarity of pore channels, in an elementary volume of a porous medium, the volumeVr of all conical roughness elements is(12)Vr=L031−ϕ.Furthermore, according to (11) and (12), the mean height of the roughness elements can be effectively expressed as(13)h¯=VrS0=L03−DT1−ϕ2−DT−DfπDfλmax2−DT1−λmin/λmax2−DT−Df.The relative mean roughnessξ¯ in a rough cylindrical capillary can then be defined as(14)ξ¯=2h¯λ¯=2L03−DT2−DT−DfDf−1πDf2λmax2−DTλmin1−λmin/λmax2−DT−Df1−ϕ,where λ¯ is given by (9).According to (14), the relative mean roughness of a porous medium can be expressed as a function of porosity, pore area fractal dimension, tortuosity fractal dimension, characteristic length, and pore size (e.g., maximum and minimum pore size). Equation (14) reveals some of the mechanisms related to roughness in rough channels/capillaries.


## 2.3. Permeability Model for Porous Media

Modifying the well-known Hagen–Poiseuille equation, Yu et al. [14] derived an expression for the flow rate through a single curved smooth cylindrical capillary:(15)qλ=π128ΔPLtλλ4μ,where μ is fluid viscosity, ΔP is pressure drop, and λ is pore diameter.Due to the roughness, the effective diameters for fluid flow in capillaries certainly are smaller than those in smooth capillaries. Hence, (15) can be refined as(16)qrλ=π128ΔPLtλλ−2h¯4μ,where h¯ is the effective height of conical roughness elements on average, as described by (13). (16) suggests that due to surface roughness in capillary, the flow rate decreases if the pressure difference maintains constant, or this pressure difference increases if the flow rate is kept constant.Therefore, the total flow in an elementary volume can be expressed as(17)Qr=∫λminλmaxqrλ−dN=π⋅ΔP⋅Df⋅λmaxDf128L0DTμ∫λminλmaxλ−2h¯4λDT−Df−2dλ,where(18)∫λminλmaxλ−2h¯4λDT−Df−2dλ=λmax3+DT−Df3+DT−Df1−λminλmax3+DT−Df−4ξ¯λ¯λmax2+DT−Df2+DT−Df1−λminλmax2+DT−Df+6ξ¯2λ¯2λmax1+DT−Df1+DT−Df1−λminλmax1+DT−Df−4ξ¯3λ¯3λmaxDT−DfDT−Df1−λminλmaxDT−Df+ξ¯4λ¯4λmaxDT−Df−1DT−Df−11−λminλmaxDT−Df−1.In fractal porous media model,1<DT<3,1<Df<3 so −4<2−DT<Df<0. λmin/λmax<0.01, and thus, according to (14) and considering the pore size of the actual porous medium, we get ξ¯2≅0, ξ¯3≅0, ξ¯4≅0.Furthermore, as3+DT−Df>1,2+DT−Df>0, we have λmin/λmax3+DT−Df≅0, λmin/λmax2+DT−Df≅0. Considering the roughness of the capillary wall, the total flow in an elementary volume can be simplified as(19)Qr=πΔPDfλmax2+DT128L0DTμλmax3+DT−Df−4ξ¯λ¯2+DT−Df.According to Darcy’s law, the absolute permeabilityK of a porous medium can be given by [14](20)K=μL0QrΔPA,where A is the cross-sectional area of an elementary volume (cube) with a side length of L0.By substituting (19) into (20) with A=L02, the permeability K∗ of the porous medium with rough pore channels can be obtained as(21)K∗=πDfλmax3+DT128L01+DT3+DT−Df1−Cξ¯,where C is a constant, which is given as(22)C=3+DT−Df2+DT−Df⋅4DfDf−1⋅λminλmax.(21) shows that the permeability of a porous medium with rough capillaries is jointly determined by pore fractal dimension Df, tortuosity fractal dimension DT, structural parameters (L0, λmin, and λmax), and the relative mean roughness ξ¯. The structural parameters are apparently not independent of each other; instead, there is a certain correlation between them. (21) also demonstrates that the larger the relative mean roughness is, the smaller the permeability becomes. This is because the flow resistance increases with the increase in roughness, which is consistent with our understanding about multiphase flow in porous media. The permeability is proportional to the fourth power (3+DT>4) of the maximum pore diameter λmax, so the maximum pore diameter has a big impact on permeability.If the inner walls of capillaries in porous media are smooth, that is,ξ¯=0, the permeability model can be expressed as(23)K=πDfλmax3+DT128L01+DT3+DT−Df.Equation (23) is identical to the fractal permeability model by Yu and Cheng [14] if only smooth capillaries are taken into account.


## 3. Results and Discussion

The tortuosity of the capillary channel aggravates the roughness of the channel surface, so it is very important to determine the tortuosity for the permeability prediction of homogeneous porous media. Therefore, we develop a new tortuosity model by incorporating geometric methods and integral median value theorem [24].


### 3.1. New Tortuosity Model

This work focuses on the porous media that are modelled by square-shaped particles (grains) separating pores. In principle, the tortuosity can be characterized by averaging the measured values for all possible configurations of flow paths in a porous medium. Similar to Yu [25], this work considers two idealized configurations, that is, solid particles being allowed to either contact or overlap. The two possible configurations of the Newtonian incompressible fluid flowing through porous media can be idealized as the cases shown in Figure 2.Figure 2
Two idealized scenarios of streamlines for creeping flow in porous media.The red square in Figure2 is a universal unit cell containing a square, representing particle. A,B, ..., J stand for different positions, and the arrow represents flow paths. Figure 2(a) shows ideal flow paths with particles allowed to overlap randomly (the overlapping is not illustrated here just for a clear visualization), where a is the side length of the square particle and d is the width of pores between grains (i.e., pore size). Since overlapping may occur at any position within a unit cell, the flow path between B and D might have different pathways. For instance, B⟶C⟶D is one of the representative paths. The position C locates on the line between D and E and can be altered arbitrarily (i.e., in Figure 2(c), C can appear anywhere on the line between D and E), which suggests that there might be many flow pathways between B and D through C. Figure 2(b) illustrates another flow path along which there is no overlapping between particles. G⟶H⟶I⟶J in Figure 2(d) is the representative flow path among the unit cells.With the configuration in Figure2(c), the porosity can be estimated by(24)ϕ=SpSt=1−aa+d2,where St is the total area of the unit cell and SP is the pore area in a unit cell.Equivalently, we can have the following formula to link porosity to pore size:(25)d=a11−ϕ−1.It is clear from (25) that as ϕ⟶1, we have d⟶∞,or oppositely when ϕ⟶0, the pore size tends to zero. Namely, when ϕ⟶1, there is no space for particles to occupy, which indicates that the pore size turns to be infinite.Considering the incompressible Newtonian fluid and the thin fluid layer on particle surfaces, the pore tortuosity of a porous medium can then be measured by averaging the tortuosities of all flow paths, that is,(26)τ¯=1n∑iτi,where n is the total number of fluid flow paths in consideration, and τi is the tortuosity of the ith flow path.The particles in natural porous media can be regarded as randomly distributed and overlapped—some particles overlap; meanwhile, others do not. Yu and Li [25] considered “unrestricted overlapping of particles” and studied streamlines’ tortuosity. They drew a conclusion that the period of streamlines may not be the same as that shown in Figure 1(a) and it may be in the scale of pore size BE; B⟶C⟶D is any one of the fluid paths, and the radial length of the flow is BE of length d, as shown in Figure 2(c). These form the fundamental point of view for modelling flow configuration in this work.In Figure2(c), since the position of point C between D and E is arbitrary, we let the height CE of a streamline path be x. Therefore, the length g(x) of streamline B⟶C⟶D can be expressed as(27)gx=BC+CD=a2φ−12+x2+a2−x,x∈0,a2,where φ≡1/1−ϕ.According to the integral median value theorem [24], considering x as a random variable, the expectation of g(x) within [0, a/2] can be computed by(28)gζ=1a/2∫0a/2gxdx=a2γ+aγ2−1/4ln1/2+γγ2−1/4+a4,where γ≡1/4+φ−12, gζ is determined by porosity ϕ only, which represents the average length of all flow paths between B and D passing through variable positions at C along the unit cell wall.The mean tortuosityτ1 in a porous medium composed of unrestricted overlapping particles can then be given by(29)τ1=gζd=2γ+14γ2−1/4+γ2−1/4⋅ln1/2+γγ2−1/4.For the streamline that has no surrounding particles overlapping, the flow paths associated can be assumed to be the idealized paths, as shown in Figure2(c).For an idealized path in Figure2(d), the streamline HI of length a/2 is perpendicular to the horizontal direction, and the whole length of the streamline G⟶H⟶I⟶J should be taken into account in tortuosity, denoted by τ2, that is,(30)τ2=HI+a+da+d=1+12φ.The representative tortuosity of a porous media can then be obtained by averaging the estimates for the two configurations. So, the mean tortuosityτ¯ of a porous media can be rewritten as(31)τ¯=τ1+τ22=12+14φ+2γ+18γ2−1/4+γ2−1/42⋅ln1/2+γγ2−1/4,where γ≡1/4+φ−12, ϕ≡1/1−ϕ, and ϕ is the porosity.Note that (31) is an approximated tortuosity over all possible flow paths in a two-dimensional porous medium of square-shaped particles. The equation describes the tortuosity as a function of porosity only, and hence, there is no extra empirical constants needed to be determined in advance. Since the pore microchannel structure is assumed to be isotropic and self-similar, the streamline characteristics within each cross section of the porous medium are similar; therefore, (31) can also be used to simulate the streamline paths of isotropic three-dimensional porous media. The proposed model of (31) may look complex, but it is powerful to boost our understanding of physical mechanisms of tortuous flow paths in porous media.


### 3.2. Permeability Prediction for Real Porous Media

To verify the rationality and accuracy of (21), nine digital core samples, including carbonate (C), sandstone (S), synthetic silica (SS), and sand-pack (SP) models, were selected for this study. The three-dimensional CT images and permeability data come from Imperial College London (Table 1 [26]). The selected digital cores have porosities ranging from 0.168 (C2) to 0.4297 (SS1) and pore diameters ranging from 1019 μm (S9) to 3000 μm (SP1).Table 1
Basic parameters of the digital core samples (pores in black) from Imperial College London.
Core sampleBinary imagesSide length (μm)Porosity (ϕ)Permeability (μm2)C111400.2331.0876C221380.1680.0714S21486.80.2463.8470S51199.10.2114.5902S615300.2410.8305S71440.90.2516.8749S91019.40.2222.1949SS111550.42977.1256SP13000.60.33149.7409In Table2, λmin and λmax denote the minimum and maximum pore diameters, respectively, for each sample listed in Table 1, obtained using the box-counting method [27]. τ¯ denotes the average tortuosity of fluid paths, calculated by (31), and Df denotes the pore fractal dimension from (6), while DT denotes the tortuosity fractal dimension by (8), and ξ¯ denotes the relative mean roughness by (14).Table 2
Fractal characterization of the digital core samples listed in Table1.
Core sampleλmin (μm)λmax (μm)τ¯DfDTξ¯C12.8596.902.65632.68491.19330.0187C25.3566.783.46052.60481.24950.0257S24.96133.812.54762.67791.19690.0222S54.00123.912.87172.63831.22390.0260S65.10173.402.59622.67161.20110.0228S74.80148.892.50892.68301.19350.0217S93.4098.542.75852.65151.21480.0247SS13.8577.001.73132.81391.11240.0105SP110.00330.072.05732.75091.14970.0156As shown in Tables1 and 2, the greater the porosity is, the larger the pore fractal dimension and the smaller the tortuosity and tortuosity fractal dimension are. This observation is in agreement with the conclusion drawn in previous sections. Table 3 compares the predicted results from different models, particularly with Yu and Cheng [14] ((21)with ξ¯=0) and Xu et al. [27]:(32)KX=132τ¯2−Df4−Dfϕ1−ϕλmax2.Table 3
Comparison of the predicted permeability values (μm2) of the digital core samples.
Core sampleKeK∗KKXC11.08761.82101.84002.0857C20.07140.07810.08070.1206S23.84703.84603.90604.5171S54.59023.96904.02804.6006S610.830510.552211.449511.8306S76.87496.86416.94557.4722S92.19492.25232.28712.6367SS17.12561.01921.03381.1661SP149.740945.256545.731949.6788In Table3, Ke denotes the permeability data from Imperial College London, K∗ the permeability predicted using our rough capillary permeability model (i.e., Equation (21)), K the value predicted by the permeability model with smooth capillary walls (Equation (23)), and KX the value predicted by Equation (32).To further evaluate the reliability of our permeability prediction model, it is necessary to analyze prediction accuracy statistically. Relative error and mean absolute percentage error (MAPE) were used to quantify the error between the models and the experiment, that is,(33)ε=xexp−xprexexp,MAPE=100n∑i=1nxiexp−xiprexiexp,where xexp denotes the experimental data, xpre denotes the predicted value, ε quantifies the relative error between experimental data and predicted values, n represents the total number of experimental data points, and xiexp and xipre represent the ith experimental value and the corresponding predicted value, respectively, i=1,2,…,n.As shown in Figure3, the relative errors between the permeability values predicted by different models and the experimental data were large (>50%) for the carbonate (e.g., C1) and synthetic silica (SS1). This results from high heterogeneity of pore structures when using the box-counting method to determine the maximum pore diameter, which imposes the most significant influence on permeability. Equations (21) and (23) predicted the permeability values reasonably well, with only three samples as exceptions that have prediction errors greater than 10%. Furthermore, Equation (21), which incorporate relative mean roughness as a major impact element, predicted a lower permeability than Equation (23). The relative errors for the five samples predicted by Equation (32) all exceed 10%. Table 4 shows that the MAPE value between the predicted permeability value using Equation (21) and the experimental value is the smallest at 0.2116. Therefore, the proposed model for rough porous media (Equation (21)) leads to more accurate permeabilities. This is due to the improvement of the average tortuosity model (Equation (31)), which improves the accuracy of the tortuosity fractal dimension in the permeability model, and relative mean roughness is additionally introduced as a structural parameter into the prediction model, which could more realistically reflect the influence on the permeability.Figure 3
Relative errors of the permeability values, whereε∗,ε,εX denote the relative errors between the predicted values K∗,K,KX and the experimental permeability Ke, respectively.Table 4
Statistical analysis of the prediction precision of the permeability listed in Table3.
MetricEquation (21)Equation (23)Equation (32)MAPE0.21160.22280.3335


## 3.1. New Tortuosity Model

This work focuses on the porous media that are modelled by square-shaped particles (grains) separating pores. In principle, the tortuosity can be characterized by averaging the measured values for all possible configurations of flow paths in a porous medium. Similar to Yu [25], this work considers two idealized configurations, that is, solid particles being allowed to either contact or overlap. The two possible configurations of the Newtonian incompressible fluid flowing through porous media can be idealized as the cases shown in Figure 2.Figure 2
Two idealized scenarios of streamlines for creeping flow in porous media.The red square in Figure2 is a universal unit cell containing a square, representing particle. A,B, ..., J stand for different positions, and the arrow represents flow paths. Figure 2(a) shows ideal flow paths with particles allowed to overlap randomly (the overlapping is not illustrated here just for a clear visualization), where a is the side length of the square particle and d is the width of pores between grains (i.e., pore size). Since overlapping may occur at any position within a unit cell, the flow path between B and D might have different pathways. For instance, B⟶C⟶D is one of the representative paths. The position C locates on the line between D and E and can be altered arbitrarily (i.e., in Figure 2(c), C can appear anywhere on the line between D and E), which suggests that there might be many flow pathways between B and D through C. Figure 2(b) illustrates another flow path along which there is no overlapping between particles. G⟶H⟶I⟶J in Figure 2(d) is the representative flow path among the unit cells.With the configuration in Figure2(c), the porosity can be estimated by(24)ϕ=SpSt=1−aa+d2,where St is the total area of the unit cell and SP is the pore area in a unit cell.Equivalently, we can have the following formula to link porosity to pore size:(25)d=a11−ϕ−1.It is clear from (25) that as ϕ⟶1, we have d⟶∞,or oppositely when ϕ⟶0, the pore size tends to zero. Namely, when ϕ⟶1, there is no space for particles to occupy, which indicates that the pore size turns to be infinite.Considering the incompressible Newtonian fluid and the thin fluid layer on particle surfaces, the pore tortuosity of a porous medium can then be measured by averaging the tortuosities of all flow paths, that is,(26)τ¯=1n∑iτi,where n is the total number of fluid flow paths in consideration, and τi is the tortuosity of the ith flow path.The particles in natural porous media can be regarded as randomly distributed and overlapped—some particles overlap; meanwhile, others do not. Yu and Li [25] considered “unrestricted overlapping of particles” and studied streamlines’ tortuosity. They drew a conclusion that the period of streamlines may not be the same as that shown in Figure 1(a) and it may be in the scale of pore size BE; B⟶C⟶D is any one of the fluid paths, and the radial length of the flow is BE of length d, as shown in Figure 2(c). These form the fundamental point of view for modelling flow configuration in this work.In Figure2(c), since the position of point C between D and E is arbitrary, we let the height CE of a streamline path be x. Therefore, the length g(x) of streamline B⟶C⟶D can be expressed as(27)gx=BC+CD=a2φ−12+x2+a2−x,x∈0,a2,where φ≡1/1−ϕ.According to the integral median value theorem [24], considering x as a random variable, the expectation of g(x) within [0, a/2] can be computed by(28)gζ=1a/2∫0a/2gxdx=a2γ+aγ2−1/4ln1/2+γγ2−1/4+a4,where γ≡1/4+φ−12, gζ is determined by porosity ϕ only, which represents the average length of all flow paths between B and D passing through variable positions at C along the unit cell wall.The mean tortuosityτ1 in a porous medium composed of unrestricted overlapping particles can then be given by(29)τ1=gζd=2γ+14γ2−1/4+γ2−1/4⋅ln1/2+γγ2−1/4.For the streamline that has no surrounding particles overlapping, the flow paths associated can be assumed to be the idealized paths, as shown in Figure2(c).For an idealized path in Figure2(d), the streamline HI of length a/2 is perpendicular to the horizontal direction, and the whole length of the streamline G⟶H⟶I⟶J should be taken into account in tortuosity, denoted by τ2, that is,(30)τ2=HI+a+da+d=1+12φ.The representative tortuosity of a porous media can then be obtained by averaging the estimates for the two configurations. So, the mean tortuosityτ¯ of a porous media can be rewritten as(31)τ¯=τ1+τ22=12+14φ+2γ+18γ2−1/4+γ2−1/42⋅ln1/2+γγ2−1/4,where γ≡1/4+φ−12, ϕ≡1/1−ϕ, and ϕ is the porosity.Note that (31) is an approximated tortuosity over all possible flow paths in a two-dimensional porous medium of square-shaped particles. The equation describes the tortuosity as a function of porosity only, and hence, there is no extra empirical constants needed to be determined in advance. Since the pore microchannel structure is assumed to be isotropic and self-similar, the streamline characteristics within each cross section of the porous medium are similar; therefore, (31) can also be used to simulate the streamline paths of isotropic three-dimensional porous media. The proposed model of (31) may look complex, but it is powerful to boost our understanding of physical mechanisms of tortuous flow paths in porous media.


## 3.2. Permeability Prediction for Real Porous Media

To verify the rationality and accuracy of (21), nine digital core samples, including carbonate (C), sandstone (S), synthetic silica (SS), and sand-pack (SP) models, were selected for this study. The three-dimensional CT images and permeability data come from Imperial College London (Table 1 [26]). The selected digital cores have porosities ranging from 0.168 (C2) to 0.4297 (SS1) and pore diameters ranging from 1019 μm (S9) to 3000 μm (SP1).Table 1
Basic parameters of the digital core samples (pores in black) from Imperial College London.
Core sampleBinary imagesSide length (μm)Porosity (ϕ)Permeability (μm2)C111400.2331.0876C221380.1680.0714S21486.80.2463.8470S51199.10.2114.5902S615300.2410.8305S71440.90.2516.8749S91019.40.2222.1949SS111550.42977.1256SP13000.60.33149.7409In Table2, λmin and λmax denote the minimum and maximum pore diameters, respectively, for each sample listed in Table 1, obtained using the box-counting method [27]. τ¯ denotes the average tortuosity of fluid paths, calculated by (31), and Df denotes the pore fractal dimension from (6), while DT denotes the tortuosity fractal dimension by (8), and ξ¯ denotes the relative mean roughness by (14).Table 2
Fractal characterization of the digital core samples listed in Table1.
Core sampleλmin (μm)λmax (μm)τ¯DfDTξ¯C12.8596.902.65632.68491.19330.0187C25.3566.783.46052.60481.24950.0257S24.96133.812.54762.67791.19690.0222S54.00123.912.87172.63831.22390.0260S65.10173.402.59622.67161.20110.0228S74.80148.892.50892.68301.19350.0217S93.4098.542.75852.65151.21480.0247SS13.8577.001.73132.81391.11240.0105SP110.00330.072.05732.75091.14970.0156As shown in Tables1 and 2, the greater the porosity is, the larger the pore fractal dimension and the smaller the tortuosity and tortuosity fractal dimension are. This observation is in agreement with the conclusion drawn in previous sections. Table 3 compares the predicted results from different models, particularly with Yu and Cheng [14] ((21)with ξ¯=0) and Xu et al. [27]:(32)KX=132τ¯2−Df4−Dfϕ1−ϕλmax2.Table 3
Comparison of the predicted permeability values (μm2) of the digital core samples.
Core sampleKeK∗KKXC11.08761.82101.84002.0857C20.07140.07810.08070.1206S23.84703.84603.90604.5171S54.59023.96904.02804.6006S610.830510.552211.449511.8306S76.87496.86416.94557.4722S92.19492.25232.28712.6367SS17.12561.01921.03381.1661SP149.740945.256545.731949.6788In Table3, Ke denotes the permeability data from Imperial College London, K∗ the permeability predicted using our rough capillary permeability model (i.e., Equation (21)), K the value predicted by the permeability model with smooth capillary walls (Equation (23)), and KX the value predicted by Equation (32).To further evaluate the reliability of our permeability prediction model, it is necessary to analyze prediction accuracy statistically. Relative error and mean absolute percentage error (MAPE) were used to quantify the error between the models and the experiment, that is,(33)ε=xexp−xprexexp,MAPE=100n∑i=1nxiexp−xiprexiexp,where xexp denotes the experimental data, xpre denotes the predicted value, ε quantifies the relative error between experimental data and predicted values, n represents the total number of experimental data points, and xiexp and xipre represent the ith experimental value and the corresponding predicted value, respectively, i=1,2,…,n.As shown in Figure3, the relative errors between the permeability values predicted by different models and the experimental data were large (>50%) for the carbonate (e.g., C1) and synthetic silica (SS1). This results from high heterogeneity of pore structures when using the box-counting method to determine the maximum pore diameter, which imposes the most significant influence on permeability. Equations (21) and (23) predicted the permeability values reasonably well, with only three samples as exceptions that have prediction errors greater than 10%. Furthermore, Equation (21), which incorporate relative mean roughness as a major impact element, predicted a lower permeability than Equation (23). The relative errors for the five samples predicted by Equation (32) all exceed 10%. Table 4 shows that the MAPE value between the predicted permeability value using Equation (21) and the experimental value is the smallest at 0.2116. Therefore, the proposed model for rough porous media (Equation (21)) leads to more accurate permeabilities. This is due to the improvement of the average tortuosity model (Equation (31)), which improves the accuracy of the tortuosity fractal dimension in the permeability model, and relative mean roughness is additionally introduced as a structural parameter into the prediction model, which could more realistically reflect the influence on the permeability.Figure 3
Relative errors of the permeability values, whereε∗,ε,εX denote the relative errors between the predicted values K∗,K,KX and the experimental permeability Ke, respectively.Table 4
Statistical analysis of the prediction precision of the permeability listed in Table3.
MetricEquation (21)Equation (23)Equation (32)MAPE0.21160.22280.3335


## 4. Conclusion

In this work, we propose a pore tortuosity model for porous media and a relative mean roughness model for rough capillary walls. Then, taking pore tortuosity and relative mean roughness as important parameters, the permeability model of porous media is derived according to fractal theory, which improves the predicted value accuracy of permeability. The major conclusions can be summarized as follows:(1)
Based on the fractal theory, a so-called relative mean roughness of pore walls in porous media is proposed in Equation (14) to characterize fractal features of the pore channels, which is mathematically correlated to many parameters, including porosity, pore area fractal dimension, tortuosity fractal dimension, digital core sample scale, and the maximum and minimum pore diameters. The relative mean roughness model shows that there is also a power-law relationship between relative mean roughness and porosity.(2)
Assuming that the porous medium is composed of square particles, under the special structural configurations (Figure2), the average tortuosity (Equation (31)) is newly expressed as a function of porosity without other empirical constants involved unlike many other formulas in the literature. The derivation of this formula relies on the integral median value theorem; consequently, the calculation accuracy of the tortuosity model is improved considerably. Moreover, using the tortuosity as a parameter leads to improve the accuracy of the tortuosity fractal dimension (Equation (8)).(3)
We established a novel fractal model of permeability (Equation (21)) based on fractal theory. With regard to the model, we can conclude that permeability is in direct relation to pore area fractal dimension, tortuosity fractal dimension, digital core sample scale, pore diameter, and relative mean roughness. The model also suggests that the larger the relative mean roughness is, the smaller the permeability becomes. Compared with other permeability models, the calculation of the tortuosity fractal dimension in this model utilizes the newly derived average tortuosity model (Equation (31)), and the parameter of relative mean roughness is newly added to the model. From the comparison of the prediction results of different permeability models, it can be seen that new permeability model has high accuracy.In future research, we will study the average tortuosity model of 3D porous media, considering the streamline situation that is closer to the actual situation. The prediction results of permeability will be improved by improving the calculation accuracy of tortuosity.



---

*Source: 8088151-2022-08-23.xml*