Comparative analysis on pore‐scale permeability prediction on micro‐CT images of rock using numerical and empirical approaches

Varieties of pore‐scale numerical and empirical approaches have been proposed to predict the rock permeability when the pore structure is known, for example, microscopic computerized tomography (micro‐CT) technology. A comparative study on these approaches is conducted in this paper. A reference dataset of nine micro‐CT images of porous rocks is generated and processed including artificial sandpacks, tight sandstone, and carbonate. Multiple numerical and empirical approaches are used to compute the absolute permeability of micro‐CT images including the image voxel‐based solver (VBS), pore network model (PNM), Lattice Boltzmann method (LBM), Kozeny‐Carman (K‐C) equation, and Thomeer relation. Computational accuracy and efficiency of different numerical approaches are investigated. The results indicate that good agreements among numerical solvers are achieved for the sample with a homogeneous structure, while the disagreement increases with an increase in heterogeneity and complexity of pore structure. The LBM and VBS solver both have a relative higher computation accuracy, whereas the PNM solver is less accurate due to simplification on the topological structure. The computation efficiency of the different solver is generally computation resources dependent, and the PNM solver is the fastest, followed by VBS and LBM solver. As expected, empirical relation can over‐estimate permeability by a magnification of 50 or more, particularly for those strong heterogeneous structures reported in this study. Nevertheless, empirical relation is still applicable for artificial rocks.

The traditional way to obtain rock physics generally relies on laboratory test on natural specimens. Due to the complex and disordered pore structure of natural rock specimens, the test results may also vary from samples under a consistent test environment. Besides, the test sample will be damaged and polluted over time, which limits the repeatable use for parallel analysis. 5,6 Thus, extraction and reconstruction of rock structure from natural rock specimens for petrophysics computation are valuable and have been regarded as an efficient approach to handle those issues. Then, digital rock petrophysics (DRP) has been adopted to obtain those macroscopic properties of rock including strength and permeability, as well as acoustic, electrical, and thermal properties, by simulating the physical process on these reconstructed models. [7][8][9][10] The most widely used DRP approaches of permeability prediction based on CT images are including pore network modeling, Lattice Boltzmann method, and image voxel-based direct N-S solver. 11 PNM, which extracts pore-throat network with topologically representative from micro-CT images, has been regarded as an effective approach for absolute permeability calculation and multi-phase flow simulation. 12,13 Generally, the Hagen-Poiseuille law is adopted to compute the conductance between the adjacent pore bodies, and the absolute permeability is derived from Darcy's law. 14 The corresponding computational methods for PNM can be classified into quasistatic and dynamic models. The marked difference lies in that the quasi-static model simulates the equilibrium states of the drainage and imbibition processes controlled by capillary force, while the dynamic model can realize the simulation of intrusion process which is time-dependent and controlled by both viscous and capillary forces. 15 Fatt 16 constructed the first network model who regarded the analogy between electricity in a random resistor network and flows in a porous medium. Afterward, PNM has been used to study the complex transport behavior in porous medium considering uniform or mixed wettability conditions, including phase change, reactive transport phenomenon, and non-Newtonian displacement. 17,18 However, PNM clearly makes approximations concerning rock pore-throats geometry which limits its application. Another tool for pore-scale transport properties prediction is direct simulation which is performed directly on the 3D-voxelized segmented image (eg, CT image). 13,19 One remarkable numerical method of this category is the LBM, which solves a discrete, mesoscale form of Boltzmann equation directly on voxel grid of segmented image, thus it is suited to solve the porous flow in complex geometries. Moreover, this method is easy to code and is ideally suited for parallel computing, whereas it is time-consuming even with a massively parallel implementation. 20 Recent researches indicate that it is possible to compute relative permeability and to compute the interfacial area in the multi-phase flow. 21,22 However, the application of relaxation factor limits the use in single-factor analysis. For more details of LBM theories and applications, the literature can be referenced. [23][24][25] The other method in this category is VBS, which solves Stokes and Naiver-Stokes equation directly on the voxel grid using finite volume method (FVM) 26 by fast Fourier transform (FFT). The VBS solver computes permeability, as well as velocity and pressure, on 3D pore space image using an adaptive grid to reduce the number of grid cells rather than a regular grid. According to previous studies, this method solves fast for high-porosity structures, but needs more iterations for low-porosity structures to reach desired accuracy. 27 The convergence speed of VBS solver, in general, depends on the complexity and heterogeneity of the porous medium, and the most challenge of this method is to simulate slip boundary conditions. Besides for the DRP approaches, the absolute permeability can also be predicted by some empirical relations which are derived from laboratory and engineering data, when the parameters of pore structure are known, for example, F I G U R E 1 CT imaging process and model reconstruction

Sample preparation
Imaging and data aquisition Reconstruction porosity, surface area, etc The most well-known is the K-C equation which is originally derived for a granular medium but is later widely applied in the geoscience. 28,29 Besides, there are also several improved empirical relations for permeability prediction. The initial form of the permeability formulas can be generally expressed as the function of pore geometry parameters and porosity, 30 which have considered more parameters to improve the accuracy and to fit specific rock types. 31,32 In this paper, nine rock samples are drilled and imaged using micro-CT technology. The absolute permeability of these samples is computed by different DRP approaches and empirical relations. The results are compared and analyzed for the accuracy and efficiency.

| Micro-CT imaging and processing
A total of nine microstructures of natural rock are extracted from micro-CT images and used as inputs for absolute permeability computation. The selected samples cover artificial sandpacks, unconsolidated sandstone, tight sandstone, and carbonate, and the porosity ranges from 10% to 31%. Considering the heterogeneity of microscopic structure of rock, a total of nine samples with different pore size distribution are adopted to make the results more reliable and representative. 33 The mini-core plugs with 3-6 mm in diameter and 7-10 mm in length are generally used for CT scanning to obtain highresolution data. The standard process of CT scanning is illustrated in Figure 1. During the scanning, the sample stage is rotated at a specific time, so that the core plug can be scanned by X-ray from a different direction. 34 Some CT datasets (S1, S5, S6, and S7) used in this study are scanned by ourselves, and the others are available from the open-access library at Imperial College London ([http://www.imper ial.ac.uk/earth-scien ce/resea rch/resea rch-group s/perm/resea rch/pore-scale-model ling/microct-images-and-netwo rks/sand-pack-f42a/] and the Digital Rocks Portal [https ://www.digit alroc kspor tal.org/proje cts/]). And the original CT images used in this paper are listed in Figure 2.

| Reconstruction of pore-scale rock model
The image processing is needed before reconstruction, including denoising, filtering, segmentation, and binarization Pore (i) P ore (j) The lattice velocity of D3Q19 model 39 of the raw images. The purpose is not only to remove the image defects, that is, noises and concentric shadows caused by system device but also to extract the target object, that is, pore space and mineral composition. 35 Then, the representative elementary volume (REV) of the rock sample is identified to optimize the computational demands. 36 The general work-flow of model reconstruction including image processing and REV selection is illustrated in Figure 3.
On the basis of extracted pore structure, pore space analysis including porosity, pore-size distribution, specific F I G U R E 6 Pore space extraction and analysis from (A) original CT image, and (B)-(D) are separated label field, extracted pore network model, and calculated pore-throat radius distribution of sample S1, respectively. The results listed subsequently are the rest of all samples. The colors of the pores are different at random in separated label field model surface area, and coordination number can be performed. These parameters are crucial to the permeability estimation using empirical relations, such as the K-C equation. The details will be given and discussed in subsequent sections.

| Theoretical basis of rock physics computation
The numerical approaches and empirical relations reported in this study are presented as follows.

| VBS solver
VBS solver solves Stokes and Navier-Stokes equations directly on a 3D-segmented voxelized image using the FVM by FFT. Combining the assumption of incompressible, Newtonian fluid within a steady-state laminar flow, the N-S equation can be simplified and given as 37 : where ⃗ ∇ and ⃗ ∇ ⋅ are the gradient operator and divergence operator;∇ 2 and μ are the Laplacian operator and dynamics viscosity; ⃗ V and P are the velocity and pressure of the fluid. The nonslip condition is adopted at the solid-fluid interface. Once the equation system above is solved, Darcy's Law is adopted to determine the absolute permeability coefficient.

| PNM solver
The PNM solver computes the absolute permeability by simulating single-phase flow in the pore network extracted from CT images. 15 The pores and throats are commonly reproduced by idealized spheres and cylinders, respectively. In the pore network model, the flow rate, Q a for single-phase flow between two connected pores i and j (as explained in Figure 4) is given by 15 : where p a i and p a j are the pressure in pores i and j, respectively, g a p,ij represents the conductance of two adjacent pores i and j of fluid a and can be derived from: where L ij is the distance from the pore-throat interface of pore i to j (also represents the total length of the pore throat, L t ), L p, i and L p, j are the radii of pore body i and j, respectively. For a given shape of the channel, the conductance g p can be derived from the Hagen-Poiseuille formula: where μ p is the dynamic viscosity of fluid a, A and G are the cross-sectional area and shape factor of the pore network model, respectively. k is a constant and for a circular, equilateral and square tube, the value is 0.6, 0.6, and 0.5623, respectively. 38 Adopting Darcy's law, the absolute permeability K a of pore network model can be derived by: where Q a is the total flow rate of fluid a through a pore network model of length L with potential pressure drop ΔP.

| LBM solver
The LBM solver solves a discrete, mesoscale Boltzmann equation and can be reduced to N-S equation in a low-Mach number limit. 25 The LBM solver used in this study is implemented by adopting the D3Q19 model, as shown in Figure 5, and Eqs. (6)-(8) constitute the iterative model of LBM. 39 The discrete direction of velocity can be expressed as: The evolution equation can be derived as follow: where f i (x,t) is the particle distribution function of lattice x at time node t in i direction; is the relaxation time, and the term f eqi which is also named equilibrium distribution function can be expressed as: Δt is the lattice velocity;Δx and Δt are lattice step and time step, respectively. And the term t which is called weight coefficient can be written as:

| Empirical relation solver
The most well-known empirical relation for permeability prediction is K-C equation, which can be written as 40 : where is the porosity of porous media (dimensionless); S represents the specific surface area which is defined as the ratio of surface area of whole pores to total volume of specimen (m −1 ); c is the Kozeny constant which depends on the geometry of porous media, for example, for cylindrical capillaries, c = 2, is the tortuosity of porous media (dimensionless). For tortuosity, the empirical relation proposed by Saxena et al 41 can be used and written as: Besides, Thomeer 42 proposed another empirical permeability model using pore size distribution and mercury intrusion capillary pressure which can be written as: where P D is the entry pressure of mercury, and G is a constant which reflects pore shape property and the impact of tortuosity (generally, G = 0.2 for siliciclastic rocks, dimensionless). The equation can also be driven in terms of pore diameter D (in μm) as: ( ± 1,0,0) (0, ± 1,0) (0,0, ± 1), i = 1, ⋯ ,6; ( ± 1, ± 1,0) ( ± 1,0, ± 1,) (0, ± 1, ± 1), i = 7, ⋯ ,18.  A, Abundant dissolution pores due to calcite dissolution (red arrows), plane-polarized light (PPL), thin section; B, micropores and inter-granular pore (red arrow) presented in "Cluster" quartz aggregates (red-dotted ellipse); C, dissolution pores associated with calcite, dissolution residue with irregular morphology (red arrows); D, secondary dissolution pore with irregular morphology observed in calcite grain (red arrow)

| Quantitative analyses of pore structure characteristics
The pore space analysis including porosity and pore size distribution is performed on the pore structure model extracted from CT images, as shown in Figure 6(A). For the purpose to reduce the impact of the exterior region on image processing and to decrease the data size, a cubic region of interest (ROI) located in the center is generally selected for model reconstruction. Then, the separated operation is performed on extracted pore space to generate a label field, as shown in Figure 6(B) which splits the whole pore space into numerous parts to generate the pore network model ( Figure 6C) and for pore size distribution analysis, as shown in Figure 6(D). The pore network extracted from the separated label field is implemented in Avizo™. 43 Compared with the pore network extracted by using the maximal ball (MB) algorithm, 15 the radius of the pore network model is computed by the equivalent volume method in Avizo, whereas the MB method is equivalent to the radius of the inscribed sphere. The main difference between these two approaches on pore-throat size computation lies in the treatment of throat. The results of pore space analysis are listed in Table 1. More details of pore type, morphology, and petrography can be seen by thin sections and scanning electron microscope (SEM) in Figure  7.

| Computed absolute permeability
The absolute permeability is calculated by simulating single-phase flow with the assuming of an incompressible, Newtonian fluid and a steady-state laminar flow.

| The first category-VBS solver
The VBS solver simulates the experimental condition by simulating fluid flow between two opposite faces while other four faces are sealed with a one-voxel-wide grid (as shown in Figure 8). In this case, the Stokes equations are solved directly on the segmented 3D images of nine rock samples. The computed permeability values using VBS solver for all nine rock samples are presented in Table 2.

| The second category-PNM solver
The extracted PNM models are shown in Figure 9, in which the spheres and cylinders represent pores and throats, respectively. The single-phase flow simulation is performed on all extracted pore networks to compute absolute permeability. In this case, the conductance between connected pores is analytically given by Hagen-Poiseuille formula, while those isolated pores do not participate in fluid flow. The absolute permeability is derived from Darcy's law, as listed in Table 3.

| The third category-LBM solver
LBM solver treats the fluid flow as a continuous motion and collision process between fluid particles with collision model. Similarity, a certain pressure drop is applied on two opposite faces to force fluid flow in and out, and the rest of four faces are sealed with a one-voxel-wide grid. In order to ensure the second-order precision, the curve boundary condition is adopted between pore and grain. 44 The computed permeability using LBM solver is listed in Table 4.
Since the LBM solver and VBS solver are both performed on image voxel grids, we visualize the normalized local velocity of both LBM and VBS simulation for further comparative analyses, as illustrated in Figure 10. According to the colormap of the velocity field, the global flow field and the main flow channel distribution of LBM (located on left side) agree well with that of VBS (right side) solver.

| The fourth categoryempirical relation
In this study, the empirical K-C equation and another empirical relation proposed by Thomeer are used to predict the absolute permeability of all nine rock samples. It is reported that the empirical relation of permeability is the least accurate approach compared with other solvers, which has been confirmed in previous studies. 30,41 The empirical relation is designed for comparisons in this study. The estimated permeability is listed in Table 5. Compared with simulation results, both of the empirical relations over-estimate the permeability of rock samples, especially the results of Thomeer relation.

| Comparative analysis
A comparative analysis is carried out to present the strategies on solver selection in DRP analyses, including computed deviation and computation time, as well as the solver variability aiming at different types of rock.
There is no laboratory-measured permeability for natural rocks reported in this study, which means the direct comparison of computation accuracy between various numerical solvers is impossible. The reasons lie in the environmental difference and scaling effect between experimental test and simulation, as well as the solid-fluid interaction which is generally ignored in single-phase flow simulations. Thus, it is meaningful for this study to focus on the comparison among various numerical approaches and empirical relations.
As is shown in Figure 11, the computed permeability values of three numerical solvers are basically in a great agreement with each other whereas the estimated values of empirical relations are relatively discretized and higher. For samples with strong sensitivity and heterogeneity, such as tight sandstone and carbonate, the estimated permeabilities by using empirical relations are less accurate. Compared with Thomeer relation, the K-C equation comprehensively considers the specific surface area and tortuosity of pore structure, and the values estimated by using the K-C equation are closer to the simulation results. For sample S2 (an artificial sandpacks with a resolution at 9.996 μm) in Figure 11, the digitally computed permeability values of three different solvers are in great agreement with K-C equation value. Considering that the K-C equation is originally derived for granular medium using the laboratory-tested results, it can be regarded as a validation for numerical solvers. However, the Thomeer relation still over-estimates the permeability for sandpacks.
Using the mean value of computed permeability as references, the ratio of absolute permeability computed using different numerical solver to the references is presented in Figure 12, because the permeability computed by different approaches agrees well with each other. It can be found that the PNM solver shows the relative larger difference with respect to the reference value (mean value), whereas the relative deviations of VBS solver are less than 0.1 except for rock S8 and S9. The results presented in Figure 12 demonstrate that the VBS solver shows the best agreement with the reference value, followed by LBM and PNM solver. It should be noted the computed values of VBS and LBM solver are close to each other. Only the values computed by PNM solver are obviously discrete. Especially for carbonate, the rocks with strong heterogeneous structure, the disagreements are large over 0.5 and a maximum value of 0.66. It can be concluded by Figures 11 and 12 that the VBS and LBM solver show relatively higher computation accuracy, whereas the PNM solver is less accurate. However, the PNM solver costs a minimum computation time, followed by VBS and LBM solvers.
Computation efficiency is also of significant interest in the image-based rock physics simulation. The numerical solvers described in this study are distinct in their governing equations, as well as the implementation algorithms. Thus, we implemented the solving process on the same platform to quantitatively analyze the time-consuming of different solvers. The solving processes of different solvers on nine rock samples reported in this study are all implemented on PCs equipped with Core i7-8700, i7-2600K and 32G RAM. The comparisons on the adopted platform, running time, boundary condition setups, and convergence criteria for various solvers are listed in detail in Table 6.

F I G U R E 9
The extracted pore network model from sample S9.
In this case, gray represents the grain region, and the colorful spheres and cylinders represent pores and throats respectively, where the color varies from the size of spheres and cylinders, which is determined by the normalized radius of pore and throat Overall, the computation time is generally depended on the model structure and size. The PNM solver is relatively faster compared with VBS and LBM solver but also relatively less accurate. In addition, the extra time spent on pore network extraction process is not considered. The VBS solver generally runs faster on GPU card compared with CPU computation, especially for the model with complex structure. The LBM solver owes the advantages of large-scale parallel implementation, whereas is the slowest solver for the given computation platform reported in this study.

| CONCLUSIONS
Image-based rock permeability computation by using numerical approaches contributes to reducing the subsurface uncertainties. Compared with laboratory-test-based measurements, the computing efficiency and accuracy of the selected numerical solver are of great concern in the DRP technique.
Thus the performance of efficiency and convergence of these numerical solver become relevant.
The work described in this study presents a comprehensive comparison between multiple numerical solvers and empirical relations on computation efficiency and accuracy of permeability. It can be concluded that the results of the three different numerical solvers agree well with each other, whereas the Constant pressure at the inlet and outlet, no-slip boundary condition at the side wall and curve boundary condition is adopted between pore and grain, this model implementation with single-relaxation time (SRT set to 1.5 in this study) 4/>24 (Depending on the processers called in the calculation) empirical models generally over-estimate the permeability by a magnification of 50 or more, especially for samples with complex and heterogeneous structure. The PNM solver needs the shortest computation time, whereas the LBM solver requires the longest running time. The PNM has a relatively less accuracy compared with the other two approaches due to the simplification on topological structure, especially for samples with strong heterogeneous structure, that is, carbonate.