# Permeability Estimation Based on the Geometry of Pore Space via Random Walk on Grids

Authors: Tongchao Nan; Jichun Wu; Kaixuan Li; Jianguo Jiang
Journal: Geofluids
 (2019)

---

## Abstract

In the literature, the mean penetration depth (MPD) calculated by “walk on spheres” or “walk on cubes” was used to quickly estimate the intrinsic permeability of digitized porous media. However, these two methods encounter difficulties such as irregular boundaries and the determination of arrivals at a boundary. In this study, an MPD method that is based on a more flexible “random walk on grid” (WOG) is explored. Moreover, the accurate MPDs for the pores of simple shapes are derived with Green’s functions to validate the WOG-based MPD. The results suggested that MPDs based on Green’s functions and WOG are consistent with each other; the factorC in the permeability expression is slightly dependent on roundness of the cross sections and is approximately 1.125 on average, according to analytical and numerical results. In a synthetic complex pore, the permeability estimated by WOG is comparable to, but greater than, the estimate based on the pore-scale dynamics simulation in COMSOL.


---

## Body

## 1. Introduction

Permeability, which is a measure of the ability of a porous material to allow fluid to pass through it, is one of the most important properties of geologic formations to hydrogeologists, petrophysicists, and other geologic fluid researchers. Since permeability (or hydraulic conductivity) greatly influences the fluid movements and flow fields of groundwater, it is generally the fundamental factor in all groundwater-related problems, such as groundwater exploitation, surface-groundwater interaction, contaminant transport, agricultural irrigation, and site selection for radioactive waste disposal [1, 2].As a macroscopic parameter of formations, permeability is closely related to pore microstructures. In addition to the direct measurement of samples, indirect measurement is also important for different environmental and technological applications. Some studies have been devoted to establishing the relationship between microscopic geometry and macroscopic permeability by theoretical models and numerical calculations [3–8], which are becoming increasingly valuable for the growing availability of microscopic tomographic data, such as nuclear magnetic resonance, X-ray microtomography, and scanning-electron-microscopic imaging ([9–12]; among many others).For pores of some of the simplest geometries, the permeability can be calculated by the exact solution to the Poiseuille flow equation. For example, in a straight tube of circle-shaped cross sections, one obtainsk=R2/8 by combining three equations: Q=πR4/8μΔP/L (Hagen-Poiseuille law), Q=KΔP/ρgLπR2 (Darcy’s law), and K=kρg/μ (relation between permeability and conductivity) (see, e.g., [13]). For tubes of elliptical cross sections, k=a2b2/4a2+b2, where a and b are the semiaxes [14]; similar results are available for several other shapes [15, 16]. For media of higher complexity, for example, k=cKn3/1−n2d2 in a packed bed of solids, where d is the mean grain diameter, n is the porosity, and cK is a parameter that is dependent on cross-sectional geometry [6, 17, 18].However, more powerful tools are required in media of general microstructure. By analyzing diffusion processes in porous media, Torquato [8] linked the permeability tensor of the media to a so-called trapping constant, which is not easy to evaluate and thus is not practical. However, the idea of using diffusion or Brownian processes to reflect the properties of porous media is useful. Based on the Brownian motion model and the Debijf-Brinkman equation (see, e.g., [19]), Hwang et al. [5] found that the mean penetration depth (MPD, denoted by D), which is the expected depth that a Brownian motion can enter before reaching solid surface for the first time, is a good measure of the permeability of the medium and suggested that k=D2. Simonov and Mascagni [7] took porosity into account and suggested that k=CSMnD2, where n is the porosity and cSM is a factor. They estimated CSM=k/nD2≈6.63 by using the Poiseuille law-based permeability k=0.0351d2 in an ideal, straight pore of the square cross section with a side length that is equal to d, and n=0.5. It is still not clear whether k=CSMnD2 is also valid for pores of other geometries. Sabelfeld [20] used a spectral projection method to avoid the numerical calculation on overlapping discs and spheres which, however, applies to simplified pore surfaces and is unsuitable for real media with rough surfaces.The equations above are all in a quadratic form; i.e.,k=CL2, where L is some type of characteristic length (e.g., hydraulic radius, mean pore size, mean grain diameter, or MPD), and the coefficient C may depend on the cross section and/or porosity ([13]; and references therein).If one uses MPD as the characteristic length, the coefficientC may depend on both cross-section shape and porosity. The influence of porosity on C was also investigated by Simonov and Mascagni [7] using “walk on cubes” and “walk on spheres.” However, “walk on cubes” is only suitable for voids of polyhedron and is inefficient for curly boundaries [21]. “Walk on spheres” has difficulty in determining exactly when and where the random walk ends at a boundary. The difficulty leads to significant error and generates some bias (Milstein and Tretyakov 1999; Deaconu and Lejay 2006). To avoid these problems, in this study, we extend the MPD method by using a more flexible “random walk on grid” (WOG), which is able to easily determine the arrival at boundaries and handle irregular surfaces, including cambering, narrow gaps, and wedge outs.To validate the effectiveness of the WOG-based MPD method and to investigate the extent ofC’s dependence on the roundness of pore cross sections, we derived the theoretical values of MPD and C using Green’s function in the pores of several basic geometries and compared them to their WOG-based counterparts. Simonov and Mascagni [7] realized the possible weak dependence of C on the cross-sectional geometry of void space, but they did not investigate the problem further. Here, the dependence is investigated through C derived by Green’s functions, as well as the C estimated by WOG.To date, the permeability calculated by the MPD method has not been verified in a porous medium against dynamic simulation or experimental results in the literature. A synthetic porous sample is constructed for permeability estimation using COMSOL simulation and the WOG-based MPD method.The paper is organized as follows: In Section2, the MPD method based on WOG and Green’s functions, as well as COMSOL simulation, is introduced. In Section 3, the results of the MPD estimation in simple pores of different geometries are reported, and the factor C is calibrated for each geometry. Furthermore, the average C factor is used to estimate the permeability of a synthetic pore, and the result is verified by a COMSOL model. In Section 4, the results are discussed and conclusions are drawn.


## 2. Methodology

### 2.1. Mean Penetration Depth (MPD) and MPD Calculation via Green’s Function

Given a porous medium, let a Brownian motion initiate atx′ (denoted by x0=x′) on its outer boundary plane (say, z=0; Figure 1) and enter the medium until the walker hits a solid surface (i.e., the inner boundary, represented by blue surfaces in Figure 1) at a random position, x∗. The depth of the hit point relative to the outer boundary plane (i.e., zω; ω is the ensemble index) is called the penetration depth of the Brownian motion, denoted by lx′,ω. The mean penetration depth (MPD) is the ensemble mean of the penetration depths in Monte Carlo experiments, starting from an area on the outer boundary (termed “departure region” Ω, i.e., the green area in Figure 1), i.e., DB=1/S∫x′∈Ωlx′,ωdx′, where  stands for expectation with respect to the ensemble index ω, DB denotes the MPD estimated by the Brownian motion (or random walk), and S is the area of Ω.Figure 1
A diagram of random walk in void space. The outer boundary is the face of the medium; the inner boundary consists of the void-solid interface. Black line shows the path of a random walk realization.In fact, the Brownian motion in void space, as mentioned above, is a diffusion process that is reflected back to void space at the outer boundary and is killed at the inner boundary. Given a starting pointx0=x′, the spatial density cx∣x′ (i.e., mean times of occurrence in space) of the diffusion satisfies the following equation:
(1a)Δcx∣x′+δx−x′=0,(1b)cx∣x′=0,x∈Γ,(1c)∂cx∣x′∂ne=0,x∈Ω,where x′=x′,y′,z′∈Ω and x=x,y,z, δ is the Dirac function, and ne is the vector normal to outer boundary Ω. Equation (1b) is the absorbing boundary at the void-solid interface Γ. Equation (1c) is the reflection (no-flux) boundary at the departure region Ω (outer boundary).One may notice that (1) is actually a Green’s function problem. Its solution cx∣x′ is the Green’s function of the following Poisson equation:
(2a)Δux=−fx,(2b)DirichletboundaryatΓ,(2c)NeumannboundaryatΩ.For some simple geometries, the analytical expressions ofcx∣x′ can be found or easily derived from the libraries of Green’s functions in the literature (e.g., Cole et al. [22]). When cx∣x′ is known, the probability distribution of the hit point x∗ at the inner boundary can be expressed as (see, e.g., Deaconu and Lejay [23])
(3)px∗∣x0=x′=−∂cx∗∣x′∂n,x∗∈Γ,where n is the vector normal to Γ and pointing to solid. Hence, the penetration depth and MPD can be calculated via cx∣x′. That is, the penetration depth of point x′ is lex′=lx′,ω=−∫Γz∂cx∣x′/∂n∣x∈Γdx, which will be further averaged on the departure region Ω to calculate the MPD, i.e.,
(4)DG=1S∫x′∈Ωlex′dx′=−1S∫x′∈Ωdx′∫Γz∂cx∣x′∂n∣x∈Γdx,where S is the area of Ω.We consider the pores of three simple geometries for which Green’s functioncx∣x′ is available, i.e., parallel planes, straight square, and circular tubes (Figures 2(a)–2(c)). Here, we used the “isoperimetric ratio” φ as a measure of the pore cross-sectional shape. φ=lc2/A is used to measure how greatly the shape differs from a circle, where lc is the length of a curve (perimeter) and A is the corresponding area. φ has a minimum value 4π for the circle and higher values for all other shapes [24]. For parallel planes φ=∞, since the length of planes is infinite and the plane distance d is finite; for the square cross section, φ=4d2/d2=16; for the equilateral triangle, φ=3d2/3d2/4=123≈20.8; for the ellipse with semiaxes a and ba≥b, 4π≤φ<∞ and φ≈15.0 as a=2b.Figure 2
Void spaces (a) between parallel planes and in straight tubes of (b) a square cross section; (c) a circular cross section; (d) an elliptic cross section; and (e) an equilateral triangular cross section. (f) An ideal porous medium sample consisting of solid and many straight pores. The green area indicates the outer boundary.
(a)
(b)
(c)
(d)
(e)
(f)For parallel planes of distanced (Figure 2(a)), the corresponding Green’s function is expressed as [22]
(5)cx,z∣x′,z′=0=2π∑n=1∞1ne−nπz/dsinnπxdsinnπx′d.Applying Equation (4), we obtain the MPD based on Green’s function, i.e.,
(6)DG=7ζ3dπ3≈0.2714d,where ζ∙ is the zeta function and the subscript “G” stands for Green’s function. For straight tubes of square cross sections (with side-length d, Figure 2(b)), Green’s function is
(7)cx,y,z∣x′,y′,z′=0=4πd∑m=1∞∑n=1∞exp−πn2+m2z/dn2+m2sinmπxdsinmπx′dsinnπydsinnπy′d,and the MPD based on Green’s function is
(8)DG=16dπ5∑m=1∞∑n=1∞1n2n2+m23/21−−1m1−−1n2≈16dπ5⋅3.2535=0.1701d.For straight tubes of circular cross sections (with radiusR, Figure 2(c)), Green’s function is written as
(9)cr,z∣r′,z′=0=1πR∑m=1∞1βmJ12βmexp−βmzRJ0βmrRJ0βmr′R.The MPD based on Green’s function is(10)DG=4R∑m=1∞1βm3≈0.3235R,where J0 and J1 are Bessel functions (of zeroth and first orders) of the first kind; βm is the m‐th root of J0βm=0.No analytical Green’s functions are found for the remaining examples in Figure2; therefore, the numerical simulation of the Brownian motion will be used to directly estimate the MPD.Taking into account the impact of porosity, let us consider a porous medium consisting of many straight tubes (Figure2(f)), whose porosity n is not equal to 1. For simplicity, assume that all of the pores are the same. Recall that k represents the permeability of a single pore, DB is the MPD with respect to a single pore, and k=CDB2. One may note that the permeability of such a medium is
(11)kb=nk=nCDB2,where “b” stands for bunch. Let Db be the effective MPD after being averaged on the whole area S; it is obvious that Db=nDB+1−n⋅0=nDB. Let kb=CbDb2. Then, nk=Cbn2DB2, i.e., Cb=k/n/DB2. Recall (C is the factor for porosity n=1), one has Cb=C/n. Thus,
(12)kb=CDb2n,which is valid for media comprising uniform straight tubes. We will check this formula in the next section to test its applicability. The impact of tortuosity on permeability is not discussed but can be considered by some empirical equations in the literature (e.g., [13]).


### 2.2. Random Walk on Grid (WOG)

Green’s functions are accessible to only limited, simple geometries and are not practical for general pores. Thus, one has to rely on the numerical simulation of the Brownian motion or random walk. To simulate a Brownian motion in the medium, here we employ a random walk on grid (WOG) due to its ease in determining hit occurrences and its flexibility of irregular inner boundaries. In homogeneous void space, WOG is equivalent to a simple lattice random walk. Here, implementation of WOG is only briefly introduced in Figure3. In step (a), the domain can be discretized with a regular/irregular block-/point-centered grid. In this study, we use regular block-centered grids; i.e., each node represents a block or cell of the same size. In step (b), given that the distances from the current node to the neighboring nodes in six directions (i.e., x±, y±, and z±) are ∆x±, ∆y±, and ∆z± in the Cartesian coordinate and that px±, py±, and pz± denote the probability that the walker moves from the current node to six neighboring nodes, the probability of the neighboring node in x+ direction can be estimated as
(13)px+=1/Δx+1/Δx++1/Δx−+1/Δy++1/Δy−+1/Δz++1/Δz−,which is similar for other directions. If the grid is nonuniform, the probabilities are location-dependent and direction-dependent. The probabilities to neighbors at all nodes are calculated and stored before random walks are simulated. In step (c), every node x′ in the departure area will be selected to serve as the starting node in the Monte Carlo loop for NMC times (Figure 3). After completion of the Monte Carlo loop, the node will be excluded from the departure area and will not be used again in the following simulation. If the original departure area contains NS nodes after discretization, a random walk will be simulated for NS×NMC times, and the number of hit nodes is NS×NMC in total. Note that each realization of the random walk has an uncertain number of steps (in the random walk loop), and the location of its hit node is also uncertain before the realization is complete. When all of the nodes in the original departure area have served as starting nodes, the simulation procedure ends and the depths of all hit nodes will be averaged to estimate the MPD. More details of WOG algorithms and the superiority of WOG over “walk on cubes” and “walk on spheres” can be found in Nan and Wu [21].Figure 3
A flowchart of the MPD estimation via WOG.We evaluated the permeability in the pores of some simple geometries, including parallel planes and void of square, circular, elliptic, and equilateral triangular cross sections (Figures3(a)–3(e)). For the void space of these shapes, the theoretical expressions of permeability also exist (Table 1). In all of these pores, the porosities are all equal to 1. Factor C can be estimated using CB=k/DB2 in each case after DB is estimated by WOG.Table 1
Permeability, MPD, andC factor estimated by WOG and Green’s function for voids of simple geometries.
Void shape
Reference length
Theoreticalk∗
MPDDB¢
MPDDGǂ
CB¢
CGǂ
φ£
Parallel plane (Figure2(a))
Plane distanced
d212
0.2763d
0.2714d
1.091
1.131
∞
Square (Figure2(b))
Side-lengthd
0.0351d2
0.1751d
0.1701d
1.145
1.123
16
Circle (Figure2(c))
RadiusR
R28
0.3318R
0.3235R
1.135
1.194
4π
Ellipse (Figure2(d))
Semiaxes:a, b
a2b24a2+b2
0.4210b (as a=2b)
—
1.128
—
4π~∞ (15 as a=2b)
Equilateral triangle (Figure2(e))
Side-lengthd
d280
0.1054d
—
1.125
—
123
∗Theoretical k is obtained according to pressure drops and theoretical flow rates in Drazin and Norman [25]. ¢DB is the MPD estimated by random walk on grid (WOG); CB=k/DB2. ǂDG is the MPD evaluated by Green’s function cx∣x′ by Equations (6), (8), and (10); CG=k/DG2. £φ according to Berger [24].


### 2.3. COMSOL Simulation

COMSOL is a general-purpose platform software for modeling engineering problems. It is able to simulate multiple physical processes and models, including the flow of fluid. In this study, we use COMSOL to solve the steady-state Stokes flow in the pore space. The void-solid interface is treated as a no-slip boundary. The prescribed pressure boundaries are imposed onto two opposite faces of the parallelepiped and no-flux boundaries onto the other four lateral faces (due to the symmetry of the medium).


## 2.1. Mean Penetration Depth (MPD) and MPD Calculation via Green’s Function

Given a porous medium, let a Brownian motion initiate atx′ (denoted by x0=x′) on its outer boundary plane (say, z=0; Figure 1) and enter the medium until the walker hits a solid surface (i.e., the inner boundary, represented by blue surfaces in Figure 1) at a random position, x∗. The depth of the hit point relative to the outer boundary plane (i.e., zω; ω is the ensemble index) is called the penetration depth of the Brownian motion, denoted by lx′,ω. The mean penetration depth (MPD) is the ensemble mean of the penetration depths in Monte Carlo experiments, starting from an area on the outer boundary (termed “departure region” Ω, i.e., the green area in Figure 1), i.e., DB=1/S∫x′∈Ωlx′,ωdx′, where  stands for expectation with respect to the ensemble index ω, DB denotes the MPD estimated by the Brownian motion (or random walk), and S is the area of Ω.Figure 1
A diagram of random walk in void space. The outer boundary is the face of the medium; the inner boundary consists of the void-solid interface. Black line shows the path of a random walk realization.In fact, the Brownian motion in void space, as mentioned above, is a diffusion process that is reflected back to void space at the outer boundary and is killed at the inner boundary. Given a starting pointx0=x′, the spatial density cx∣x′ (i.e., mean times of occurrence in space) of the diffusion satisfies the following equation:
(1a)Δcx∣x′+δx−x′=0,(1b)cx∣x′=0,x∈Γ,(1c)∂cx∣x′∂ne=0,x∈Ω,where x′=x′,y′,z′∈Ω and x=x,y,z, δ is the Dirac function, and ne is the vector normal to outer boundary Ω. Equation (1b) is the absorbing boundary at the void-solid interface Γ. Equation (1c) is the reflection (no-flux) boundary at the departure region Ω (outer boundary).One may notice that (1) is actually a Green’s function problem. Its solution cx∣x′ is the Green’s function of the following Poisson equation:
(2a)Δux=−fx,(2b)DirichletboundaryatΓ,(2c)NeumannboundaryatΩ.For some simple geometries, the analytical expressions ofcx∣x′ can be found or easily derived from the libraries of Green’s functions in the literature (e.g., Cole et al. [22]). When cx∣x′ is known, the probability distribution of the hit point x∗ at the inner boundary can be expressed as (see, e.g., Deaconu and Lejay [23])
(3)px∗∣x0=x′=−∂cx∗∣x′∂n,x∗∈Γ,where n is the vector normal to Γ and pointing to solid. Hence, the penetration depth and MPD can be calculated via cx∣x′. That is, the penetration depth of point x′ is lex′=lx′,ω=−∫Γz∂cx∣x′/∂n∣x∈Γdx, which will be further averaged on the departure region Ω to calculate the MPD, i.e.,
(4)DG=1S∫x′∈Ωlex′dx′=−1S∫x′∈Ωdx′∫Γz∂cx∣x′∂n∣x∈Γdx,where S is the area of Ω.We consider the pores of three simple geometries for which Green’s functioncx∣x′ is available, i.e., parallel planes, straight square, and circular tubes (Figures 2(a)–2(c)). Here, we used the “isoperimetric ratio” φ as a measure of the pore cross-sectional shape. φ=lc2/A is used to measure how greatly the shape differs from a circle, where lc is the length of a curve (perimeter) and A is the corresponding area. φ has a minimum value 4π for the circle and higher values for all other shapes [24]. For parallel planes φ=∞, since the length of planes is infinite and the plane distance d is finite; for the square cross section, φ=4d2/d2=16; for the equilateral triangle, φ=3d2/3d2/4=123≈20.8; for the ellipse with semiaxes a and ba≥b, 4π≤φ<∞ and φ≈15.0 as a=2b.Figure 2
Void spaces (a) between parallel planes and in straight tubes of (b) a square cross section; (c) a circular cross section; (d) an elliptic cross section; and (e) an equilateral triangular cross section. (f) An ideal porous medium sample consisting of solid and many straight pores. The green area indicates the outer boundary.
(a)
(b)
(c)
(d)
(e)
(f)For parallel planes of distanced (Figure 2(a)), the corresponding Green’s function is expressed as [22]
(5)cx,z∣x′,z′=0=2π∑n=1∞1ne−nπz/dsinnπxdsinnπx′d.Applying Equation (4), we obtain the MPD based on Green’s function, i.e.,
(6)DG=7ζ3dπ3≈0.2714d,where ζ∙ is the zeta function and the subscript “G” stands for Green’s function. For straight tubes of square cross sections (with side-length d, Figure 2(b)), Green’s function is
(7)cx,y,z∣x′,y′,z′=0=4πd∑m=1∞∑n=1∞exp−πn2+m2z/dn2+m2sinmπxdsinmπx′dsinnπydsinnπy′d,and the MPD based on Green’s function is
(8)DG=16dπ5∑m=1∞∑n=1∞1n2n2+m23/21−−1m1−−1n2≈16dπ5⋅3.2535=0.1701d.For straight tubes of circular cross sections (with radiusR, Figure 2(c)), Green’s function is written as
(9)cr,z∣r′,z′=0=1πR∑m=1∞1βmJ12βmexp−βmzRJ0βmrRJ0βmr′R.The MPD based on Green’s function is(10)DG=4R∑m=1∞1βm3≈0.3235R,where J0 and J1 are Bessel functions (of zeroth and first orders) of the first kind; βm is the m‐th root of J0βm=0.No analytical Green’s functions are found for the remaining examples in Figure2; therefore, the numerical simulation of the Brownian motion will be used to directly estimate the MPD.Taking into account the impact of porosity, let us consider a porous medium consisting of many straight tubes (Figure2(f)), whose porosity n is not equal to 1. For simplicity, assume that all of the pores are the same. Recall that k represents the permeability of a single pore, DB is the MPD with respect to a single pore, and k=CDB2. One may note that the permeability of such a medium is
(11)kb=nk=nCDB2,where “b” stands for bunch. Let Db be the effective MPD after being averaged on the whole area S; it is obvious that Db=nDB+1−n⋅0=nDB. Let kb=CbDb2. Then, nk=Cbn2DB2, i.e., Cb=k/n/DB2. Recall (C is the factor for porosity n=1), one has Cb=C/n. Thus,
(12)kb=CDb2n,which is valid for media comprising uniform straight tubes. We will check this formula in the next section to test its applicability. The impact of tortuosity on permeability is not discussed but can be considered by some empirical equations in the literature (e.g., [13]).


## 2.2. Random Walk on Grid (WOG)

Green’s functions are accessible to only limited, simple geometries and are not practical for general pores. Thus, one has to rely on the numerical simulation of the Brownian motion or random walk. To simulate a Brownian motion in the medium, here we employ a random walk on grid (WOG) due to its ease in determining hit occurrences and its flexibility of irregular inner boundaries. In homogeneous void space, WOG is equivalent to a simple lattice random walk. Here, implementation of WOG is only briefly introduced in Figure3. In step (a), the domain can be discretized with a regular/irregular block-/point-centered grid. In this study, we use regular block-centered grids; i.e., each node represents a block or cell of the same size. In step (b), given that the distances from the current node to the neighboring nodes in six directions (i.e., x±, y±, and z±) are ∆x±, ∆y±, and ∆z± in the Cartesian coordinate and that px±, py±, and pz± denote the probability that the walker moves from the current node to six neighboring nodes, the probability of the neighboring node in x+ direction can be estimated as
(13)px+=1/Δx+1/Δx++1/Δx−+1/Δy++1/Δy−+1/Δz++1/Δz−,which is similar for other directions. If the grid is nonuniform, the probabilities are location-dependent and direction-dependent. The probabilities to neighbors at all nodes are calculated and stored before random walks are simulated. In step (c), every node x′ in the departure area will be selected to serve as the starting node in the Monte Carlo loop for NMC times (Figure 3). After completion of the Monte Carlo loop, the node will be excluded from the departure area and will not be used again in the following simulation. If the original departure area contains NS nodes after discretization, a random walk will be simulated for NS×NMC times, and the number of hit nodes is NS×NMC in total. Note that each realization of the random walk has an uncertain number of steps (in the random walk loop), and the location of its hit node is also uncertain before the realization is complete. When all of the nodes in the original departure area have served as starting nodes, the simulation procedure ends and the depths of all hit nodes will be averaged to estimate the MPD. More details of WOG algorithms and the superiority of WOG over “walk on cubes” and “walk on spheres” can be found in Nan and Wu [21].Figure 3
A flowchart of the MPD estimation via WOG.We evaluated the permeability in the pores of some simple geometries, including parallel planes and void of square, circular, elliptic, and equilateral triangular cross sections (Figures3(a)–3(e)). For the void space of these shapes, the theoretical expressions of permeability also exist (Table 1). In all of these pores, the porosities are all equal to 1. Factor C can be estimated using CB=k/DB2 in each case after DB is estimated by WOG.Table 1
Permeability, MPD, andC factor estimated by WOG and Green’s function for voids of simple geometries.
Void shape
Reference length
Theoreticalk∗
MPDDB¢
MPDDGǂ
CB¢
CGǂ
φ£
Parallel plane (Figure2(a))
Plane distanced
d212
0.2763d
0.2714d
1.091
1.131
∞
Square (Figure2(b))
Side-lengthd
0.0351d2
0.1751d
0.1701d
1.145
1.123
16
Circle (Figure2(c))
RadiusR
R28
0.3318R
0.3235R
1.135
1.194
4π
Ellipse (Figure2(d))
Semiaxes:a, b
a2b24a2+b2
0.4210b (as a=2b)
—
1.128
—
4π~∞ (15 as a=2b)
Equilateral triangle (Figure2(e))
Side-lengthd
d280
0.1054d
—
1.125
—
123
∗Theoretical k is obtained according to pressure drops and theoretical flow rates in Drazin and Norman [25]. ¢DB is the MPD estimated by random walk on grid (WOG); CB=k/DB2. ǂDG is the MPD evaluated by Green’s function cx∣x′ by Equations (6), (8), and (10); CG=k/DG2. £φ according to Berger [24].


## 2.3. COMSOL Simulation

COMSOL is a general-purpose platform software for modeling engineering problems. It is able to simulate multiple physical processes and models, including the flow of fluid. In this study, we use COMSOL to solve the steady-state Stokes flow in the pore space. The void-solid interface is treated as a no-slip boundary. The prescribed pressure boundaries are imposed onto two opposite faces of the parallelepiped and no-flux boundaries onto the other four lateral faces (due to the symmetry of the medium).


## 3. Results

### 3.1. MPD andC Factor Estimation

Table1 reports the summary results in our investigation of simple pores, including the theoretical permeability k, the MPD and factor C estimated by WOG (DB and CB), and the MPD and C estimated by the three Green’s functions (DG and CG). First, it can be seen that DB is slightly larger than DG (if available) in all three cases. This suggests that DB tends to slightly overestimate the real MPD (DG), which may lead to a limited overestimation of the permeability (6%). Such a bias may result from spatial discretization, which is further explored in the sensitivity tests below. Second, the largest C value (1.194) is only 9% larger than the smallest C value (1.091). That is, the estimated values of C are more or less the same. The average value for five cross-sectional shapes is approximately 1.125. Figure 4 shows the C value versus the isoperimetric ratio φ. The smaller φ is, the closer the shape is to being circular. It seems that there is no remarkable relationship between φ and C, which means that the dependence of C on the roundness of the cross section is quite weak.Figure 4
The isoperimetric rationφ versus factor C (a=2b for ellipse case).In Hwang et al. [5], C was simply assumed to be 1, leading to the relationship k=D2. However, our results from both numerical simulation and Green’s function, as shown above, indicate that the factor C≈1.125 on average rather than 1. Thus, the previous relation k=D2 underestimates the permeability by 10% or so compared to that of k=1.125D2.To investigate the sensitivity of the MPD estimates to the ensemble size,NMC in Monte Carlo experiments, and grid discretization length Δ, we estimate the values of MPD for parallel planes using five different values of Δ (10−3 mm, 3 × 10−3 mm, 10−2 mm, 3 × 10−2 mm, and 10−1 mm) and nine different NMC values (102, 3 × 102, 103, 3 × 103, 104, 3 × 104, 105, 3 × 105, and 106). NMC=104 when Δ is changed. Furthermore, Δ=10−2mm when NMC is changed. The resulting errors (e=DB−DG) of MPD are reported in Figure 5. In Figure 5(a), it is found that the error e has the same order of magnitude as that of Δ (i.e., e∝O∆) and continues decreasing as the grid becomes finer. Figure 5(b) illustrates the two-regime behavior of the error e with respect to NMC. That is, when NMC<104, the error decreases with the increase in NMC; when NMC>104, the error stabilizes at approximately 0.0076 (close to the error of Δ=10−2mm, NMC=104), which means that when Δ=10−2mm, the increase of realizations will not improve the accuracy of estimates remarkably once NMC>104. The discretization error is dominant compared to the sampling error in Monte Carlo simulation. The MPD estimates seem to be sensitive to the discretization length.Figure 5
Sensitivity tests of the MPD estimates to discretization the length Δ and ensemble sizeNMC.
(a)
NMC=104
(b)
Δ=10−2mm


### 3.2. Permeability Validation in a More Complex Medium

In this section, we predict the permeability of a given porous medium using a WOG-based MPD method. Meanwhile, COMSOL is used to solve the steady-state Stokes flow in the pore and estimate the permeability using the dynamics method to verify the MPD method.Although there are geometric models available for isotropic random porous media (see, e.g., [26]), a synthetic ideal porous medium (Figure 6) is used here so that we can reconstruct the void space precisely in both the WOG and COMSOL simulations and minimize the model error. The void space is generated by removing tightly placed spheres of same size from a cuboid. To reduce the computational cost, we cut a representative volume (red-edge prism in Figure 6) from the medium, which contained a pore with enough length (Figure 7). The radius of each sphere is d/2. d is sufficiently small so that the requirement of Stokes flow can be easily satisfied in the following dynamic simulation. In the COMSOL simulation, we let d=1mm. The porosity of the medium was n=8−4π/3/8=0.476.Figure 6
A synthetic porous medium consisting of solid spheres of equal radius (1 mm) and void.Figure 7
A void space formed among equal-sized solid spheres. The top face (red) is the outlet, and the bottom face (hidden) is the inlet; all the lateral planes (blue) are no-flow boundaries; the spherical faces (gray) are interfaces between the void and solid and are treated as a no-slip boundary in dynamic simulation.For simulating the Brownian motion with WOG, we established a grid of size100×100×300 in the space. The number of Monte Carlo realizations was NMC=104 for each starting point. Figure 8 shows the locations where the Brownian motion stops (note that the Brownian motion starts at the bottom surface). The starting points within the solid immediately lead to stopping and, thus, are excluded. The estimated MPD is DB=0.07591d. Here, we use the average value of C=1.125. The predicted permeability by MPD is kb=nCDB2=0.476×1.125×0.075912d2≈3.09×10−3d2.Figure 8
Locations where the simulated Brownian motions stop, i.e., points on the surfaces of the void-solid interface on the100×100×300 grid. The colors from blue to green-yellow indicate the increasing z values.In a dynamic simulation via COMSOL software, we letd=1mm, the length in z direction Lz=3mm, the pressure difference ∆P=0.1Pa between the top and bottom faces, the dynamic viscosity of water μ=1.002×10−3Pa·s at 20°C, and the sample cross-sectional area S=d2. The void-solid interface is set to be a no-slip boundary (gray faces in Figure 7). The blue lateral faces are cut sections of the void, which are treated as a no-flux boundary due to the void symmetry within the entire medium.The velocity field is briefly shown in Figure9, with two slices at planes of x=0.5 and z=0.5. The maximum velocity is approximately 8.4 × 10−5 m/s, and the Reynolds number is 0.093, which is much smaller than 1 (satisfying the requirement of Stokes flow). For the simulated steady-state flow, we found that the flowrate at the outlet was Q=8.3552×10−11m3/s. Since Q=kb∆P/μLzS, then kb=Q/SμLz/∆P=2.52×10−9m2=2.52×10−3d2. Recall that the permeability estimated by MPD, as calculated above, is kb=3.09×10−3d2, which is approximately 23% larger than that estimated by the COMSOL simulation and means that the MPD method here also overestimates the permeability by 1/4.Figure 9
The magnitude of velocity on planes ofx=0.5mm and z=0.5mm solved by COMSOL.To explore the impact of Reynolds number to permeability estimated by COMSOL, we calculated permeability estimates using different pressure drops (∆P=10−4,10−3, 10−2, 10−1, 1, 10, and 100 Pa). The results are shown in Figure 10. It was found that when ∆P<1Pa (Reynolds number<0.94), the permeability estimate stabilized and was insensitive to the Reynolds number.Figure 10
The permeability (normalized byd2) for the different Reynolds numbers. The corresponding values of pressure drop ∆P are also shown in red (in Pa).


## 3.1. MPD andC Factor Estimation

Table1 reports the summary results in our investigation of simple pores, including the theoretical permeability k, the MPD and factor C estimated by WOG (DB and CB), and the MPD and C estimated by the three Green’s functions (DG and CG). First, it can be seen that DB is slightly larger than DG (if available) in all three cases. This suggests that DB tends to slightly overestimate the real MPD (DG), which may lead to a limited overestimation of the permeability (6%). Such a bias may result from spatial discretization, which is further explored in the sensitivity tests below. Second, the largest C value (1.194) is only 9% larger than the smallest C value (1.091). That is, the estimated values of C are more or less the same. The average value for five cross-sectional shapes is approximately 1.125. Figure 4 shows the C value versus the isoperimetric ratio φ. The smaller φ is, the closer the shape is to being circular. It seems that there is no remarkable relationship between φ and C, which means that the dependence of C on the roundness of the cross section is quite weak.Figure 4
The isoperimetric rationφ versus factor C (a=2b for ellipse case).In Hwang et al. [5], C was simply assumed to be 1, leading to the relationship k=D2. However, our results from both numerical simulation and Green’s function, as shown above, indicate that the factor C≈1.125 on average rather than 1. Thus, the previous relation k=D2 underestimates the permeability by 10% or so compared to that of k=1.125D2.To investigate the sensitivity of the MPD estimates to the ensemble size,NMC in Monte Carlo experiments, and grid discretization length Δ, we estimate the values of MPD for parallel planes using five different values of Δ (10−3 mm, 3 × 10−3 mm, 10−2 mm, 3 × 10−2 mm, and 10−1 mm) and nine different NMC values (102, 3 × 102, 103, 3 × 103, 104, 3 × 104, 105, 3 × 105, and 106). NMC=104 when Δ is changed. Furthermore, Δ=10−2mm when NMC is changed. The resulting errors (e=DB−DG) of MPD are reported in Figure 5. In Figure 5(a), it is found that the error e has the same order of magnitude as that of Δ (i.e., e∝O∆) and continues decreasing as the grid becomes finer. Figure 5(b) illustrates the two-regime behavior of the error e with respect to NMC. That is, when NMC<104, the error decreases with the increase in NMC; when NMC>104, the error stabilizes at approximately 0.0076 (close to the error of Δ=10−2mm, NMC=104), which means that when Δ=10−2mm, the increase of realizations will not improve the accuracy of estimates remarkably once NMC>104. The discretization error is dominant compared to the sampling error in Monte Carlo simulation. The MPD estimates seem to be sensitive to the discretization length.Figure 5
Sensitivity tests of the MPD estimates to discretization the length Δ and ensemble sizeNMC.
(a)
NMC=104
(b)
Δ=10−2mm


## 3.2. Permeability Validation in a More Complex Medium

In this section, we predict the permeability of a given porous medium using a WOG-based MPD method. Meanwhile, COMSOL is used to solve the steady-state Stokes flow in the pore and estimate the permeability using the dynamics method to verify the MPD method.Although there are geometric models available for isotropic random porous media (see, e.g., [26]), a synthetic ideal porous medium (Figure 6) is used here so that we can reconstruct the void space precisely in both the WOG and COMSOL simulations and minimize the model error. The void space is generated by removing tightly placed spheres of same size from a cuboid. To reduce the computational cost, we cut a representative volume (red-edge prism in Figure 6) from the medium, which contained a pore with enough length (Figure 7). The radius of each sphere is d/2. d is sufficiently small so that the requirement of Stokes flow can be easily satisfied in the following dynamic simulation. In the COMSOL simulation, we let d=1mm. The porosity of the medium was n=8−4π/3/8=0.476.Figure 6
A synthetic porous medium consisting of solid spheres of equal radius (1 mm) and void.Figure 7
A void space formed among equal-sized solid spheres. The top face (red) is the outlet, and the bottom face (hidden) is the inlet; all the lateral planes (blue) are no-flow boundaries; the spherical faces (gray) are interfaces between the void and solid and are treated as a no-slip boundary in dynamic simulation.For simulating the Brownian motion with WOG, we established a grid of size100×100×300 in the space. The number of Monte Carlo realizations was NMC=104 for each starting point. Figure 8 shows the locations where the Brownian motion stops (note that the Brownian motion starts at the bottom surface). The starting points within the solid immediately lead to stopping and, thus, are excluded. The estimated MPD is DB=0.07591d. Here, we use the average value of C=1.125. The predicted permeability by MPD is kb=nCDB2=0.476×1.125×0.075912d2≈3.09×10−3d2.Figure 8
Locations where the simulated Brownian motions stop, i.e., points on the surfaces of the void-solid interface on the100×100×300 grid. The colors from blue to green-yellow indicate the increasing z values.In a dynamic simulation via COMSOL software, we letd=1mm, the length in z direction Lz=3mm, the pressure difference ∆P=0.1Pa between the top and bottom faces, the dynamic viscosity of water μ=1.002×10−3Pa·s at 20°C, and the sample cross-sectional area S=d2. The void-solid interface is set to be a no-slip boundary (gray faces in Figure 7). The blue lateral faces are cut sections of the void, which are treated as a no-flux boundary due to the void symmetry within the entire medium.The velocity field is briefly shown in Figure9, with two slices at planes of x=0.5 and z=0.5. The maximum velocity is approximately 8.4 × 10−5 m/s, and the Reynolds number is 0.093, which is much smaller than 1 (satisfying the requirement of Stokes flow). For the simulated steady-state flow, we found that the flowrate at the outlet was Q=8.3552×10−11m3/s. Since Q=kb∆P/μLzS, then kb=Q/SμLz/∆P=2.52×10−9m2=2.52×10−3d2. Recall that the permeability estimated by MPD, as calculated above, is kb=3.09×10−3d2, which is approximately 23% larger than that estimated by the COMSOL simulation and means that the MPD method here also overestimates the permeability by 1/4.Figure 9
The magnitude of velocity on planes ofx=0.5mm and z=0.5mm solved by COMSOL.To explore the impact of Reynolds number to permeability estimated by COMSOL, we calculated permeability estimates using different pressure drops (∆P=10−4,10−3, 10−2, 10−1, 1, 10, and 100 Pa). The results are shown in Figure 10. It was found that when ∆P<1Pa (Reynolds number<0.94), the permeability estimate stabilized and was insensitive to the Reynolds number.Figure 10
The permeability (normalized byd2) for the different Reynolds numbers. The corresponding values of pressure drop ∆P are also shown in red (in Pa).


## 4. Discussion and Conclusions

The effectiveness of the WOG-based MPD method was tested in this study via two approaches, i.e., comparing the MPD to Green’s function results and comparing the permeability to the results of dynamic simulation. First, the results of the MPD estimated by WOG are quite close to those based on Green’s functions in three cases (Figures2(a)–2(c)) and its relative errors were less than 5%. We found that WOG estimates the MPD with good accuracy when compared to the theoretical MPD in the simple pores (although WOG still overestimates the MPD slightly, <6%). Second, the permeability of the synthetic porous medium as estimated by the WOG-based MPD method was compared with that estimated by the COMSOL simulation. To the best of our knowledge, this is the first time that an MPD estimate of the permeability was compared with the permeability estimated by a dynamic model in a “complex” porous medium. Our results suggested that the MPD-based estimate of permeability was approximately 1/4 larger than the dynamics-based estimate, which contrasts the accuracy in the pores of the simple geometries above. The deviation may result from a dramatic change in the void cross-sectional area and in the singularities generated by the tangent spheres.While factorC was taken as 1 in the literature, we found that in pores of several geometries, the factor C≈1.125 on average using both the WOG simulation and Green’s function. Factor C depends on the cross-sectional shape very weakly. The largest C value (1.194) is only 9% larger than the smallest C value (1.091). After analyzing the plot of the C value against a shape parameter φ (isoperimetric ratio), it was found that there is no remarkable relationship between φ and C, which means that the dependence of C on the cross-sectional shape is quite weak.In sensitivity tests of MPD to discretization of the length Δ and ensemble sizeNMC, it was found that the error of MPD estimates has the same order of magnitude as that of Δ. If Δ>10−3mm, the sampling error is negligible as long as NMC>104.When the pressure drop is sufficiently small (∆P<1Pa, Reynolds number<0.94), the flow in the dynamic simulation is Stokes flow, and the permeability estimate stabilizes and is insensitive to the Reynolds number.In summary, the effectiveness of the permeability estimation using the WOG-based MPD was validated using Green’s functions in simple pores and dynamic simulation in a complex pore. It was found that the WOG-based MPD method is quite accurate in simple pores and overestimates the permeability in the complex pore. It was shown that the dependence of the factorC on cross-sectional geometry is very limited, and the average C value can be used. The discretization error and the Monte Carlo sampling error both result in errors in MPD estimation.



---

*Source: 9240203-2019-01-08.xml*