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

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 factor C 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.


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][4][5][6][7][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][10][11][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 obtains k = R 2 /8 by combining three equations: Q = πR 4 /8μ ΔP/L (Hagen-Poiseuille law), Q = K ΔP/ρgL πR 2 (Darcy's law), and K = k ρg/μ (relation between permeability and conductivity) (see, e.g., [13]). For tubes of elliptical cross sections, k = a 2 b 2 / 4 a 2 + b 2 , 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 = c K n 3 / 1 − n 2 d 2 in a packed bed of solids, where d is the mean grain diameter, n is the porosity, and c K 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 = D 2 . Simonov and Mascagni [7] took porosity into account and suggested that k = C SM nD 2 , where n is the porosity and c SM is a factor. They estimated C SM = k/nD 2 ≈ 6 63 by using the Poiseuille lawbased permeability k = 0 0351d 2 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 = C SM nD 2 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 = CL 2 , 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 coefficient C 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 of C'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 crosssectional 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 Section 2, 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.

Methodology
2.1. Mean Penetration Depth (MPD) and MPD Calculation via Green's Function. Given a porous medium, let a Brownian motion initiate at x′ (denoted by x 0 = 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 l x′, ω . 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., D B = 1/S x′∈Ω l x′, ω dx′ , where stands for expectation with respect to the ensemble index ω, D B denotes the MPD estimated by the Brownian motion (or random walk), and S is the area of Ω.
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 point x 0 = x ′ , the spatial density c x | x ′ 2 Geofluids (i.e., mean times of occurrence in space) of the diffusion satisfies the following equation: where x ′ = x ′ , y ′ , z ′ ∈ Ω and x = x, y, z , δ is the Dirac function, and n e 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 c x | x ′ is the Green's function of the following Poisson equation: Dirichlet boundary at Γ, 2b Neumann boundary at Ω 2c For some simple geometries, the analytical expressions of c x | x ′ can be found or easily derived from the libraries of Green's functions in the literature (e.g., Cole et al. [22]). When c x | 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]) where n is the vector normal to Γ and pointing to solid. Hence, the penetration depth and MPD can be calculated via c x | x′ . That is, the penetration depth of point x′ is l e x′ = l x′, ω = − Γ z ∂c x | x′ /∂n | x∈Γ dx, which will be further averaged on the departure region Ω to calculate the MPD, i.e., where S is the area of Ω.
We consider the pores of three simple geometries for which Green's function c x | 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. φ = l 2 c /A is used to measure how greatly the shape differs from a circle, where l c 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; 3 Geofluids for the square cross section, φ = 4d 2 /d 2 = 16; for the equilateral triangle, φ = 3d 2 / 3d 2 /4 = 12 3 ≈ 20 8; for the ellipse with semiaxes a and b a ≥ b , 4π ≤ φ < ∞ and φ ≈ 15 0 as a = 2b.
For parallel planes of distance d (Figure 2(a)), the corresponding Green's function is expressed as [22] Applying Equation (4), we obtain the MPD based on Green's function, i.e., 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 and the MPD based on Green's function is For straight tubes of circular cross sections (with radius R, Figure 2(c)), Green's function is written as The MPD based on Green's function is where J 0 and J 1 are Bessel functions (of zeroth and first orders) of the first kind; β m is the m-th root of J 0 β m = 0.
No analytical Green's functions are found for the remaining examples in Figure 2; 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 (Figure 2(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, D B is the MPD with respect to a single pore, and k = CD 2 B . One may note that the permeability of such a medium is where "b" stands for bunch. Let D b be the effective MPD after being averaged on the whole area S; it is obvious that Recall (C is the factor for porosity n = 1), one has C b = C/n. Thus, 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]).

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 Figure 3. 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 p x± , p y± , and p z± 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 p x+ = 1/Δx + 1/Δx + + 1/Δx − + 1/Δy + + 1/Δy − + 1/Δz + + 1/Δz − , 13 which is similar for other directions. If the grid is nonuniform, the probabilities are location-dependent and 4 Geofluids 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 N MC 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 N S nodes after discretization, a random walk will be simulated for N S × N MC times, and the number of hit nodes is N S × N MC in total. Note that each realization of the random walk has an uncertain number of steps (in the random walk loop),  5 Geofluids 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].
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 (Figures 3(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 C B = k/D 2 B in each case after D B is estimated by WOG.

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). Table 1 reports the summary results in our investigation of simple pores, including the theoretical permeability k, the MPD and factor C estimated by WOG (D B and C B ), and the MPD and C estimated by the three Green's functions (D G and C G ). First, it can be seen that D B is slightly larger than D G (if available) in all three cases. This suggests that D B tends to slightly overestimate the real MPD (D G ), 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.

MPD and C Factor Estimation.
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.
In Hwang et al. [5], C was simply assumed to be 1, leading to the relationship k = D 2 . 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 = D 2 underestimates the permeability by 10% or so compared to that of k = 1 125D 2 .

Void shape
Reference length Theoretical Parallel plane (Figure 2(a)) Plane distance d d 2 (6), (8), and (10); C G = k/D 2 G . £ φ according to Berger [24].  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 N MC . That is, when N MC < 10 4 , the error decreases with the increase in N MC ; when N MC > 10 4 , the error stabilizes at approximately 0.0076 (close to the error of Δ = 10 −2 mm, N MC = 10 4 ), which means that when Δ = 10 −2 mm, the increase of realizations will not improve the accuracy of estimates remarkably once N MC > 10 4 . 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.

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 = 1 mm. The porosity of the medium was n = 8 − 4π/3 /8 = 0 476.
For simulating the Brownian motion with WOG, we established a grid of size 100 × 100 × 300 in the space. The number of Monte Carlo realizations was N MC = 10 4 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 D B = 0 07591d. Here, we use the average value of C = 1 125. The predicted permeability by MPD is k b = nCD 2 B = 0 476 × 1 125 × 0 07591 2 d 2 ≈ 3 09 × 10 −3 d 2 . In a dynamic simulation via COMSOL software, we let d = 1 mm, the length in z direction L z = 3 mm, the pressure difference ΔP = 0 1 Pa between the top and bottom faces, the dynamic viscosity of water μ = 1 002 × 10 −3 Pa · s at 20°C, and the sample cross-sectional area S = d 2 . The voidsolid interface is set to be a no-slip boundary (gray faces in  The velocity field is briefly shown in Figure 9, 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 −11 m 3 /s. Since Q = k b ΔP/μL z S, then k b = Q/S μL z / ΔP = 2 52 × 10 −9 m 2 = 2 52 × 10 −3 d 2 . Recall that the permeability estimated by MPD, as calculated above, is k b = 3 09 × 10 −3 d 2 , 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.
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 < 1 Pa (Reynolds number < 0 94), the permeability estimate stabilized and was insensitive to the Reynolds number.

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 (Figures 2(a)-2(c)) and its 8 Geofluids 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 factor C 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 size N MC , it was found that the error of MPD estimates has the same order of magnitude as that of Δ. If Δ > 10 −3 mm, the sampling error is negligible as long as N MC > 10 4 .
When the pressure drop is sufficiently small (ΔP < 1 Pa, 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 factor C 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.

Data Availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.