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

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 [38], 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 ([912]; 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 by combining three equations: (Hagen-Poiseuille law), (Darcy’s law), and (relation between permeability and conductivity) (see, e.g., [13]). For tubes of elliptical cross sections, , where and are the semiaxes [14]; similar results are available for several other shapes [15, 16]. For media of higher complexity, for example, in a packed bed of solids, where is the mean grain diameter, is the porosity, and 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 ), 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 . Simonov and Mascagni [7] took porosity into account and suggested that , where is the porosity and is a factor. They estimated by using the Poiseuille law-based permeability in an ideal, straight pore of the square cross section with a side length that is equal to , and . It is still not clear whether 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., , where is some type of characteristic length (e.g., hydraulic radius, mean pore size, mean grain diameter, or MPD), and the coefficient may depend on the cross section and/or porosity ([13]; and references therein).

If one uses MPD as the characteristic length, the coefficient may depend on both cross-section shape and porosity. The influence of porosity on 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 ’s dependence on the roundness of pore cross sections, we derived the theoretical values of MPD and 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 on the cross-sectional geometry of void space, but they did not investigate the problem further. Here, the dependence is investigated through derived by Green’s functions, as well as the 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 is calibrated for each geometry. Furthermore, the average 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 at (denoted by ) on its outer boundary plane (say, ; 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, . The depth of the hit point relative to the outer boundary plane (i.e., ; is the ensemble index) is called the penetration depth of the Brownian motion, denoted by . 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., , where stands for expectation with respect to the ensemble index , denotes the MPD estimated by the Brownian motion (or random walk), and 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 , the spatial density (i.e., mean times of occurrence in space) of the diffusion satisfies the following equation: where and , is the Dirac function, and 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 is the Green’s function of the following Poisson equation:

For some simple geometries, the analytical expressions of can be found or easily derived from the libraries of Green’s functions in the literature (e.g., Cole et al. [22]). When is known, the probability distribution of the hit point at the inner boundary can be expressed as (see, e.g., Deaconu and Lejay [23]) where is the vector normal to and pointing to solid. Hence, the penetration depth and MPD can be calculated via . That is, the penetration depth of point is , which will be further averaged on the departure region to calculate the MPD, i.e., where is the area of .

We consider the pores of three simple geometries for which Green’s function 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. is used to measure how greatly the shape differs from a circle, where is the length of a curve (perimeter) and 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 is finite; for the square cross section, ; for the equilateral triangle, ; for the ellipse with semiaxes and , and as .

For parallel planes of distance (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 , Figure 2(b)), Green’s function is and the MPD based on Green’s function is

For straight tubes of circular cross sections (with radius , Figure 2(c)), Green’s function is written as

The MPD based on Green’s function is where and are Bessel functions (of zeroth and first orders) of the first kind; is the root of .

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 is not equal to 1. For simplicity, assume that all of the pores are the same. Recall that represents the permeability of a single pore, is the MPD with respect to a single pore, and . One may note that the permeability of such a medium is where “b” stands for bunch. Let be the effective MPD after being averaged on the whole area ; it is obvious that . Let . Then, , i.e., . Recall ( is the factor for porosity ), one has . 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]).

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 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., , , and ) are , , and in the Cartesian coordinate and that , , and denote the probability that the walker moves from the current node to six neighboring nodes, the probability of the neighboring node in direction can be estimated as 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 in the departure area will be selected to serve as the starting node in the Monte Carlo loop for 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 nodes after discretization, a random walk will be simulated for times, and the number of hit nodes is 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].

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 can be estimated using in each case after is estimated by WOG.

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 and Factor Estimation

Table 1 reports the summary results in our investigation of simple pores, including the theoretical permeability , the MPD and factor estimated by WOG ( and ), and the MPD and estimated by the three Green’s functions ( and ). First, it can be seen that is slightly larger than (if available) in all three cases. This suggests that tends to slightly overestimate the real MPD (), 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 value (1.194) is only 9% larger than the smallest value (1.091). That is, the estimated values of are more or less the same. The average value for five cross-sectional shapes is approximately 1.125. Figure 4 shows the 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 , which means that the dependence of on the roundness of the cross section is quite weak.

In Hwang et al. [5], was simply assumed to be 1, leading to the relationship . However, our results from both numerical simulation and Green’s function, as shown above, indicate that the factor on average rather than 1. Thus, the previous relation underestimates the permeability by 10% or so compared to that of .

To investigate the sensitivity of the MPD estimates to the ensemble size, 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 values (102, 3 × 102, 103, 3 × 103, 104, 3 × 104, 105, 3 × 105, and 106). when Δ is changed. Furthermore, when is changed. The resulting errors () of MPD are reported in Figure 5. In Figure 5(a), it is found that the error has the same order of magnitude as that of Δ (i.e., ) and continues decreasing as the grid becomes finer. Figure 5(b) illustrates the two-regime behavior of the error with respect to . That is, when , the error decreases with the increase in ; when , the error stabilizes at approximately 0.0076 (close to the error of , ), which means that when , the increase of realizations will not improve the accuracy of estimates remarkably once . 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.

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 . 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 . The porosity of the medium was .

For simulating the Brownian motion with WOG, we established a grid of size in the space. The number of Monte Carlo realizations was 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 . Here, we use the average value of . The predicted permeability by MPD is .

In a dynamic simulation via COMSOL software, we let , the length in direction , the pressure difference between the top and bottom faces, the dynamic viscosity of water at 20°C, and the sample cross-sectional area . 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 Figure 9, with two slices at planes of and . 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 . Since , then . Recall that the permeability estimated by MPD, as calculated above, is , 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 (, 10−2, 10−1, 1, 10, and 100 Pa). The results are shown in Figure 10. It was found that when (Reynolds ), the permeability estimate stabilized and was insensitive to the Reynolds number.

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 (Figures 2(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 factor was taken as 1 in the literature, we found that in pores of several geometries, the factor on average using both the WOG simulation and Green’s function. Factor depends on the cross-sectional shape very weakly. The largest value (1.194) is only 9% larger than the smallest value (1.091). After analyzing the plot of the value against a shape parameter (isoperimetric ratio), it was found that there is no remarkable relationship between and , which means that the dependence of on the cross-sectional shape is quite weak.

In sensitivity tests of MPD to discretization of the length Δ and ensemble size , it was found that the error of MPD estimates has the same order of magnitude as that of Δ. If , the sampling error is negligible as long as .

When the pressure drop is sufficiently small (, Reynolds ), 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 on cross-sectional geometry is very limited, and the average 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.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work is funded by the National Key Research and Development Program of China (Grant No. 2016YFC040281000), the National Natural Science Foundation of China (Grant No. 41602250), and the Fundamental Research Funds for the Central Universities (Grant No. 0206-14380032).