Reliable estimation of hydraulic permeability from 3D X-ray CT images of porous rock

The hydraulic permeability is a key parameter for simulating the flow-related phenomenon so that its accurate estimation is crucial in both experimental and numerical simulation studies. 3D pore structure can be readily taken by X-ray computed tomography (CT) and it often serves as a flow domain for pore-scale simulation. However, one encounters the challenges in segmenting the authentic pore structure owing to the finite size of image resolution and segmentation methods. Therefore, the loss of structural information in pore space seems unavoidable to result in the unreliable estimation of permeability. In this study, we propose a novel framework to overcome these limitations by using a flexible ternary segmentation scheme. Given the pore size distribution curve and porosity, three phases of pore, solid, and gray regions are segmented by considering the partial volume effect which holds the composition information of unresolved objects. The resolved objects such as solid and pore phases are taken to equivalently solve Stokes equation while the fluid flow through unresolved objects is simultaneously solved by Stokes-Brinkmann equation. The proposed numerical scheme to obtain the permeability is applied to Indiana limestone and Navajo sandstone. The results show that the computed hydraulic permeability is similar to the experimentally obtained value without being affected by image resolution. This approach has advantages of achieving consistent permeability values, less influenced by segmentation methods.


Introduction
3D X-ray μ-CT images are often used as input domain for direct numerical simulations. Conventional approaches prepare the domain as binary pore structure using segmentation schemes such as Otsu's method, watershed method, visual inspection decision, etc. However, these approaches lead to the unavoidable loss of spatial information embedded in X-ray CT images because they neglect the partial volume effect of X-ray CT.
A CT number of a voxel in X-ray CT images represents the linear attenuation coefficient (LAC) of materials constituting the voxel. It is known that when a voxel is composed of multiple end members, the LAC of the voxel ( mix μ ) can be described as the combination of the LAC of each end member ( i μ ), which is weighted by the volume fraction of the end member. This is called the partial volume effect (Johns et al., 1993).
In other words, a CT number of a voxel, which is composed of unresolved pore and aggregate of solid particles, resides within the range between both LACs of air and solid particles. Thus, a CT number is proportional to the voxel porosity owing to the volume fraction of the unresolved pore. Accordingly, the partial volume effect majorly affects images as the resolution of images becomes low, because the probability that a voxel is composed of single material decreases. Consequently, the histogram of X-ray CT images with low resolution becomes unimodal shape as the amount of the voxels in the middle range of CT number histogram increases (i.e., the amount of the voxels affected by the partial volume effect increases). Compared to the image with low resolution, the typical shape of histogram for the high-resolution X-ray CT images is bimodal (Fig. 1). The binary segmentation schemes can make the accurate estimation of transport properties difficult as they classify the mixture voxel either into pure solid (so that the flow path of the voxel is neglected), or pure pore (so that the flow path is overestimated). In this study, by using pore size distribution (PSD) obtained from mercury intrusion porosimetry test, we presented a flexible ternary segmentation classifying X-ray CT images into apparent pore, gray pore, and apparent solid: apparent pore stands for the voxel composed of pure void, gray pore stands for the mixture voxel composed of unresolved pore and solid particles, and apparent solid means the voxel composed of almost pure solid. Secondly, we assigned a hydraulic diameter on the gray pore by correlating the volume information of PSD with the volume information of X-ray CT images. We implemented the proposed methods with Brinkman-force LBM (BF-LBM) and validated the proposed methods by confirming the accurate estimation of permeability of two rock specimens.

Rock specimens
Two types of rock specimens were prepared for this study; Navajo sandstone and Indiana limestone. X-ray diffraction test has shown that Navajo sandstone is mainly composed quartz, K-feldspar, and dolomite. Clay minerals such as montmorillonite, and illite are barely found to be neglected. For Indiana limestone, it is also confirmed that the specimen is composed of dolomite and the clay minerals are rarely found. Using mercury intrusion porosimetry test, the PSDs and porosity of specimens are obtained. The porosity of Navajo sandstone is 13.10 % and the porosity of Indiana limestone is 14.39 %. The permeability of Navajo sandstone is found to be in the range of 0.1 -100 mD (Dalrymple, 2007). The permeability of Indiana limestone is reported to be 9 -16 mD by the manufacturer (Kocurek Inderstries INC.)

Flexible ternary segmentation
After the acquisition of X-ray CT images, the images are segmented as three phases: apparent pore, gray pore, and apparent solid. Apparent pore is regarded as voxel porosity ( v φ ) of 1, and apparent solid has voxel porosity of 0. Gray pore has voxel porosity according to its CT number because of the partial volume effect, and the range of voxel porosity for the gray pore is between 0 -1 ( Fig. 2). To decide threshold for each phase, we correlate the histogram of X-ray CT and PSD. The apparent pore voxel has pore volume of voxel size 3 , because the voxel porosity of apparent pore phase is 1. Similarly, the gray pore voxel has pore volume of the combination of its voxel porosity and voxel size 3 . Based on this assumption, the total amount of pore volume obtained from CT number histogram should be equal to the pore volume obtained from mercury intrusion porosimetry test. Additionally, separating diameter (ds) is adopted. ds stands for the boundary of apparent pore and gray pore in the view of PSD. If ds is determined, pore with a diameter larger than ds is regarded as apparent pore volume. By matching the apparent pore volume obtained from CT number histogram with the apparent pore volume from PSD, it is possible to determine the threshold classifying apparent pore phase and gray pore  phase. The rest of pore volume is equal to the gray pore volume. By matching the rest pore volume from both PSD and CT number histogram again, the threshold distinguishing the gray pore phase from the apparent solid phase is determined.

Determination of hydraulic diameter gp d
To verify the sensitivity of the proposed method, the results corresponding to ds (= 1, 2 dx) are presented. Also, we attempted to assign hydraulic diameter on the gray pore voxel to consider the contribution of the unresolved pore. To determine the hydraulic diameter of the gray pore, it is assumed that the size of hydraulic diameter has linear relationship with CT number. When the CT number of the gray pore voxel is larger, the volume fraction of the unresolved aggregate of solid particles in the gray pore voxel is larger so that the hydraulic diameter of gray pore is smaller. Thus, the gray pore voxel with the largest CT number has the smallest hydraulic diameter. From this gray pore voxel, the corresponding pore volume is calculated, and the hydraulic diameter is determined by matching the pore volume of CT number histogram to that of PSD.
Finally, the gray pore volume of CT number histogram is correlated with the pore volume of PSD, and consequently, each gray pore voxel has hydraulic diameter from PSD according to its CT number.

Changes in transport properties according to separating diameter
Using the proposed methods to define the threshold for the three phases and define the hydraulic diameter for each gray pore voxel, voxel porosity, and voxel permeability ( v k ) for 3D voxel map of X-ray CT images are determined (Fig. 3). Two columns of X-ray CT images for each specimen show the 2D slice images at the different depth of the subdomain.
ds is 1 dx (4.17 μm) and 2 dx (8.34 μm) for Navajo sandstone. X-ray CT images shows that part of pore structure is filled with materials having very bright grayscale. Through SEM with EDS inspection, it is confirmed that this part is composed of dolomite minerals as the same as observed in X-ray diffraction test. The dolomite minerals have higher CT number than quartz, which is major component of sandstone. Voxel porosity map of Navajo sandstone shows that the solid phase composed of dolomite, quartz, etc is well classified as apparent solid phase having voxel porosity of 0.
Compared to voxel porosity map of Indiana limestone, voxel porosity map of Navajo sandstone shows clear separation of phases. The most of gray pore phases, having voxel porosity between 0 -1, is observed at the boundary of pore structure, known as most affected by the partial volume effect. From PSD obtained by mercury intrusion porosimetry test, the mean pore diameter of Navajo sandstone is estimated by 7.67 μm. Thus, the resolution of X-ray CT images is enough to resolve the major body of pore structure.
When ds increases from 1 dx to 2 dx, slight increases in the thickness of gray pore phases are observed, but the distribution of voxel porosity map doesn't change dramatically. Subdomains for 1 dx and 2 dx has the same bulk porosity because the bulk porosity is consistent regardless of separating diameter. The difference is the ratio of apparent pore phase and gray pore phase. Subdomains for this analysis have bulk porosity of 12.58%. We applied the proposed method to the cylinder-shaped bulk domain including most of X-ray CT images, and the cube-shaped subdomain is utilized as the domain for the flow simulation (500 3 voxels). Thus, the bulk porosity of the subdomain is not the exactly same as the measured porosity. We confirmed that the percentage of apparent pore volume accounting for the total pore volume over the entire domain is 60.25 % for 1 dx, and 47.46 % for 2 dx. This observation is consistent over the voxel permeability map. Changes in voxel permeability due to separating diameter doesn't seem noticeable. Voxel permeability in this study is determined by modified Kozeny-Carman equation for two parallel plate model (Eq. 2).
The voxel permeability is calculated by using the voxel porosity and hydraulic diameter introduced in the previous section.
ds for Indiana limestone is 1 dx (11.87 μm) and 2 dx (23.74 μm). The bulk porosity of subdomain used in this study is 13.53 %. The mean pore diameter is estimated by 2.10 μm. Although the X-ray CT images of Indiana limestone show distinguishable pores having a size of approximately 0.2 -0.3 mm, a considerable amount of pore is unresolved in the images. Thus, the voxel porosity map of Indiana limestone shows the unclear separation of phases compared to Navajo sandstone. In other words, gray pore region is observed in the solid region, as well as boundaries of distinguishable pore in the X-ray CT images.
It is noticeable that the area of gray pore phases in the voxel porosity map increases due to changes in ds, contrary to the subdomain of Navajo sandstone, The area of apparent pore phases decreases. It is confirmed that the percentage of apparent pore accounting for the total pore volume is 31.33 % for 1 dx and it decreases as 17.00 % for 2 dx.
The changes seem drastic through the observation on the voxel permeability map. The yellow area corresponding to apparent pore phase decreases. But the light blue area of gray pore phase increases so that the overall permeable area increases.

Changes in connectivity of pores according to separating diameter
We investigated the composition of connected pore affected by flexible ternary segmentation. Pore is classified as 5 groups: connected apparent pore (CAP), semi-connected apparent pore (SAP), disconnected apparent pore (DAP), connected gray pore (CGP), and disconnected gray pore (DGP). CAP means that apparent pore is connected through the flow direction from the inlet boundary to the outlet boundary. SAP means that apparent pore is not connected by itself, but the apparent pore is connected by the interplay with gray pore. CGP means that gray pore is connected either by itself or by the interplay with apparent pore (Fig. 4).
Navajo sandstone (NS) has CAP, while Indiana limestone (IL) doesn't have CAP. This can be caused by the difference in the resolution of images. However, considering the mean pore diameter of IL (2.10 μm), it suggests that most of pore structure of IL is composed of large pores connected by the small pores which are unresolved in the images. This shows that flexible ternary segmentation ensures the conservation of the connectivity of pores well, regardless of the rock types. As ds increases, the portion of SAP for IL decreases. Changes in ds decrease the apparent pore volume so that the amount of decreased apparent pore volume is converted into the gray pore volume For NS, small portion of DAP and DGP is observed. DAP and DGP mean isolated apparent pore and isolated gray pore, respectively. As ds increases, the portion of CAP decreases but that of CGP increases, so that the connectivity of overall pore structure increases. Thus, the increases of separating diameter cause the decreases of DAP and DGP.

Estimation permeability by using BF-LBM and its sensitivity to separating diameter
To verify the feasibility of the proposed method, we performed the flow simulation using Brinkman-force lattice Boltzmann model (Kang et al., 2018). Brinkmanforce lattice Boltzmann model (BF-LBM) solves different governing equations in both apparent pore region and gray pore region (Fig. 2). In apparent pore region, which has voxel porosity of 1, the conventional Stokes equation is solved. (4) where p is fluid pressure, μ is dynamic viscosity, and u is fluid velocity. In gray pore region, which has voxel porosity between 0-1 and voxel permeability calculated from voxel porosity and hydraulic diameter, Brinkman-Stokes equation is solved.
where μe is effective dynamic viscosity, and kv is voxel permeability.
The estimated permeability of NS is 70.72 mD for ds of 1 dx and 53.38 mD for ds of 2 dx. The estimated permeability of IL is 6.44 mD for ds of 1 dx and 8.98 mD for ds of 2 dx. It is confirmed that the results based on the transport properties estimated by the proposed methods are in a reasonable agreement with the expected range of the permeability (0 -100 mD for NS, 9-16 mD for IL). The decrease of CAP for NS due to the increase of ds leads to a decrease of permeability, despite the increase of CGP affected by the drag force corresponding to the last term of Eq. 5. For IL, the increase of ds induces the increase of permeability. This suggests that, in the case of the pore structure with no CAP, the increase of ds majorly contributes to the increase of the connectivity of pore due to CGP. However, despite the variation of estimation due to ds, the proposed methods and framework over this study still present a reasonable estimation, because the variation of apparent pore volume is automatically compensated by the changes in gray pore volume and the improvement of connectivity due to the gray pore.

Conclusions
In this study, we propose a method to assign the transport properties on the voxels of X-ray CT images of rocks by using PSD. For the implementation, we present flexible ternary segmentation with separating diameter (ds), and the procedure to determine the hydraulic diameter of the gray voxel. We confirmed that the proposed method gives reliable estimation of permeability, although the composition of three phases varies with the separating diameter. This suggests that, by employing the proposed method, one can perform analysis of materials characterization flexibly without the concerns about the accurate determination of threshold for the image segmentation. High-resolution 3D X-ray CT has been widely utilized over the industrial and academic fields, but the acquisition of X-ray CT images with good quality is still challenging because the settings for the imaging should be varied with the materials density and its composition. This study presents the method to overcome the limited size of resolution by gray pore phase having hydraulic diameter to conserve the connectivity of unresolved pore.