An intercomparison of the pore network to the Navier–Stokes modeling approach applied for saturated conductivity estimation from X-ray CT images

Different modeling techniques can be used to estimate the saturated conductivity of a porous medium based on computed tomography (CT) images. In this research, two methods are intercompared: direct modeling using the Navier–Stokes (NS) approach and simplified geometry pore network (PN) modeling. Both modeling approaches rely on pore media geometry which was determined using an X-ray CT scans with voxel size 2 μm. An estimate of the saturated conductivity using both methods was calculated for 20 samples prepared from sand with diverse particle size distributions. PN-estimated saturated conductivity was found to be statistically equivalent to the NS-determined saturated conductivity values. The average value of the ratio of the PN-determined conductivity to the NS-determined conductivity (KsatPN/NS) was equal to 0.927. In addition to the NS and PN modeling approaches, a simple Kozeny-Carman (KC) equation-based estimate was made. The comparison showed that the KC estimate overestimated saturated conductivity by more than double (2.624) the NS estimate. A relationship was observed between the porous media specific surface and the KsatPN/NS ratio. The tortuosity of analyzed samples was estimated, the correlation between the porous media tortuosity and the specific surface of the samples was observed. In case of NS modelling approach the difference between pore media total porosity and total porosity of meshes, which were lower, generated for simulations were observed. The average value of the differences between them was 0.01. The method of NS saturated conductivity error estimation related to pore media porosity underestimation by numerical meshes was proposed. The error was on the average 10% for analyzed samples. The minimum value of the error was 4.6% and maximum 19%.


Results and discussion
Pore space porosity analysis. The saturated conductivity for all examined samples was estimated using three different modeling approaches. Thresholded images of analyzed pore space gathered from CT scans were the initial source of information for all of the intercompared modeling approaches. However, the information about pore space geometry was used differently by each of the modeling techniques. The NS modeling methodology considers the complete physical description of the transport phenomena and relies on true pore space geometry in the numerical mesh generation procedure. Meanwhile, pore network modeling relies on a simplified PN consisting of spherical pores and cylindrical throats.
In the KC approach, pore space geometry was considered only in the form of two global measures: total porosity and specific surface (presented in Table 1). Specific surface and total porosity values for each of the pairs of samples, prepared from the same sand (Table 2), were similar, which indicates that the specimens used for CT scanning were successfully homogeneously prepared. The total porosity of the thresholded 3D CT image (Φ IMG ) used for the NS mesh and PN network generation was treated as a reference for comparisons with the total porosities of the NS meshes, Φ NS , and PN networks, Φ PN , which should be the same- Fig. 1 presents this dependence. Pore network porosity, Φ PN , has approximately the same values as image porosity, Φ IMG . As can be seen in Fig. 1, Φ PN followed almost a 1:1 relationship with Φ IMG . We may conclude that the PN network generation procedure worked well in this context.
In contrast, a comparison between Φ NS and Φ IMG reveals systematic underestimation. The average value of the differences between Φ IMG and Φ NS was 0.01 [m 3 m −3 ]. The total porosity of NS meshes was underestimated by 1.5% to 5.8%. Such differences in the total porosities of the numerical meshes used for these simulations and total porosities of pore space images cannot be simply neglected, as they could potentially impact the calculation results.
To deal with this issue, the error of saturated conductivity determination using the NS approach was estimated using proposed error estimation approach (Eq. 11). This approach assumes that the NS modeling error is related to the decrease in total porosity of a pore space following the same dependence as in the KC modeling approach. The absolute and relative values of these errors are presented together with estimated K sat values in Table 1. The error due to numerical mesh total porosity underestimation introduced via K sat estimation using the NS approach was in the range of 4.6-19% of the estimated K sat value. The error was on the average 10% for analyzed samples. The relative error seem to be higher for fine grained samples. The lowest values of the relative error was achieved for samples prepared from coarsest sand material.
Technical aspects of simulation and modeling. Table 1 also provides information about the generated numerical meshes and pore networks. The generation of the NS meshes and PN networks was the most computationally demanding task in this study. A multicore workstation with 512 GB RAM was used for this purpose. The size of the RAM was a limiting parameter that prevented the generation of more detailed NS meshes. The NS meshes that were generated consisted of 20 million cells in the case of the simplest pore space geometries and up to almost 90 million cells in the case of the samples with the highest specific surface. The files containing NS meshes have sizes ranging from ~ 3.5 to 11.5 GB. NS mesh generation using parallel mesh generation software took anywhere from five to 55 h. Pore network generation took a similarly long time, from 10 to 30 h each. Pore network generation time was substantially shortened from initial runs by parallelization of the initial ball searching step of the original 17 pore network generation code. Compared to the NS mesh and PN network generation Table 1. Pore space geometry, numerical mesh, and pore network stats of samples. During the automated procedure of mesh creation based on pore space geometry, despite quality constraints enforced by the mesh creation utility, some of the cells were of bad quality. These cells frequently displayed a large skewness that could impair the stability and quality of the calculations. As a final step in the mesh creation procedure, these bad cells were simply removed from the mesh. The number of cells removed from the NS meshes due to their failure to meet quality constraints is also presented in Table 1. Generally, the number of cells removed was minuscule and could be safely neglected. It is worth mentioning that removing low-quality cells had nothing to do with the NS mesh total porosity underestimation. On average, 0.015% of cells were removed from each NS mesh, while the level of total average porosity underestimation was 3.6%.

Name of sample
The total porosity underestimation for NS meshes could be attributed to several different sources. First, due to the constraints imposed on the size of the smallest possible cell size, the mesh generation utility could not generate a mesh in narrow regions where sand grains touch or are close to each other. Fortunately, these narrow regions should not greatly impact fluid flow, as in the case of a full saturation flow, it is a known phenomenon that the main flowpaths to transport dominate contributions to the overall pore space flow. Second, a possible cause of total porosity underestimation is the impossibility of mesh generation in the detached areas of a pore space. However, in this case, this should not be a problem, as these regions naturally do not take fluid flow into account.
When the numerical mesh was generated, NS equations together with boundary conditions allowed for the determination of the pressure and velocity fields in the simulated porous media flow. Figure 2 presents an example of the longitudinal cross-section along the direction of water flow in sample 1d. The uniformly changing pressure can be observed in Fig. 2a and has a pressure gradient enforced by the NS model boundary conditions. The velocity magnitude field is presented in Fig. 2b, in which areas of higher velocity, caused by the narrow pores, may be observed.
Saturated conductivity-the dependence on the modeling approach used for estimation. The NS modeling methodology, unlike PN modeling, accounts for full physical descriptions of the transport phenomena. Additionally, NS modeling relies on true pore space geometry and is presented here as a reference.
Other modeling approaches (e.g., PN and KC) were compared with the NS model. The results of this comparison are shown in Fig. 3 together with linear models linking saturated conductivities determined by the PN and KC methodologies with their NS-based counterparts. Based on these results one can conclude that, compared to the NS model, the KC-based modeling approach systematically overestimates saturated conductivity values, as the slope parameter of the linear model (LM) is greater than one (2.169) and the intercept is positive. The average ratio of K satKC to K satNS (K satKC/NS ) is even larger than the slope parameter and is equal to 2.624. On average, samples analyzed using the KC equation-based estimations are almost three times greater than NS-based K sat estimations. The difference between K satNS and K satKC can also be observed in the subplots of the boxplots presented in Fig. 3.
The KC equation is known to be only a rough estimation 40 of conductivity, though some discrepancy is not unexpected. Different studies have shown that KC underestimates 40,41 or overestimates 42 saturated conductivity depending on the porous medium for which the KC equation is used. In the case of Kuang et al. 40 , the authors www.nature.com/scientificreports/ used a tortuosity estimation method based on CT porous media images, which generated high tortuosity values of ~ 4.5 that could lead to conductivity underestimation. Our results, as they relate to the KC saturated conductivity estimation, are not as overestimated as those reported by Mostaghimi et al. 42 , in which KC estimation was also compared with NS estimation. This difference might be due to the different pore space characteristics of the porous material studied. The greatest overestimation was reported by Mostaghimi et al. 42 for samples with porosities lower than 0.25, which were not observed in the samples used in this study. Taheri et al. 43 reported high levels of agreement between KC and NS saturated conductivity estimations, but their work was related to a virtual homogeneous sand porous medium built from spherical particles. This medium very closely resembled the KC equation assumptions. The accuracy of KC estimations may be problematic, but considering how easily this estimation can be performed in comparison to the application of NS and PN models, the KC method is suitable in some situations. Estimations of saturated conductivity via PN modeling are also presented in Fig. 3. In the case of PN estimations, the LM model slope parameter is equal to 1.593, but the intercept is negative. When the average ratio of the K satPN to K satNS (K satPN/NS ) is calculated, it equals 0.927. This result shows that PN-based estimation of the saturated conductivity gives, on average, practically the same results as NS-based estimations. Boxplots (Fig. 3) confirm that these two sets (i.e., K satPN and K satNS ) are comparable. This conclusion is also confirmed by a two-sided paired     26 analyzed one artificial sphere pack porous medium and found high agreement between PN-and NS-estimated saturated conductivity. Finally, real 3D porous media were also analyzed by Jiang et al. 3 which demonstrated similarly good agreement when three low resolution (256 3 voxels @ 6 µm) rock samples were studied. Our study, which includes larger and higher-resolution real porous media geometry (~ 1400 3 voxels @ 2 µm), confirmed the equivalence of NS and PN estimations. The relationship of NS-determined saturated conductivity to its PN and KC counterparts, expressed by the average of the aforementioned ratios (i.e. K satKC/NS and K satPN/NS ), does not completely describe the situation. Regression analysis demonstrates the dependence of both properties, K satKC/NS (r = 0.72) and K satPN/NS (r = -0.86), on the sample's specific surface, determined based on thresholded images-the same images which were used for the PN network and NS mesh generation. This dependence is presented in Fig. 4, where K satKC/NS falls within the range of 2.2-3.5 and K satPN/NS is in the range of 0.4-1.5.
It is important to note that if K satKC/NS increases, but K satPN/NS decreases with the increase of the specific surface σ [m -1 ]. This means that these opposite trends in the ratios are not caused by the dependence between K satNS and σ, as if that were the case, the trends would follow a similar pattern. To ensure the possible influence of the K satNS estimation error on the values of K satKC/NS and K satPN/NS , their errors (presented as error-bars in Fig. 4) were calculated based on error propagation theory. Since the increasing specific surface is related to the particle size distribution of the material from which the samples were prepared, which in turn influences pore space complexity, we hypothesized that the observed trend of K satKC/NS might be connected to the changing tortuosity of the pore space.
To validate this hypothesis, the tortuosities of the samples were estimated (see Table 1) based on the streamlines analysis. The lowest values of tortuosity for the bed of spherical grains, close to the theoretical value of 1.414, were estimated for samples 1u and 1d, which were made of the coarser fraction of sieved sand grains. However, the majority of samples were made of the milled sand, which consisted mostly of non-spherical grains where one might expect higher tortuosity to occur. As a result, tortuosity is positively correlated (r = 0.71) with the image specific surface and the highest values of tortuosity ~ 1.6-1.7 were achieved in samples with a higher fraction of finely-milled sand.
When these tortuosity estimations were used in the KC Eq. (4) instead of the constant value 1.414, the KC estimations no longer depended on the value of σ (see data series K satKC/NS (tc) in Fig. 4). This confirms the hypothesis about tortuosity dependence on the pore space complexity.
The results obtained from numerical estimations and from Kozeny-Carman approach were also subjected to experimental verification (Fig. 5). As per one specimen, for which K sat was measured, two samples belongs simulated and measured values could not be compared directly. Saturated conductivities estimated for the samples were recalculated to the specimen's values according to Eq. (2). In case of NS estimations, error using proposed www.nature.com/scientificreports/ approach (Eq. 11) was also for all samples determined. On the figure (Fig. 5) for all specimens K satNS errors are presented-they were recalculated based on error propagation theory from appropriate samples' error estimates.

Conclusions
The PN-estimated saturated conductivity (K satPN ) was found to be statistically equivalent to NS-determined saturated conductivity values (K satNS ). The average value of the K satPN /K satNS ratio was equal to 0.927. The KC equation-based estimation (K satKC ) overestimated saturated conductivity by more than double (2.624) the K satNS estimate. However, despite high levels of overestimation, both sets of values are well-correlated with R 2 = 0.99.
The total porosities of the generated PN were the same as the total porosities of images used for their creation. Total porosities of numerical meshes were on average lower by 0.01 [m 3 m −3 ] than the porosities of the corresponding reference images.
Proposed error estimation approach was used for relative error determination of the NS saturated conductivity modelling. The error was on the average 10% for analyzed samples. The minimum value of the error was 4.6% and maximum 19%. The relative error seem to be higher for fine grained samples. The lowest values of the relative error was achieved for samples prepared from coarsest sand material.
The tortuosity of 20 samples was estimated. The value of the tortuosity determined for the samples prepared from coarser fractions of sand approximated the theoretical tortuosity value of the porous medium built from spherical grains, 1.4. The tortuosities of samples prepared from milled sand (non-spherical grain material) were higher. The value of tortuosity was correlated (r = 0.71) with the porous medium-specific surface, reaching the highest values of 1.6-1.7 for samples prepared from the finest sand material.

Materials and methods
Specimens preparation. Specimens were prepared by mixing different unimpaired sieved and milled sand fractions. The first step in specimen preparation was to sieve the sand using a 0.5 mm sieve. This raw sand was then milled in the planetary mill (Pulverisette 6 classic line, Fritsch, Germany, Idar-Oberstein). The raw sand was milled for 5, 10, and 20 min to achieve milled sand fractions with different particle size distributions.
10 specimens were then prepared from a mixture of sieved sand and milled sand fractions. Specimens s1, s2, and s3 were prepared by sieving raw sand through 0.32, 0.16, and 0.08 mm sieves, respectively. Table 2 presents the mixing scheme used to create material with diverse particle size distributions (PSD) for the other specimens. Each prepared specimen's material prior to the specimen preparation was measured using a laser diffractometry method (LDM) to verify the PSD that was achieved (Fig. 6). A Malvern Mastersizer 2000 with a measurement range of 0.02 µm to 2 mm was used for PSD determination. To obtain homogeneity in the sand suspension, a Hydro G dispersion unit was used, with a pump speed of 1750 rpm and stirrer speed of 700 rpm. Light intensity, measured with the detectors, was recalculated into PSD according to the Mie theory ( ISO 13,320:2009ISO 13,320: , 2009). The Mie model parameters include an absorption coefficient of 0.1 and a refractive coefficient of 1.52 45 . Due to this specimen preparation procedure, it was possible to evaluate the potential impact of the pore space geometry on saturated conductivity modeling.  www.nature.com/scientificreports/ The specimen was prepared in a polypropylene tube (4 mm internal diameter) of low X-ray absorbance (Fig. 7). An adequate amount of the specimen's material was poured into the cylinder to achieve a 10 mm stack inside the container. The specimen was subsequently compacted in the pipe by vibration. At the bottom and top of the specimen's material stack bronze meshes were placed to keep the material in place. In the case of the specimens (s4, s5, s6, s7, s8, s9, and s10) additionally, the circle cut from filter paper was placed between the mesh and the material to avoid loss of the fine particles from the specimen. Bronze meshes were used as they were stiff enough to support properly the specimen's material, but there was a drawback of this choice. As the metal is far less penetrated by X-rays than sand, strong streaking artifacts are observed in the reconstructed scans near the mesh. As a result, the pore space in the neighborhood of the bronze meshes can't be observed. That was also a reason for the preparation of relatively high (10 mm) specimens as we would like to ensure a big enough part of the specimen could be CT analyzed to ensure the reasonable comparison between modeling and measurement of the saturated conductivity. As the height of the area which is examined using CT is roughly equal to its diameter, in the 4 mm wide and 10 mm high specimen two samples area were defined.
X-ray CT sample analysis. An X-ray CT scanner (GE Nanotom 180S) was used to examine the specimens with a 180 kV/15 W microfocus X-ray tube. For each specimen, two scans were made, one in the upper part of the specimen and the second in the lower part (Fig. 7). Two regions from each specimen were scanned to achieve the highest possible CT resolution, which limited the dimensions of the scanned object. This also increased the number of analyzed samples. As a result, 20 scans were made. They will be referred to as "10d" and "10u, " for example, where "10d" identifies the lower (downward) part of the 10th specimen and "10u" is the upper part of the 10 th specimen.
For each scan, a series of 1200 2D radiograms were collected, with the specimen rotating through the full angular range using a rotation step of 0.3°. Each 2D radiogram was averaged over 15 images to reduce noise. The images were recorded using a detector with a resolution of 2284 × 2304 pixels registering images at the 14-bit gray-level depth. The X-ray source operated at 90 kV with a 120 μA cathode current and tungsten exit window. Immediately before each scan, a short pre-scan was performed, which lasted for 30 min to preheat the specimen and minimize the impact of thermal expansion during the primary scan for three-dimensional (3D) image reconstruction. The voxel size achieved in these CT scans was 2 µm. 3D image reconstruction and processing. Tomographic reconstruction was performed using a series of radiograms (DatosX version 2.0.1, General Electric). The beam hardening correction was not necessary, although the exit window filter was not used during the X-ray scan. From the obtained 3D volume representing the sample and a part of the tube, a cubic region of interest (ROI) of 1400 × 1400 × 1700 voxels was chosen for further processing. Cross-sections of the ROIs of selected samples are presented on the Fig. 8. A median filter with a three-pixel kernel diameter was applied to the ROI. The next step was thresholding using an IsoData thresholding method 46 . Image processing was conducted with Fiji and VGStudio MAX 2.1 (Volume Graphics). The thresholded images serve as a common starting point for both the PN and NS modeling approaches.

Measurement of the specimen's saturated conductivity.
For each of ten prepared specimens, saturated water conductivity was measured. The measurement was based on the static pressure head principle (Fig. 9). The pressure difference was fixed, and water flux was measured, then saturated conductivity K [m s −1 ] was calculated using Darcy's law: where q [m s −1 ] is fluid flux, Δp is the pressure difference [Pa] along the sample length Δl [m], g is the gravitational constant [m s −2 ], and ρ is fluid density [kg m −3 ].
The pressure difference was adjusted by changing the water level in the water tank. The value of the difference in the pressure head was read from the water level difference in the two burettes mounted before and after the specimen. Water flux was measured based on the scale readings on which the container collecting the water flowing through the specimen was placed.
The idea of the measurement is simple, but due to very low water fluxes-especially for fine-grained specimens-some additional factors had to been taken into account. To avoid the capillary forces' influence on the scale reading, the water container located on the scale was initially filled up to the level little higher than the end of the glass pipe placed in it. Some of the measurements were long enough to force the water weight loss correction due to water evaporation from the container placed on the scale. The evaporation rate from the container in the constant air temperature kept in the laboratory was determined upfront. The last factor that needed correction was the influence of the bronze meshes and paper filter placed at the ends of the specimen's material. Five empty (consisting only from meshes and filter paper) specimens were prepared, and dependence between water flux and pressure head difference was determined for them. Based on averaged results calibrating function linking flux with pressure head difference was prepared. Then during regular measurements, based on registered water flux, appropriate pressure head difference was subtracted from the value read from the burettes.
The results of measured specimen's saturated conductivities were compared with values of samples conductivities estimated using different methods. There values cannot be compared directly. But taking into account that two subsequent samples are defined in one specimen (Fig. 7) one may calculate effective water conductivity K eff of two stacked samples with conductivities K 1 and K 2 , as a harmonic average: www.nature.com/scientificreports/ Pore network saturated conductivity modeling. The thresholded 3D images were processed using the network extraction code after applying a maximal ball algorithm 17 . The result was a dataset containing information about the spherical pores and cylindrical throats including dimensions, spatial coordinates, and connections 47 . A simplified pore network (Fig. 10c) was then used for permeability calculations based on the assumption of laminar flow, the Hagen-Poiseuille formula validity for flow in the throats, and mass conservation for each pore-throat connection 48 . The boundary condition pressure difference was set between the PN input and output nodes, resulting in fluid flow throughout the network.
Navier-Stokes saturated conductivity modeling. The NS FVM approach is used to model a singlephase, laminar, steady-state flow of incompressible fluid through a porous medium. The set of equations used for the model contains the following momentum balance equation: (2) To solve the NS equations, simpleFoam software was used, which is a part of the OpenFOAM CFD toolbox. The solver implements a semi-implicit method for pressure-linked equations (SIMPLE) algorithm to solve the momentum balance equation.
Appropriate boundary conditions were established in the simulation. The fluid velocity was fixed to 10 -5 m s −1 , and the pressure was fixed to 0 Pa on the input patch and zero gradient Neuman conditions on the output patch. Input and output patches were equivalent to the top and bottom sides of the sample. No-slip boundary conditions were applied to the pore walls. The saturated conductivity was calculated from Darcy's law using the pressure drop and fluid velocity at the input boundary.
The first step of an NS simulation using FVM is mesh preparation. The surface of the solid fraction was determined from the thresholded 3D images and approximated using a triangulated surface mesh (Fig. 10a), which was saved in an STL format. For pore space surface triangulation, Fiji software was used along with an isosurface module from the BoneJ plugin 49 . The mesh (Fig. 10b) was then created based on the STL file using the snappyHexMesh utility (a part of the OpenFOAM toolkit). Next, the generated mesh was checked for highly skewed cells and other deficiencies. Cells that did not meet quality constraints, in particular those with skewness greater than four, were removed to avoid impairing the quality or numerical stability of the simulations. A very small number, not greater than 0.0282% of the total number of cells, had to be removed.   . Different values of C 0 and τ may be used for conductivity estimation. In this study, we used the most commonly used values, which are valid for equally-sized spherical granular material: C 0 = 2.0 and τ = 1.414. Total porosity and specific area were determined based on the thresholded 3D images of the samples using an image analysis approach. Total porosity was determined using the voxel counting method. The specific surface was determined based on the surface area of the triangulated solid-void interface approximation.
Saturated conductivity estimation errors. Saturated conductivity estimation relies on the modeling of numerical meshes. In principle, numerical meshes that are generated based on images of a porous medium should match its overall geometric characteristics (e.g., total porosity). Any discrepancies may lead to errors in the saturated conductivity estimation 50 . To account for this, some method of saturated conductivity error estimation is needed. Ideally, error estimation should be based solely on numerical modeling methodology, but it is not possible. As a an error estimation method, following KC equation-based approach can be used, as it relates total porosity and specific surface to saturated conductivity. Proposed approach assumes correlation between K NS and K KC . This correlation in case of samples analyzed here was 0.99. And of course validity of Kozeny-Carman estimation for pore media in question. Correlation between K NS and K KC means that the following dependence is assumed: where A is some constant and explicitly dependence on mesh porosity φ was noted. Based on error estimation theory, the error ΔK NS related to porosity underestimation �φ will be calculated using following equation: But taking into account Eq. (6) following equation will be derived: Based on Eq. (5), using the dependence of KC saturated conductivity estimation on total porosity, Eq. (8) become: Which, using Eq. (5) again could be rewritten as: which finally, after taking into account Eq. (6), becomes the formula for relative K NS error estimation: The relative and absolute errors estimated based on this formula are presented in Table 1.

Tortuosity estimation.
Porous medium tortuosity is one parameter involved in saturated conductivity estimation using the KC approach. Tortuosity is related to pore space geometry and pore space complexity. It is also defined as the ratio of the length of a path between two points to the distance between them.
The tortuosity of porous media can be estimated using either direct analysis of possible paths in 3D porous media images 51 or analysis of a simulated velocity field in a pore space. The velocity field analysis approach may be used to deduce tortuosity based on the averaged velocity magnitude to longitudinal velocity component ratio 52 or streamlines length and distance analysis 53 . The latter approach was used in this work.
For all samples, fluid velocity fields were simulated using an NS modeling approach. Based on the information about pore geometry and the velocity field, streamlines were generated. For each sample, 1000 randomly Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.