Correct estimation of permeability using experiment and simulation

Permeability (or the related permeance, hydraulic conductivity, conductance, drag, or friction coefficient) quantifies the ability of a porous geometry to conduct fluid. Permeability estimations are relevant for the understanding of soil infiltration,^1,2 sea ice evolution,³ durability of concrete,⁴ carbon dioxide sequestration,⁵ hydrocarbon recovery,^6–8 drug delivery,⁹ groundwater flow,¹⁰ and fluid mechanics of tumors.¹¹ Knowledge of permeability is essential in the design of fuel cells,¹² lithium-ion batteries,¹³ paper-based microfluidic biosensors,¹⁴ face masks,^15,16 hydrogels,¹⁷ engineered bones,¹⁸ and textiles.¹⁹ 

Hydraulic conductance K was introduced by Darcy²⁰ in an empirical relation between volumetric flow rate Q, imposed hydraulic head loss $Δh$ over length L of a vertical tube of cross-sectional area A: $Q=KAΔh/L$⁠. This equation was generalized by Hubbert²¹ to a physical law. For the horizonal, isothermal, steady-state discharge of incompressible fluid into the atmospheric pressure P₀, the intrinsic permeability $k*$ is

where μ is the dynamic fluid viscosity and $Pin$ is the inlet pressure. Later, the theoretical links between the Darcy and differential Stokes equations were provided by, e.g., Whitaker.²² Here, we employ a dimensionless form of intrinsic permeability, k, which is $k*$ divided by the square of a characteristic length. Applying Darcy's relation (1) is the simplest, standard way for estimating permeability of a porous sample.

Despite apparent simplicity of permeability measurements, the reported values are so scattered that the relative difference within a decade (i.e., 1000%) for similar samples is not uncommon. At the same time, estimations of pressure, flow rate, sample dimensions, and viscosity all have 1% accuracies or better. In some cases, the differences in measured permeabilities are genuine and caused by natural variation of similar samples. Even in this case, the problem we face is when permeability is measured, it is difficult to assess its value accurately, and tenfold differences, real or imaginary, seem to be acceptable in some fields.^3,23–26

Computer simulations can serve as an independent tool for estimating permeability of real samples. The simulations of flow are performed directly on a three-dimensional (3D) image of sample pore space. After scanning the interior structure of a 3D sample with, e.g., x-ray computed tomography (CT), magnetic resonance imaging, or confocal laser scanning microscopy, the output 3D image is processed and segmented, resulting in a set of binary (solid or fluid) voxels, which are used as an input to flow simulations. Here, “simulation of flow” is defined as finding a fluid velocity vector for each fluid voxel, and up to billions of such vectors that fully describe the 3D flow field within the sample [e.g., Fig. 1(e)]. The 3D computer models are capable of providing sub-percent accuracy. However, several-fold differences between simulations and experiments are not uncommon.^23,27–33

Comparison of the simulated permeability with experimental data includes the following four steps: (i) measuring the sample permeability, (ii) scanning the sample, (iii) processing and segmenting the resulting 3D image, and (iv) performing flow simulations on the segmented 3D image. This procedure is highly sensitive to all of the aforementioned steps. The picture is further complicated by the different signs of errors originating from each step, allowing them to cancel one another. Non-optimal design, improper execution, or incorrect interpretation of any single step ruins the final agreement and precipitates attempts to find the origin of the disagreement or acceptance of the unexplained difference.

In this work, we follow the aforementioned four steps and obtain an accurate match between experimentally measured permeability and its simulated counterpart. Simulations are based on the solution of Stokes equation provided by the lattice Boltzmann method and are free from discretization resolution impact, computed-tomography (CT) artifacts, and operator-dependent segmentation. Matching experiment with simulation reveals a flaw originating from the out-of-sample pressure port placements and improper pressure gradient calculation that also affected our own measurements. Detailed investigation of these disagreements suggests a similar flaw in the current industry standards, such as ASTM D5084 (D5856).

This study is organized as follows. We briefly describe the employed experimental, imaging, image processing, and simulation approaches. We then share our initial experience, where the unexplained difference (flaw) is observed between experiment and simulations using a representative volume of the sample. We compare our initial results with the existing literature data. Next, we provide an updated simulation approach resulting in an accurate match between experiment and simulations. Hereafter, we introduce a simulation setup to mimic our experimental setup and to reproduce the observed flaw. We use this simulation setup to study the influence of its geometrical parameters on the flaw magnitude. Finally, we provide a simple empirical equation for the flaw magnitude and discuss its minimization. The main text is followed by the  Appendixes.

To justify our findings, a large number of technical details needs to be explicitly stated and quantified due to the multidisciplinary nature of this work. Therefore, we support this study with extensive  Appendixes describing our experimental, imaging, and simulation procedures. These  Appendixes are summarized in this paper.

We use dense irregular packings of glass beads as porous samples for permeability measurements. Packed glass beads (1) create a complex network of pores between the impermeable glass phase, (2) allow control of the discretization resolution and Reynolds number by choosing an appropriate bead diameter, and (3) enable comparison with the previous studies. We employ two types of commercially available ∼0.5 mm glass beads referred to as “beads1” and “beads2.” We use the Sauter mean diameter (⁠$d32=∑idi3/∑idi2$⁠) as a characteristic length for each bead type for permeability nondimensionalization,³⁴ $d32beads1=470.3$ and $d32beads2=541.1$ micron. We note that a particular choice of the characteristic length establishes comparison of our permeability values with other studies but has no impact on the experiment–simulation comparisons in this work. Glass beads are packed in 9 mm glass tubes under ultrasonic vibration resulting in samples P3 (beads1) and P4 (beads2). The corresponding length and porosity of each sample are about 6 cm and 0.355, respectively. Please see Table I and Appendix A 3 for more details on the samples and packing procedure.

Glass beads are held inside glass tubes using custom-designed plastic plugs with 1.89 mm-diameter inner holes enabling in- and outflow through the bead packing. To prevent beads from entering the plug hole, we use custom-designed silicon meshes (Appendix A 2). Each sample is connected to a syringe pump with flexible tubing of 1.73 mm inner diameter [Figs. 1(a) and 1(c)]. The total length of plug holes and tubing is about 15.2 cm. Pressure values are recorded using a gauge-type pressure transducer connected to the tubing before the sample. The outlet is open to the atmosphere [Figs. 1(a) and 1(c)], and this results in the minor pressure oscillations related to droplet formation and detachment from the outlet [Fig. 1(b)]. Please see more detailed analysis of the open outlet impact on the accuracy of the recorded pressure drop in Appendix A 5.

We flow three glycerol–water solutions, which allow us to measure absolute permeability at different viscosities (21, 38, and 55 mPa s), confirming reproducibility of the experimental permeability values. Viscosity of each solution was determined using a rotational rheometer (MCR 302 by Anton Paar, Graz, Austria). More details of the viscosity measurements can be found in Appendix A 7.

Special attention is given to the stabilization of the experimental permeability values. As a result, we reduce the initial $∼15%$ scatter to $∼1%$ (Appendix A 6), and the measured permeability values have no outliers.

Using the x-ray CT scanner, the 3D images of sample pore space, including the confining wall, are acquired at 18 different discretization resolutions r for each of P3 and P4 samples. We define r as the number of voxels per Sauter bead diameter. All voxels outside of the confining wall are marked as wall. Hereafter, the acquired gray-scale images are segmented using a global threshold. The value of the global threshold is taken from the laboratory-measured porosity of each sample. Based on its gray value, each voxel is initialized as solid or fluid in order to match image porosity with its experimental value. Such segmentation procedure avoids operator-dependence of the output binary images (Appendix B 2).

Simulations of incompressible flow are performed from first principles and based on the solution of the continuity equation, $∇·v→(r→)=0$⁠, and the Stokes equation, $μ∇2v→(r→)=∇p(r→)−B→$⁠. The latter is a balance between viscous forces, pressure $p(r→)$⁠, and body force $B→$ acting on a small fluid element with dynamic viscosity μ and velocity $v→(r→)$ located at $r→$⁠. We solve for the velocity field $v→(r→)$ using the lattice-Boltzmann method.³⁵ We apply bounce-back boundary condition at solid–fluid and wall–fluid interfaces as well as periodic boundary conditions at the external boundaries of the computational domain.

Permeability values obtained from the simulations on CT images are impacted by the image resolution, r, and image contrast. Influence of the image contrast on the simulated permeability values decreases with the increase in r. To better understand and minimize the impact of the image resolution and contrast, we use computer-generated unconfined packings (with beads1 and beads2 diameters) in addition to the CT images and vary resolution for all these geometries. The computer-generated geometries are by definition free from the contrast impact. First, for each geometry, we determine $∼0.1$%-accurate reference permeability value (⁠$kP3,P4ref$ for CT images and $kbeads1,beads2ref$ for computer-generated geometries) using two approaches: extrapolation³⁶ and relative error correction.³⁷ Each reference value $kP3,P4ref$ is free from the resolution and contrast impact (Appendixes B 2 and C 3). We, thus, obtain the dependency of the relative error in permeability on the discretization resolution. The corresponding result for beads1 (computer-generated packing and P3 sample) is shown in Fig. 2(a), while the case of beads2 and more details can be found in  Appendix C.

For the fixed resolution of $r≈20$⁠, we performed several CT scans with different settings resulting in the 3D images of different contrast [red circle and all except red crosses on the vertical dashed line in Fig. 2(a)]. Robust segmentation procedure allows one to observe systematic decrease in permeability with the image contrast improvement (see also Appendix B 2). For some “optimal” CT settings we found [red circles in Fig. 2(a), Appendix B 2], the contrast impact on permeability vanishes at higher resolutions (r > 50), and the relative errors of the real and computer geometries overlap. After applying downsampling procedure to the high-resolution (r > 50) CT images, we obtained lower-resolution (r < 20) downsampled CT geometries of ultimate contrast, as demonstrated by their permeability error–resolution behavior overlapping with the corresponding computer-generated geometry [green circles and red cross in Fig. 2(a)]. Please note that P3 and the corresponding computer-generated geometry have a 20% difference in their reference permeabilities (⁠$5.82×10−4$ vs $6.85×10−4$⁠, Appendix C 3), but their relative errors overlap with the increasing resolution. For beads2, the behavior is identical, as shown in Appendix C 3. We are not aware of other studies showing such a clear match in the error–resolution behavior of real and computer-generated geometries.

According to Fig. 2(a), the relative error in permeability does not converge to zero with increasing resolution but crosses zero and remains negative due to slow convergence. This fundamental behavior of the discretization error can be observed elsewhere,^38–40 and it originates from the interplay between the convergence rates of the boundary condition (−1) and the numerical scheme itself (−2).³⁶ Therefore, the error in permeability is smaller at $r≈20$ compared with $r≈60$⁠, if we consider downsampled CT image of ultimate contrast [red cross in Fig. 2(a)]. We exploit this fact later in this study, where for the simulations on the sample+plugs+tubing system [Fig. 1(d)], we use the geometries of ultimate contrast at $r≈20$⁠.

Our first attempt to compare simulation with experiment is illustrated in Fig. 3(a). In Fig. 1(c), the permeability of the entire system between black dots with the P3 sample is measured, while simulations are performed on 80% of the sample length $Lsample$⁠. Here, we assume that solely the sample contributes to the pressure drop. Changing 80% to 60%, 70%, or 90% has unnoticeable impact on the simulated permeability, which confirms that 80% sample length, which is about 100 Sauter diameters, is representative. Experimental values are calculated according to (1) using $Lsample$ for the pressure gradient $(Pin−P0)/Lsample$—as far as we know the only accepted interpretation of the experimental setup in Fig. 1(c). This results in a ∼35% overestimation of the experiment by the simulations, Fig. 3(a).

The observed disagreement can originate from the pressure drop in the tubing (1.73 mm i.d.) and plugs (1.89 mm i.d.). To estimate it, we simulated flow in the empty system (i.e., without the sample installed), which is a series of coaxial hollow cylinders mimicking tubing, plugs, and 9 mm sample holder. Using fixed cross section area for this geometry (see Sec. III C), we estimate permeability for this empty system: the corresponding pressure loss is about 6% of the system with the sample installed. The similar result is obtained if we calculate the pressure drop in the empty tube of 1.8 mm i.d. and 15.2 cm length using the Hagen–Poiseuille equation, and compare it with the experimental pressure drop. The remaining 30% of the permeability difference between experiment and simulation is still unexplained, and it cannot be attributed to the precision of experimental or simulation procedures.

Let us assess our results by comparing them with literature. Figure 3 shows previous studies based on a) the flow simulations in computer-generated packings of equal spheres,⁴¹ experiments by (b) Fand et al.,⁴² and Macini et al.,⁴³ and (c) Wyllie and Gregory,⁴⁴ James et al.,⁴⁵ and Loudon.⁴⁶ Studies (c) placed the pressure ports inside of the sample, while studies (b) outside, connecting the pressure ports to the inflow and effluent tubes, see Fig. 3(b). According to our own experience and an interdisciplinary literature survey, external pressure ports prevail. Adding the theoretical Kozeny–Carman⁴⁷ estimation suggests that internal pressure ports, previous simulations, and the current simulations produce similar results, and that the current experimental value appears to be abnormal. Figure 3(a) also suggests that external pressure ports may or may not lead to a significant underestimation of permeability, e.g., Macini vs Fand's data (discussed below). All the mentioned experimental setups are depicted in  Appendix D.

We emphasize that placing pressure ports outside of the sample can be found already in the original Darcy's experiment [lower pressure port in Fig. 3(c)]. This raises the formal question of the impact of external pressure port location on the permeability measurements, which we address here.

The only option left to match the experiment and simulations is to simulate the full system between the black dots in Fig. 1(c), including the sample, plugs, meshes, and tubes. For this purpose, we attach hollow cylinders mimicking tubing and plugs to the CT sample image of $r≈20$ and ultimate contrast. We refer to the full system with the sample P3 as S3 and with the sample P4 as S4. The corresponding experimental permeability values are now obtained by not using $Lsample$ but $Lsystem$⁠: the pressure gradient is calculated as $(Pin−P0)/Lsystem$⁠, which changes experimental dimensionless permeability values from about $0.4×10−3$ to $1.5×10−3$ (compare Y-axes in Figs. 3 and 4). Hereafter, we obtained, in principle, an exact match between both approaches, Fig. 4. The viscosity and flow rate vary significantly (2.5 and 4 times, respectively) compared with the level of obtained mismatch (<1%). Each group consists of 16 experiments performed one-after-another and has no excluded outliers. Our goal to overlap experimental and simulation precisions is achieved. We emphasize that the absolute dimensional permeability values of both systems, $kS3=3.4×10−10$ $m2$ and $kS4=4.7×10−10$ $m2$⁠, differ by 38%. We are not aware of a similar level of agreement between the measured and simulated permeabilities when details sufficient to verify that the match is not a coincidence are also provided.

In order to match the experimental and simulated permeabilities, we need to calculate the seepage or superficial velocity from the simulated three-dimensional flow field, which is $εv→avg$⁠, where $v→avg$ is the velocity averaged over all fluid voxels and ε is the porosity. For this purpose, we need to know the porosity of the system. While for the sample region, its calculation is obvious (porosity is just the ratio of fluid to the total number of fluid and solid voxels inside the sample), the tube region requires more clarification. According to Darcy's law (1), the flow rate is calculated as Q/A using the single value of sample cross section A. In order to satisfy this definition, we need to “enclose” the tube into the same cross section area as the sample and consider the tube as a long single “pore” in such an enclosure. The porosity of the tube region is not 100% but $Atube/Asample$⁠. The tube with a 1.73 mm diameter enclosed in a 9 mm cylinder has the porosity of $(1.73/9)2≈3.7$%. Similarly, for the plug region, $(1.89/9)2≈4.4$%, Fig. 1(d). Using the relative length fractions of the tube, plug, and sample regions, the total system porosity becomes 13% (see also Table I). After connecting pressure ports via tubing with in- and effluent to the sample, one measures permeability of the virtual sample highlighted as the red solid but not as the white dashed contour in Fig. 1(c) or the black dashed contour in the following Fig. 5(a).

The match between experiment and the Stokes equation-based computer model with no fitting parameters (Fig. 4) implies that the Stokes physics is sufficient to mimic the key processes in the considered experimental system. From this point on, it is straightforward to use the full power of computer simulations to better understand the impact of tubes on the permeability measurements. A virtual sample can be created, and its correct permeability estimated and compared with the permeability of the geometry with tubes, also addressing the case of incorrect interpretation of the applied pressure gradient.

To address the impact of the out-of-sample pressure port connections on permeability measurement, we have designed the following simulation setup [Figs. 5(a) and 5(b)]: the sample is a cylinder of diameter $Dsmpl$ and length $Lsmpl$⁠, packed with spheres of diameter $dsp$ at the porosity of 40% (except for the “membrane” case below). During packing generation, the spheres are confined by the cylinder, and the hard walls at both cylinder ends. Each generated packing (sample) is (i) discretized at a resolution between 10 and 20 voxels per $dsp$ and (ii) mirrored in flow direction to obtain its reference permeability. Thereafter, a tube of length $Ltube$ is attached to the original (not-mirrored) sample as another hollow cylinder with $Dtube$ internal and $Dsmpl$ external diameters. The tube and sample define the system of length $Lsyst$ and permeability $ksyst$⁠. A special case is the membrane case, where $Lsmpl/Dsmpl=0.1, Dtube=Dsmpl$⁠, and $Ltube≫Lsmpl$ [Fig. 5(d)].

First, we determine the reference permeability for each sample $ksmpl$⁠, and then we find the permeability $ksyst$ for each (sample+tubes) system. Finally, we estimate the permeability using the incorrect pressure gradient as $kincorrect=ksyst(Lsmpl/Lsyst)$⁠. We found that $kincorrect$ is determined by the $Dsmpl/Dtube, Ltube/Lsmpl, Dsmpl/dsp, Lsmpl/Dsmpl$ ratios which we address in Fig. 6(a). For the justification of our choices, see below.

The expected outcome of the employed model is that $kincorrect$ always underestimates $ksmpl$ [Fig. 6(a)]. Tubing of larger diameter and shorter length, as mimicked by the smaller $Dsmpl/Dtube$ and $Ltube/Lsmpl$ ratios, results in smaller errors due to the smaller pressure drop in the tubes.

The unexpected result is that for the fixed tube diameter $Dsmpl/Dtube$⁠, the increase in $Dsmpl/dsp$ has almost no impact on $kincorrect$ for $Dsmpl/dsp>30$ [Fig. 6(a)], even for very short tubing $Ltube/Lsmpl≪1$⁠. Gradual increase in $Dsmpl/dsp$ implies the reduction of the sample permeability (or increase in the pressure drop inside the sample) as (D_smpl/d_sp)⁻² relative to the empty system. Analysis of the corresponding pressure fields in Fig. 5(c) reveals that (i) for the empty system, the pressure drop occurs only in the tubes [the inclined vs horizontal black dashed line segments in Fig. 5(c)], while (ii) adding a sample introduces a sharp pressure drop at the tube–sample contact, and the magnitude of this drop is mainly defined by $Dsmpl/Dtube$ but not $Dsmpl/dsp$⁠. Increasing the sample length (⁠$Lsmpl/Dsmpl$⁠) at fixed $Dsmpl/Dtube$ reduces the relative contribution of this sharp pressure drop to the full system pressure drop, and $kincorrect$ approaches $ksmpl$ [Fig. 6(a), filled squares]. The existence of this sharp pressure drop also violates the common assumption of equality between the total pressure drop over the system and the sum of pressure drops of its parts (the sample and tubes).

The sharp pressure drop at the tube–sample contact has the following origin. Based on the continuity equation and the divergence theorem, the integral longitudinal flow velocity (or flux) in each transverse cross section is constant along the system. This includes cross sections in the tube region, inside the sample itself, as well as at the tube–sample contact. As such, the fluid entering/leaving the sample from/to the tube has to “squeeze” the constant integral velocity using only a fraction of the pores from the full sample cross section. This leads to a significant energy loss if $Dtube$ is small compared to $Dsmpl$⁠. If the sample has very small pores compared to $Dtube$⁠, still only their fraction $∼(Dtube/Dsmpl)2$ will be used to enter the sample, and the ratio of $Dtube$ to characteristic pore size is not of key importance.

The membrane case addresses the ratio $Dsmpl/Dtube=1$ with a thin sample $Lsmpl/Dsmpl≪1$ placed inside a long tube $Ltube/Lsmpl≫1$ [Fig. 5(d)]. Such ratios are common for the permeability measurements of membranes.^12,48,49 For example, a 100-micron membrane placed between pressure ports 10 cm-apart leads to $Ltube/Lsmpl=100$⁠. Figure 6(a) reveals unsurprisingly a stronger tube length impact for the higher permeability (porosity) samples. Often the tube has variable cross section, as schematically shown in Fig. 5(d) by barely visible gray dashed lines. However, the results in Fig. 6(a) should not be taken for granted in practice: real measurements may deal with the significantly higher (tube diameter)/(pore size) ratios, which we do not address due to the prohibitive computational costs.

Let us qualitatively assess the choice of the parameters describing the geometry in Figs. 5(a) and 5(b). The normalized pressure in Fig. 5(c) always drops from 1 to 0 while fluid passes the tube, tube–sample contact, and the sample (and tube–sample contact with the tubing again). Altering the geometrical configuration redistributes the pressure drops among these individual contributions. Namely, increase in the tube length relative to the rest of the geometry (⁠$Ltube/Lsmpl$⁠) increases its relative contribution to the total pressure drop. If the sample has higher permeability or smaller pressure drop, the relative contribution of the tube to the total pressure drop increases [e.g., the slope of orange dots vs orange diamonds in Fig. 6(a), membrane case 50% vs 80% porosities]. Similarly, increasing the tubing diameter decreases the corresponding relative pressure drop in tubes: compare the slopes of the orange and blue diamonds in Fig. 6(a). Once the tube and sample diameters are fixed (⁠$Dtube/Dsmpl$⁠), their ratio defines the magnitude of the pressure drop at the tube–sample contact [orange vs green vs blue symbols in Fig. 6(a)]. After that, the increase in the sample length (⁠$Lsmpl/Dsmpl$⁠) reduces the relative contribution of the tube–sample contact and tubing to the total pressure drop, and vice versa. This is noticeable, for example, by comparing the orange filled squares + dashed line vs the open squares vs the filled squares + solid line.

The employed simulation setup captures well our initial experimental error, as indicated by the black crosses and green squares nearby [Fig. 6(a)]. Compared to the model with $Dsmpl/Dtube=5,Lsmpl/Dsmpl=6$⁠, experiments have lower influence of tubes: $Dsmpl/Dtube=9 mm/1.89$ mm $≈ 4.8, Lsmpl/Dsmpl=60 mm/9$ mm $≈ 6.7$⁠, and their underestimation is smaller. Figures 5(d) and 6(a) explain why the Fand et al. data in Fig. 3(a) are correct, despite external pressure port connection. They are a membrane case with $Dsmpl=Dtube, Ltube/Lsmpl≪1$⁠, and an incorrectly calculated pressure gradient provides a correct result. It also reveals why Macini et al. underestimated permeability more than our results: these authors have $Lsmpl/Dsmpl≈2$⁠. Their exact $Dsmpl/Dtube$ is unknown, but from a private communication, we assume it to exceed 10. The corresponding experimental setups are shown in  Appendix D.

The current misunderstanding is that connecting the highly permeable tubes (say, 1 mm-diameter) to a low-permeability sample (with micrometer pores) has minimal impact on permeability measurements. This misunderstanding can be traced to the industry standards such as ASTM D5084 (D5856) and D2434. Section 1.4 of ASTM D5084 (or 1.2 in D5856) suggests to perform measurements of highly permeable samples (⁠$k>10−12 m2$⁠) using internal pressure port connections (D2434), while the lower permeability samples can be measured with external ports (D5084). In addition, Sec. 5.2.3 of ASTM D5084-16a (5.2.2 of D5856-15) suggests estimating “parasitic head losses,” which are pressure drops in the experimental system without the sample installed.^50–52 According to Fig. 5(c), the sharp pressure drops at the tube–sample contact appear only when the sample is installed, and cannot be measured without it.

For the sake of completeness, we note that in geotechnical community, porous disks (end pieces) are commonly used to confine a granular sample inside the sample holder. The experimental system without the sample but with porous disks is used to estimate parasitic head losses. It appears that such porous disks may mostly or partially offset the pressure loss at the sample–tube contact we report here. However, such estimations are impacted by the thickness of porous disks (a mesh confining the sample from both sides is actually a thin porous disk) as well as by the ratio of the disk pore size to the sample pore size. Also, a tight placement of each porous disk against the corresponding tube end plays a crucial role.

While it is very difficult to make definitive statements about previous studies with limited available information, in our opinion, some traces of the reported flaw could probably be observed during permeability estimations of ice (Fig. 12 in Ref. 3), sandstone (Fig. 11 in Ref. 53), sintered glass (Fig. B1 in Ref. 32), cement (Figs. 13 and 14 in Ref. 54), fractures (Fig. 7(b) in Ref. 55), and microfluidic device (Fig. 11 in Ref. 33).

Figure 6(b) summarizes some data points from Fig. 6(a) and provides a very approximate and simple empirical correlation between geometrical parameters of the model system and permeability estimation error (⁠$E=ksmpl/kincorrect−1$⁠) for the relatively short tubing length, $Ltube/Lsmpl=0.1$⁠. The empirical correlation takes the form of

for $Dsmpl/Dtube=2,…,10$ and $Lsmpl/Dsmpl=1,…,6$⁠. It is based on the observation that increasing $Dsmpl/Dtube$ leads to the larger error while increasing $Lsmpl/Dsmpl$ reduces it. Please note that this empirical correlation is very approximate and it already fails if $Dsmpl/Dtube=1$⁠.

To test the presented empirical correlation, we employ a commercial equipment for the estimation of the permeability of geological samples. The sandstone sample is installed inside a rubber sleeve [Fig. 7(a)] between flat endcaps with $Dsmpl/Dtube≈24$⁠, and then it is exposed to the confining pressure. The equipment operates with uncommonly long samples (⁠$Lsmpl/Dsmpl≈8$⁠) and allows simultaneous connection of the five internal pressure ports, in addition to the external ones [Fig. 7(a), schematic]. Figure 7(c) shows the brine permeability values calculated using different pairs of the pressure ports. We measured pressure drop for the empty equipment (without the sample) and use these values to correct the permeability obtained from the external pressure ports.

On average, permeability from the internal pressure ports 1–5 overestimates the external one by +56% (or by +93% without the empty equipment correction), Fig. 7(c), red vs black crosses. The empirical correlation suggests $0.5(Dsmpl/Dtube)(Dsmpl/Lsmpl)=0.5·24/8=1.5$ or permeability overestimation by $+150%$⁠. During installation of the solid (non-granular) sandstone sample, we noticed a visible gap between the sample and one of the flat endcaps, Fig. 7(b), bottom. The gap width is comparable with a 1.6 mm endcap tube diameter. The gap results in the reduction of the pressure loss at tube–sample contact and lower permeability error. However, it also has the opposite effect: taking another sample with better base alignment against the endcaps will lead to the larger error. This also suggests a simple scenario for the scatter in permeability of similar sub-samples extracted from a large homogeneous sample, measured on the same equipment.

Here, we consider the base case of a uniform hollow cylinder with a constant cross section attached directly to the sample. In practice, it is common to see deviations from this case. Namely, Fig. 8 shows a part of laboratory equipment (Hassler cell⁵⁶) used for permeability estimation of geological samples, which has grooves on its surface to enable more uniform inflow. (Please note that the experimental setup in Fig. 7 has flat endcaps without grooves.) Our laboratory experiments on rock samples suggest that adding area of the grooves to the tube area allows (very) approximate estimations of $Dsmpl/Dtube$ and the corresponding permeability error. In some fields, it is common to use porous disks to confine a granular sample inside the holder and to have gaps between the tube openings and the disks. The ends of a solid sample can be nonparallel to the surface of the holder end pieces with tubes, tubes can have non-uniform cross section, and so on. Such modifications will affect $ksyst$ and $kincorrect$⁠. However, what will stay unchanged is the measurement of permeability of a virtual geometry with a constant cross section over its entire length [the red in Figs. 1(c) and 5(a)].

The flaw affecting permeability measurements originates from the external placement of the pressure ports and calculating incorrect pressure gradient using length of the sample but not of the system. This can lead both to correct and incorrect permeability estimations, depending on the sample-to-tube diameter ratio $Dsmpl/Dtube$⁠. Larger $Dsmpl/Dtube$ can lead to significant pressure loss at the sample–tube contact and to permeability estimation error. The pressure loss holds for the samples, in principle, with any pore size compared to the tube diameter [Fig. 5(a)]. The error magnitude can be very approximately estimated using $ksmpl/kincorrect−1=0.5(Dsmpl/Dtube)/(Lsmpl/Dsmpl)$ empirical correlation. We posit that this pressure loss may have significant impact not only on the single-phase permeability but also on multi-phase flows measurements, such as capillary pressure estimation.

The flaw can be addressed in several ways. One way is to use internal pressure ports [Fig. 3(b)], but they are more difficult to engineer, and they capture only a part of the sample geometry, which may increase the scatter in permeability estimation. Another way is to have identical cross-section area of the sample and tubes, placing the pressure ports close to the sample, and calculating incorrect pressure gradient using $Lsmpl$⁠. Introducing gaps, intentionally or unintentionally, between the tube opening and the sample end surface reduces the pressure loss at the sample–tube contact. It appears that the gap impact will be the strongest for $ksmpl≪ktubes$⁠. Grooves on the surface of endcap in Fig. 8 act similarly to gap. Increasing the sample length (⁠$Lsmpl/Dsmpl$⁠) also helps, as suggested by Fig. 6(b). One more solution is to follow the strategy described here: to mimic the measured geometry (sample itself, confining container, tubes) in simulations and extract $ksmpl$ from the subsequent analysis. The last solution is time-consuming but can be used to correct the already measured permeabilities if sufficient experimental details are available.

Connecting reality with its virtual counter-part using only first principles is a very time-consuming process. However, at the same time, this connection reveals the full predictive power of computer models that improve understanding and control of real-life processes around us.

In the current study, we operate with dimensionless permeability, making our results quantitatively universal as long as Stokes physics coupled with no-slip boundary condition is valid. We assume that deviations from it, such as occurrence of inertial or Klinkenberg effects, will not change qualitatively the reported results.

This study illustrates the importance of having independent assessments for such scalable and non-trivial quantity as permeability: using only one method, either experiment or simulation, leads to the lack of clarity, which may persist for decades (not to say, centuries). We believe that the presented results will help to remove the discrepancies between gas and liquid permeability measurements, between falling and constant head methods, between permeabilities of similar samples measured by different groups, or between subsamples of different dimensions extracted from a larger homogeneous sample. The starting point is to keep in mind the virtual geometry of constant cross section over the full system length [the red in Figs. 5(a) and 1(c)].

S.K. acknowledges Alexander Kurilik and Vasily Belokhin for numerous discussions on the fundamentals of computed tomography; Evstigneev Nikolay and Dmitry Logashenko for discussions on theoretical aspects of the Stokes equation; Kenneth Kennedy and the cleanroom staff of the KAUST Nanofabrication Core Lab for co-design and fabrication of the silicon meshes; Ernest Neil Davison for accurate glass blowing of the sample tubes and KAUST Workshop for fabrication of PVC plugs and Teflon CT holder; Marijn Boon, Manuel Dierick, and Wesley De Boever (XRE Tescan) for discussions on the usage of XRE Tescan hardware and software; and Paolo Macini for sharing details on their experimental study. We are grateful for the allocation of computational resources by the Supercomputing Laboratory at King Abdullah University of Science and Technology (KAUST) in Thuwal, Saudi Arabia.

This work was supported through funding to the Ali Al-Naimi Petroleum Engineering Research Center (ANPERC) at KAUST.

The authors have no conflicts to disclose.

S.K. initiated the study, designed experiments using equipment available at KAUST, surveyed the literature, identified the flaw, and performed experiments, CT scanning, image processing, and LBM simulations. S.K. wrote the manuscript and prepared figures. M.Y. introduced and consulted S.K. on flow-related laboratory facilities at ANPERC, performed permeability measurements on rock samples, designed and assembled the first experimental setup in room1 together with S.K.; suggested the flow system degassing approach, participated in technical discussions, corrected and discussed the manuscript with S.K., and proofread the manuscript. T.P. supervised and funded the project, provided feedback and corrections, and edited the paper.

Siarhei Khirevich: Conceptualization (lead); Visualization (lead). Maxim Yutkin: Methodology (supporting). Tadeusz W. Patzek: Funding acquisition (lead); Supervision (supporting).

3D gray-scale reconstructions and binary (segmented) images: https://doi.org/10.25452/figshare.plus.16821412

2D projections and 3D gray-scale reconstructions: https://doi.org/10.25452/figshare.plus.16850653

Pressure and viscosity measurement logs, experimental permeability values, 2D optical scans of beads: https://doi.org/10.25452/figshare.plus.16867417

The appropriate choice of bead diameter played a key role in this study. There are two competing effects: discretization (scanning) resolution and Reynolds number, $Re=ρdspvavg/μ$⁠, where $dsp$ is the bead diameter, $vavg$ is the interstitial flow velocity, μ is the fluid viscosity, and ρ is the fluid density. Decrease in $dsp$ reduces both the discretization resolution (in voxels per $dsp$⁠) and Reynolds number, and vice versa. Our original plan was to obtain the maximum resolution of r = 100 voxels per bead diameter and Re < 1, ideally $Re≪1$⁠. The available CT scanning equipment (see below) advertises the minimum voxel size of 3 μm, but taking into account the usual contrast loss at the minimum declared voxel sizes, we targeted the voxel size of 5 μm or larger. This leads to the average bead diameter of 500 μm. Using water as working fluid (μ = 1 mPa·s at 20 °C), such a diameter choice results in the pressure drop of ∼300 Pa for $Re=1$ when using the sample dimensions and permeability from this study. The pressure transducer we used has the working range of 0–34.5 kPa (or 0–5 psi), and it is already the lowest range in the product line of up to 10 000 psi. This simple estimation suggests that water was not suitable as the working fluid for the available equipment. The only choice left is to use another fluid with higher viscosity. Increase of viscosity is an efficient tool to reduce Reynolds number: for a given pressure drop, $Re∼μ−2$ because increase in μ simultaneously reduces $vavg$ and $Re$⁠. We chose glycerol–water solutions as the working fluid for two reasons. First, glycerol has excellent solubility in water, allowing us to switch between different glycerol–water solutions by flushing an “old” solution with a “new” one. We also could replace the dissolved glycerol with water, then dry the sample, and obtain air–glass phases for better CT imaging. Second, compared with water, glycerol has a significantly higher viscosity (∼1400 mPa·s at 20 °C). The only drawback is its strong viscosity dependence on temperature, about 10 times higher than that for water (about 10% change in μ for 1 °C at 20 °C). Therefore, we mixed glycerol with water to obtain the viscosities of only about 20–50 mPa·s to balance the lower Reynolds number and weaker viscosity–temperature dependence. The bead diameter was fixed at ∼500 μm, and this resulted in the actual $Re$ values between 0.0077 and 0.26 in this study.

We used two commercially available glass beads, provided by Kramer Industries (Piscataway, United States) with the US mesh 30–40 grade (beads1) and “Cell disruption media, 0.5 mm” by Scientific Industries (Bohemia, United States) referred to as beads2. Both types of glass beads are made of soda-lime glass, as reported by manufacturers.

For a two-dimensional characterization of the beads, we spread them in a single layer over an adhesive tape and acquired high-resolution images using an optical microscope DM6 FS by Leica Microsystems (Wetzlar, Germany), Figs. 9(a)–9(c). The resulting images were processed with the circle finding method by Atherton and Kerbyson,⁵⁷ implemented in the function imfindcircles in the MATLAB (Mathworks, Portola Valley, United States) software package, Fig. 9(d). Each adhesive tape fragment contained about 4000 beads diameters. We collected four fragments for each beads type, and the corresponding beads diameter distribution are shown in Fig. 9(e). Figures 9(b)–9(d) suggest that beads1 are less homogeneous and contain more defects than beads2, which is also noticeable in the larger distribution width and skewness [Fig. 9(e)], lower density values [Fig. 10(c)], and the sample color change after CT scans [Fig. 12(g)]. The beads1 diameters below 350 μm were excluded from the Sauter diameter calculation because they are not separate beads [red circles in Fig. 9(d)]. We again note that particular choice of Sauter diameter as characteristic length has no impact on experiment–simulation permeability comparisons within this study, and its main purpose is to compare the present results with past studies [Fig. 3(a)].

Bead density was determined using the Archimedes law [Fig. 10(a)],

where $Wair,wat$ is beads weight either in air or de-ionized water (reference ultrapure water (type 1) from the Milli-Q Purification System, Millipore, Merck, Kenilworth) of the corresponding density $ρair,wat$⁠. We used the values of $ρair=0.0012$ and $ρwat=0.997 66 g·cm−3$⁠. For some measurements, we observed a slow change in $Wwat$ [Fig. 10(b)] that we attribute to the dissolution of bubbles or trapped air film on the surface of a plastic Petri dish (with or without beads) after its immersion. Therefore, for the density calculations, we took the average value of $Wwat$ after it was observed stable for at least 12 h. Weights $Wair,wat$ were determined using two analytical balances with a 0.1 mg readability, XS204 (room1), and XSE104 (room2), by Mettler Toledo (Columbus, United States). Each density measurement involved about 2–5 g of beads, which is 12 000–30 000 beads assuming 0.16 mg weight of a single 500 μm bead. Each 2–5 g sample was disposed of after the measurement. Weight oscillations in Fig. 10(b) observed for all measurements in room1 are due to temperature (and water density) oscillations [Fig. 18(i)]. These oscillations give an idea about accuracy of the employed analytical balances relative to the detected weight change. We determined and used the values of $ρbeads1=2.490$ and $ρbeads2=2.515 g·cm−3$ in all calculations. Normalized standard deviation of density is ±0.1% for beads1 and ±0.05% for beads2. We do not consider this uncertainty in the calculation of total uncertainties ( Appendix E) because the porosity uncertainty mainly depends on the uncertainty in confining container volume, originating from the voxel size uncertainty.

To confine glass beads, we used glass tubes made by Schott (Mainz, Germany) from borosilicate glass. The tubes have the reported density of $2.23 g·cm−3$⁠. This is noticeable on CT images as darker gray value of wall voxels compared to solid voxels of $2.5 g·cm−3$ soda-lime glass beads [Fig. 21(e)]. One tube end was welded [Fig. 11(a)] and has a sightly wavy shape. The sample with packed beads did not overlap with the welded end. For the porosity calculation, we need to know the internal tube diameter. The manufacturer specifies only the outer diameter as 12 ± 0.16 mm and the wall thickness as 1.5 ± 0.07 mm, resulting in 8.7–9.3 mm (600 μm) range for the possible inner tube diameters. Such uncertainty is totally unacceptable for this study, and, therefore, internal tube diameters for each sample (P3 or P4) were determined from CT images. We note that 600 μm range is applicable for different samples, while each particular tube sample demonstrated excellent internal diameter uniformity over its length as well as roundness, at least on a ∼10 μm scale relevant for this study. For each packed tube sample, several gray high-resolution CT images with voxel sizes between 9.8 and 7.3 μm were used to manually find a circle fitting the glass tube wall. The obtained diameters were averaged to estimate the internal diameter for each tube. Our estimations resulted in $DP3tube=8984$ and $DP4tube=9044$ μm. In a similar way, the length of each tube segment between plugs was estimated from CT images as $LP3tube=59 649$ and $LP4tube=60 330$ μm. For uncertainties in these values, see below.

Plugs [Fig. 11(b), top] were fabricated from polyvinyl chloride by a KAUST workshop, according to a custom design. Four rubber O-rings restricted leaks between plug and glass tube and also oriented the plugs coaxially with the tubes. Each plug has a hole of about 1.89 mm in diameter determined from CT images with the 6-micron voxel size (the 3D images are available in the supplementary material). Initially, we used nylon mesh to prevent the beads from entering the plug hole. Later, we switched to silicon meshes for better segmented image quality [Fig. 11(b), bottom].

Glass beads were packed using tools shown in Fig. 12(a). In general, we followed previous studies⁵⁸ to reduce the packing porosity as much as possible in order to obtain the mechanically stable samples. The goal here is to avoid beads displacement during permeability measurements or sample movements within lab. We packed beads in water under ultrasonic vibration [USC-THD/HF Ultrasonic Cleaner by VWR (Radnor)] and added beads in batches [Fig. 12(d)] with about 10 batches in total, for P3 or P4. After adding the next batch of beads into the glass tube, we applied vibration for about 30 s [Fig. 12(e)] and also used glass stick to gently tap the beads from the top [Fig. 12(f)]. The necessity of the tapping step is unclear, but we kept it for its simplicity. USC-THD/HF Ultrasonic Cleaner has the reported ultrasonic frequency of 132 kHz, which significantly exceeds the frequency used by Yu et al.⁵⁸ Plastic Petri dish [Figs. 12(b) and 12(c)] was used to track the beads accidentally passing through the mesh, but none were observed. Original beads1 and beads2 with the corresponding packed samples are shown in [Fig. 12(g)]. Color change was observed later for both beads types due to extensive x-ray CT scanning, and beads1 tends to change their color more due to larger fraction of defects (test sample was exposed to x-rays for shorter time compared to P3). The porosity of both samples was determined using mass of the packed beads, their density, diameter of each glass tube and its length between plugs. The tube was assumed to be an ideal cylinder, with the top and bottom bases shown in insets in Fig. 24 as the “left” and “right” meshes. The half-mass of ∼50 beads crossing each cylinder base in the small gap between plug and tube wall was estimated. This mass (0.12–0.14%) was subtracted from the total mass of glass beads for samples P3 and P4. After that, the volume of glass beads was calculated and increased by the approximate volume (∼0.5%) of the silicon meshes for P3 and P4 samples. This allowed us to obtain a single value for the volume of glass (either beads or meshes) within P3 or P4, and segment air from glass phases applying a global threshold to gray CT images. The porosities were estimated as $εP3=35.414¯$% and $εP4=35.689¯$%, where underlined digits indicate uncertainty (discussed below).

Porosity distributions along P3 and P4 samples are shown in Fig. 24. At both ends of each sample, there are oscillations originating from ordering of the beads near a hard wall (the so-called geometrical wall effect, e.g., Figs. 1 and 2 in Khirevich et al.⁵⁹). One can see that at the welded tube end, which was at the bottom during the packing process [Fig. 12(b)], the porosity oscillations have lower minima. This suggests that the welded end of each sample is densely packed.

Let us assess the obtained porosity value of ∼35.5% for the tube-to-bead-diameter ratio of about 18. The similar porosity for the similar diameters ratio was observed by Yu et al.⁵⁸ (their Fig. 1, bottom) for samples packed in batches using one-dimensional vibration. The value of $εP4$ exceeds $εP3$ due to the narrower distribution of beads2 diameters,⁶⁰ Fig. 9(e). All porosity values are below so-called “random close packing limit” of ∼36.6% reported for packings of equal spheres,^60,61 which supports the assumption about mechanical stability of the prepared samples. From the operational experience, there was no evidence of beads displacement within the sample.

Despite the packing process with immersion, all three samples contained some air bubbles inside, visible through the transparent confining wall. Placing the sample under vacuum allowed dissolved air and bubbles to escape the sample through the holes in plugs.

Each sample was connected to a syringe installed on the PHD ULTRA 4400 pump (Harvard Apparatus, Holliston, USA) using flexible tubing, Fig. 13. The internal diameter of the tubing was estimated from thin slices using optical microscopy (1.3 μm pixel size) and x-ray CT imaging (6-μm voxel size). Both methods resulted in 1.73 mm, which is noticeably different from the specified 1.65 mm (0.065 in.). The flexible tubing length was estimated using the digimatic caliper 500–197–30 by Mitutoyo (Kawasaki, Japan) with the reported maximum permissible error of ±0.02 mm (±0.04 mm). However, taking into account the non-straight tube cuts and slight bending, we use the subjective uncertainty estimation of ±0.5 mm, assuming that the reported system length is accurate to within 1 mm. The internal tube end location inside each plug was found from CT images with the 6-micron voxel size. The total system length for S3 and S4 of about 211 mm results in the relative uncertainty of ±0.24%.

We employed four glass gastight syringes [Fig. 16(a)]: two 10 ml by SGE Analytical Science Pty Ltd (Ringwood, Australia), one 1000 series 10 ml, and one 1000 series 25 ml by Hamilton (Reno). The rationale here was to address possible flow rate variation due to the variability in dimensions of each syringe. Each series of 16 permeability measurements for each viscosity value and each P3 or P4 sample started with the 10 ml syringe until its depletion. Hereafter, we switched to the 25 ml syringe refilling it when necessary to finalize the collection of 16 permeability values. Each 10 ml syringe was used for a single viscosity value, while the 25 ml syringe was employed for all three viscosities, as indicated in Fig. 16(a). Gastight syringes by both Hamilton and SGE have the declared accuracy of $±1%$ of the dispensed volume, while for the PHD ULTRA 4400 pump, the accuracy is $±0.35%$⁠. We tested all four syringes with this pump and found an average relative error of $±0.2%$ while comparing the specified and actual dispensed volumes of de-ionized water. The volume was calculated by weighing the dispensed sample on XSE104 balance and assuming the water density of $0.99 766 g·cm−3$⁠. The infusion time was between 20 and 240 s, and no systematic dependence of the error in the dispensed volume on time was observed. Therefore, we consider the dispensed volume error of $±0.2%$ as the total uncertainty in the discharge rate during permeability measurements.

Inlet pressure $Pin$ was measured using the gauge-type DXD pressure transducer (Heise, Stratford, USA) with a working range 0–5 psi (0–34.5 kPa) and averaged its values as indicated in Fig. 14(a) during non-zero flow. The transducer was connected to tubes, and the corresponding hydrostatic region [dashed green in Fig. 1(c)] was filled with the Fluorinert Electronic Liquid FC-40 oil by 3M (Saint Paul) with the declared density of $1.85 g·cm−3$⁠. For each permeability measurement, the operator opened the right door of the glass box [Fig. 18(f)], turned on the pump, started infusion, kept the right door almost closed while waiting for infusion to finish, opened the door again, turned off the pump, and locked the door closed. The vibration from door opening–closing and from turning on–off the pump is visible as spikes in the recorded pressures [Fig. 14(b)].

In the absence of flow, the transducer recorded pressure P₀, Fig. 14(b). We found that the value of P₀ is determined by the shape of meniscus at the system outlet: a concave meniscus reduces P₀, while a convex one increases it [subpanels “i” and “ii” in Fig. 14(b)]. Figure 14(b), subpanel “iii,” shows also the situation when the drop forming a convex shape reduces P₀ due to gravity. The reported accuracy for the DXD pressure transducer is 0.02% of the full scale, meaning that it is able to resolve pressures at least 34.5 kPa/5000 $≈$ 7 Pa while the capillary pressure for 1.73 mm inner diameter tube and the employed glycerol–water solutions is about 160 Pa. Wiping off the drop “iii” [Fig. 14(c)] results in an unpredictable increase or decrease in P₀ [Fig. 14(d)], depending on the shape of the formed meniscus.

Figures 14(e) and 14(f) show P₀ before and after each permeability measurement for all 96 experiments. We observed that after many hours at rest, the system tends to form meniscus “i” at the outlet, which reduces P₀ by $∼160$ Pa and explains why left parts of Figs. 14(e) and 14(f) show lower values. Figures 14(e) and 14(f) also includes P₀ values obtained for the short time intervals, $∼20$ min, between consequent experiments. Here, P₀ before and after a given measurement may be similar. Additionally, the value of P₀ is also affected by horizontal leveling of the system, which could slightly change every time the syringe was refilled. Therefore, we took $P0=1.3$ kPa and considered its uncertainty to be ±0.05 kPa [the gray rectangles in Figs. 14(e) and 14(f)]. The pressure drop $Pin−P0$ was always between 7 and 37 kPa or 22 kPa on average, and, therefore, we calculated the relative uncertainty in pressure estimations as $±0.05/22≈±0.2%$⁠, exceeding the uncertainty of DXD pressure transducer tenfold.

In order to address the impact of horizontal alignment of the system on the obtained permeability results, we inclined the sample significantly up and down, which changed P₀ to ∼1.7 and ∼0.9 kPa, respectively. Unsurprisingly, the permeability values estimated for the horizontal and both inclined cases were identical when $Pin−P0$ was used.

Initially, we installed the sample, syringe pump, and rheometer next to each other [Figs. 18(b) and 18(c)] in order to maintain similar temperature for all of them. However, this led to ∼15% oscillations of the measured permeability values, as demonstrated by the black circles in Fig. 15. Such oscillations are typical for permeability measurements, and they also can be observed in previous studies [Fig. 12(b) in Ref. 28, the data by Loudon⁴⁶ or James and McLaren⁴⁵ in Fig. 3(a)]. We identified the following sources for the permeability oscillations: scatter due to the viscosity measurements, ambient room temperature, and viscous dissipation. Below, we discuss these sources and provide the ways for their minimization, resulting in ∼1%-stable permeability values (green and red dots in Fig. 15).

We mixed 99.5% glycerol (Fisher BioReagents, Fisher Scientific, Waltham, United States) and DI water in approximate mass ratios 119:51, 76.5:23.5, and 80:20 [Fig. 16(a)], resulting in solutions with viscosities of $μ0≈21, μ1≈38, and μ2≈55$ mPa $·$ s. For each sample P3 and P4, we mixed $μ0,…,2$ separately and, therefore, obtained six water–glycerol solutions with $μ0,…,2P3≈μ0,…,2P4$⁠. Viscosity of each solution was determined using MCR 302 modular compact rheometer by Anton Paar (Graz, Austria). We used two measurement systems together with MCR 302: “double gap” (DG) and “concentric cylinder” (CC), shown in Fig. 16(b). The DG system requires 5 times less sample volume and has larger surface contact area allowing for a more accurate torque measurement. However, tests with DI water revealed larger scatter for this system. At low shear rates, this scatter increases probably due to the dominant surface tension effects at the sample free surface (e.g., Figs. 1 and 5 in Johnston and Ewoldt⁶²). With increasing the shear rate, the relative contribution of free surface tension to the viscosity values measured in the DG system decreases. Also, increasing shear rate both for the DG and CC leads to an increase in apparent viscosity values due to transition to the nonlinear flow. For the CC, this transition occurs at significantly lower rates probably due to the developing nonlinearity in flow behavior at the system's bottom. Another source of the scatter in viscosity values for the DG system is the uneven fluid levels in the inner and outer gaps. This level gap forms when the rotating bob enters the cup filled with a fluid sample. It is reported in DG specification and appears to be inherent to the design of this measurement system.

Initially, we collected the fluid samples directly at the outlet, Fig. 16(c). In room1, where the rheometer is located, the effluent viscosity was measured without active temperature control of the rheometer, assuming that all parts of the experimental setup were in thermal equilibrium. In room2, the temperature values were recorded by a digital thermometer Fluke 54 II B (Everett) and then used as an input to the active rheometer temperature control. Working fluids were switched, and the flushing quality was approximately estimated from the pressure drop change. By doing so, we attempted to measure the viscosity of a fluid emanating from the porous sample and for which the pressure drop was just measured. More careful analysis revealed that the sample pore volume was about 1.4 ml, and we needed to collect about three pore volumes to fill the DG system sufficiently (3.8 ml). Filling the CC system would require about 14 pore volumes. This approach requires significant fluid sample consumption and contains an unclear contribution from the old solution, at least for the first measurements after solution change. The results for the DG viscosity values and the lack of proper sample washing are indicated in Fig. 15 as “control of viscosity measurement: no.”

Later, for proper sample flushing control, we switched to a different approach. First, the viscosity–temperature dependence for each solution was determined using the CC system after taking 19 ml samples from a stock solution of a given viscosity. Next, the corresponding linear calibration curve was found (see the next paragraph). Then, the sample was completely flushed with 50 ml of a “new” viscosity solution. To accelerate diffusion-controlled mixing during flushing, a heating jacket [50 °C, Fig. 16(d)] was mounted on the sample. Assuming glycerol and water diffusion coefficients to be at least $∼10−10$ m²/s,^63,64 and the porous sample radius of 0.0045 m, the time needed to diffuse from the sample wall to its center is about $0.00452/(2×10−10)≈105$ s or 28 h. In practice, we found that continuous infusion of 24 ml of a new solution during 16 h was already sufficient to obtain stable permeability values. However, we also added four cycles of 4 ml infusion at 1 ml/min with 2 h of soaking to ensure complete flushing. The corresponding permeability values are referred to as “control of viscosity measurement: yes” in Fig. 15.

Figures 17(b)–17(d) show the viscosity–temperature dependencies for the “control of viscosity measurement: yes” case determined at four temperatures (21.0, 21.5, 22.0, and 22.5 °C). For each solution, we took three 19 ml samples from the stock solution [Fig. 16(a)]. Each obtained viscosity–temperature curve was fit with a line to obtain a calibration curve, and later, this curve was used to calculate the viscosity values at the temperatures from a digital thermometer attached to the sample [Figs. 1(a) and 1(c)] inside room2. The digital thermometer readings were offset from the MCR302 temperatures by 0.1 °C, for which we corrected. Compared with DI water, the viscosity values for glycerol–water solutions obtained with the CC system show reduced scatter [panels (a) vs (b)–(d) in Fig. 17]. We relate this to the lower contribution of the free surface tension to the measured viscosity values. Also, inertial effects reduced, and the glycerol–water viscosities show no difference between shear rates of 20 and $50/$s. For DI water, the difference was about 5%. Surprisingly, glycerol–water viscosities measured with the DG system still demonstrate ∼5% scatter, despite reduced contribution from the free surface tension. The average DG value is below the CC viscosity values. We attribute this scatter to the aforementioned uneven levels of the fluid sample in the inner and outer gaps. This explains why, during the stabilization of permeability values in Fig. 15, the average value between “control of viscosity measurement: no” and “yes” increased. For the final permeability measurements, we used viscosities obtained from linear fits of the CC values at the shear rate of $20/$s. All logs from the MCR302 rheometer are available as the supplementary material.

All viscosity measurements were performed for torque values of $∼10−6$ N·m, which is well within the MCR302 torque operating range of $10−8,…,10−1$ N $·$ m. We did not find explicit values for the accuracy of the MCR302 rheometer, probably due to the variety of scenarios this device and its measurement systems can be used for. However, the calibration procedure (while using a viscosity standard 2000 AW, ν = 1459 mPa $·$ s at $20 °$C) suggests to use ±1.5% as the maximum permissible error for the shear rates of ∼10/s. For our purpose, this error seems to be reasonable taking into account the CC results for DI water [Fig. 17(a)], which is more challenging compared to glycerol–water solutions. Therefore, we assume that the error for the measured viscosity values is ±1.5%.

For the sake of completeness, we estimate the shear rate within sample pores, $γpores$⁠, during permeability measurements using equation for the shear rate in pipe flow (⁠$γpipe=8vlinear/dpipe$⁠, where $vlinear$ is the linear fluid velocity and $dpipe$ is the pipe diameter). If we assume the pore diameter to be 1/3 of a bead diameter or 150 μm, and the linear flow velocity of 1 mm/s, $γpores≈$ 50/s.

Our initial experimental setup, referred to as “room1,” is shown in Figs. 18(a)–18(c). We placed all key parts (syringe, pump, and sample) close to the rheometer, where viscosity of the working solution was measured without rheometer temperature control. Our assumption was that the temperature of the rheometer and experimental parts will be similar. The obtained permeability values exhibited the ∼15% oscillations (Fig. 15), quite typical for standard permeability measurements. The experimental setup in room1 was directly exposed to the Heating, Ventilation, and Air Conditioning (HVAC) system shown in Fig. 18(c). Using an undocumented feature of the DXD pressure transducer to measure ambient temperature $TDXD$⁠, we found periodic oscillations of about $1 °$C in amplitude [red in Fig. 18(i)]. These oscillations were also noticeable in the weight–time dependencies during the glass-bead density measurements in room1, Fig. 10(b). Later, we installed a glass box [Figs. 18(f)–18(h)] in another part of the laboratory with different temperature control and without periodic temperature oscillations. Metal surfaces of the glass box (ceiling and rear wall) and the front doors exposed to HVAC were covered with additional thermal insulation to damp ambient temperature oscillations. Tracking $TDXD$ values inside this glass box during 12 h and 14 days revealed a complex temperature–time dependence without clear patterns, see the blue dots in Figs. 18(i) and 18(j). However, sometimes we observed time intervals of several hours, when the temperature was stable within $0.05 °$C [Fig. 18(k)]. This motivated us to measure temperatures of individual parts of the experimental setup in order to estimate their equilibration time with ambient temperature.

We employed the digital thermometer Fluke 54 II B with two thermocouple junctions referred to as “T₁” and “T₂” (or “T_i”) in this paper. Compared with $TDXD$ with $0.01 °$C readability, the digital thermometer values are only $0.1 °$C readable, but this was found to be sufficient for the purposes of this work. In addition, $TDXD$ allowed us to track the temperature readings remotely and log them into a file, while the digital thermometer lacked this capability. First, we measured air temperature inside the glass box with T₁, T₂, and $TDXD$ (top row in Fig. 19). The values of T_i were found to be almost identical except one time interval between 300 and 500 min, where T₁ was 21.7, while T₂ was 21.6. We attribute this to the ambient temperature of 21.65 and slightly different location of the thermocouple junctions. Temperatures $TDXD$ were found to be ∼0.45 °C above T_i due to the lack of calibration (recall that $TDXD$ denotes temperature values from the pressure transducer, which is not a thermometer). Also, $TDXD$ values follow the ambient air temperature with a slight time delay relative to T_i because of the casing of the pressure transducer.

We placed the thermocouple junction T₁ inside a dummy sample or inside an installed syringe, while T₂ was left in air, as shown in “sample + air” and “syringe + air” rows in Fig. 19. When the temperature gradients were not very steep, T₁ approximately followed T₂ within ∼60 min (the gray rectangles with the indicated time in minutes). Such a temperature behavior suggests that if the air temperature inside the glass box is $0.05 °$C-stable within 60 min or more, the syringe and the sample will also have a similar temperature. Therefore, we measured permeability after ambient air temperature $TDXD$ within the glass box was $0.05 °$C-stable for at least 60 min. The corresponding permeability values are named as “control of ambient temperature: yes” in Fig. 15.

Due to the simplicity of implementation, we also placed T₁ inside the sample while T₂ on its surface. A thermocouple junction on the surface would allow us to track sample surface temperature during the experiments, and the corresponding temperature readings demonstrated a time lag between 25 and a few minutes (“sample + surface” rows in Fig. 19). To our surprise, during actual permeability measurements, we found that, despite $0.05 °$C-stable ambient temperature for > 60 min, both T_i installed on the sample surface showed 0.2–$0.3 °$C temperature increase within 1–2 min after the experiment start (Fig. 20). This increase was observed in all experiments, and T₁ was never below T₂. T₁ was installed near the entry point of fluid into the porous sample while T₂ at the sample end, as indicated in Fig. 1(c). Swapping T₁ and T₂ did not change this picture. For all employed flow rates, the increase in T_i never exceeded $0.3 °$C, but for lower flow rates, the temperature increase was slower. We attribute this temperature increase to viscous dissipation, which is the irreversible conversion of the fluid kinetic energy to heat. The rate of heat generation due to viscous dissipation is proportional to the first power of viscosity and the second power of the flow rate.^65,66 In a test experiment, we used two discharge rates, 1 and 12 ml/min. This rate increase translates into 12 times larger flow velocity and fluid enthalpy efflux, and ∼150 times larger heat generation rate. However, both experiments resulted in $0.2 °$C temperature increase after 1 and 2 min for 12 and 1 ml/min, respectively (Fig. 20). It seems that the generated heat increases the sample temperature and the resulting heat and enthalpy fluxes from the sample into the air. These fluxes cool the sample and counteract the viscous dissipation heating. Therefore, the sample always reaches some temperature above ambient air to establish a constant heat flux and dissipate the generated heat. In other words, for all employed flow rates, and given the sample geometry, the viscous dissipation always heats the sample until its surface reaches 0.2–$0.3 °$C above ambient air temperature. Having this hypothesis in mind, we decided to reduce the duration of permeability experiments from several minutes to 30 s while collecting pressure logs in a high-frequency mode (about 40 readings per second; in this mode, the DXD pressure transducer does not allow recording ambient temperature). From below, the duration of each experiment was limited by pressure oscillations due to drop formation and detachment at the outlet [Fig. 14(a)], and the necessity to log several such oscillations for averaging.

We used two time intervals between the consecutive experiments: about 20 min (short) and above 60 min (long). After performing each experiment, $TDXD$ was monitored during a short interval, and if it was $0.05 °$C-stable, we started the next permeability measurement. We never used more than two consecutive short intervals even if the ambient temperature remained stable. It appears that some excess heat was still in the sample, and we needed to cool the sample down while waiting for a long interval. Such an approach resulted in excellent stabilization of permeability values (“control of viscous dissipation: yes” in Fig. 15). Each series of 16 experiments in Fig. 4 for a given glycerol–water solution and sample was performed in a row, without observing outliers. The actual wait times between each 30 s experiment are specified in Figs. 4(a) and 4(b). The DXD pressure transducer logs for all 96 experiments are available in the supplementary material.

The bead packing was flushed with DI water and then dried to obtain air–glass phases inside the sample. For image acquisition, we use the CoreTOM X-Ray CT scanner (XRE Tescan, Ghent, Belgium). Figures 21(a) and 21(b) show the scanner, x-ray source, sample, rotating stage, and detector. The sample installed on the rotating stage is penetrated by x-rays from a number of angles resulting in 250–1800 of 2D projections [Fig. 21(d)], which are recorded with the detector [Fig. 21(b)]. After application of the inverse Radon transform, 3D image of the sample [Fig. 21(e)] is reconstructed from its 2D projections. This last step is implemented using the proprietary Tescan software.

Each sample P3 or P4 was scanned at 18 different resolutions, r, between $rmin≈3$ and $rmax≈64$ voxels per Sauter sphere diameter of beads1 or beads2. The corresponding voxel sizes were between 180 and 7.3 μm. It was found that, for the employed system, better image contrast is observed for the source–detector distance (SDD) of about 850 mm, within the possible range of ∼400–1000 mm. For example, changing SDD to $≈ 400$ mm substantially reduced scan time but also degraded image contrast [Fig. 22(a), top]. The tube voltage of 60 kV and power of 15 W were found to be optimal for the employed x-ray system while looking for the balance between scan time and image quality. The exposure time was about 4.5 s per projection to suppress image noise to below the level at which the use of noise reduction algorithms is unnecessary. No beam-hardening correction was applied, and no ring artifacts were removed.

When the discretization resolution exceeded ten voxels per bead diameter, the whole sample scan must be split into several segments [Fig. 21(c)], and for the highest resolution, 8 such segments with slight overlap were obtained. Later, these segments were merged together based on motor positions. We found very minimal visual misalignment between the segments, with no impact on the simulated permeability values. No image brightness normalization between the segments was necessary due to the identical scanning conditions.

During long scans, the state of x-ray source (tube) changes, and the tube must be refreshed. The scan process itself is 360° of continuous rotation of the sample at a constant angular velocity, which has angular acceleration only at the beginning and end of the scan. Interruption of the scan process may result in the slightly misaligned sets of projections. We found that the image quality degraded if a single 3D reconstruction is obtained using projections from different sets. In other words, if the tube is refreshed at least once during a given 360° scan, image quality decreases. In addition, we found that the tube needs to be refreshed after ∼4 h of scan (otherwise image quality again decreases). Four hours per single scan of the whole sample (or each segment) is, therefore, the upper time limit. Taking into account the choice of other settings (the tube voltage, power, and exposure) and difficulties in the allocation of CT scanner time, the originally planned $rmax=100$ scan was reduced to $rmax≈65$⁠. All 18 resolutions for each sample P3 or P4 were acquired without removing the sample from the stage, and the total scan time was about 6–7 days per sample.

As mentioned above, the internal wall diameter was determined from gray CT images. With this in mind, for 2–3 random gray slices of each high-resolution scan, we manually found a circle diameter, which fits best the internal tube diameter. The circle diameters were recorded with a 0.5 voxel accuracy and multiplied by the voxel size of a given scan. Then, the obtained tube diameters in micrometers were averaged over four high-resolution scans. The final average internal tube diameter was $8984 ± 1.6$ μm for P3 and $9044 ± 1.6$ μm for P4, where the indicated uncertainty is ±(one standard deviation). Using a similar procedure, the length of P3 was estimated as 59 649 ± 8 and 60 330 ± 9 μm for P4. However, here, the indicated uncertainty is ±(one average voxel size) which actually means ±(one slice No.).

After recording 2D projections and calculating a 3D reconstruction, we obtained a set of 16-bit .tif images with the gray values from 0 to 65 535, completely describing the internal 3D pore space of each sample at a given resolution. The next step was to mark wall voxels in order to exclude them from segmentation and segment only two phases: glass and air. For this purpose, we employed the circle-finding MATLAB function imfindcircles. Initially, only ±5%-approximate circle diameter was used with imfindcircles, and the function provided a pair of (X, Y) coordinates for the each .tif image. The black dots in Figs. 21(g) and 21(h) depict the whole set of (X, Y) coordinates after imfindcircles for one of the P3 high-resolution scans. Then, each set of the X or Y coordinates was fit with a separate polynomial [the red lines in Figs. 21g and 21(h)], and these two polynomials together with accurate tube diameter were used to initialize wall voxels, Fig. 21(f). Figure 21(i) shows four zoomed locations for P3 and P4, demonstrating the quality of wall placement for the described procedure and available glass tubes.

Black dots with the corresponding polynomial fits in Figs. 21(g) and 21(h) allow to assess several aspects: quality of the scan segment merging, straightness of the employed glass tube, and the sample alignment along the vertical rotation axis. First, there are no visible discontinuities between individual scan segments. Second, the axial deviation of X and Y circle coordinates is excellent, about 1.5 pixels for 7000 slices or 15 μm for 60 000 μm length. Third, the same deviation already includes the quality of the sample alignment along the rotation axis, which by itself took several hours before start of each sample scan. The accurate sample alignment along the rotation axis was necessary to employ later the periodic boundary conditions in flow simulations.

After marking the wall voxels, the rest of the 3D image was segmented using a global segmentation threshold. The gray voxels below (darker than) threshold were marked as fluid, while above—as solid [Figs. 23(a) and 23(b)]. The threshold value was chosen to match the voxel ratio fluid/(fluid+solid) to the porosity value determined from the mass of glass beads, their density, and volume of the confining container (⁠$εP3=35.414¯$%, $εP4=35.689¯$%). The confining container (glass tube) was assumed to be a perfect cylinder. We intentionally kept this step as simple as possible to avoid operator-dependence of the obtained results.

Figures 22(a)–22(c) and 22(o) show the examples of suboptimal and optimal choice of certain scan parameters. It appears that the loss of contrast in the image can be caused by various reasons, not just by the suboptimal x-ray source settings. However, the final result is always the same: reduction of contrast results in the formation of “bridges” between beads in the segmented image. The formation of bridges has simple origin. Because of the fixed target porosity of the segmented image, the numbers of fluid and solid voxels are fixed, too. With the decrease in contrast, the darker air voxels near bead contacts become brighter, and their larger fraction is marked as solid. However, because of the fixed total number of solid voxels, their number must be lowered somewhere else. First candidates to be marked as fluid are the voxels at bead boundaries, but away from contacts. This also means that the formation of bridges at lower permeability locations is accompanied by the opening of the rest of the pores at higher permeability locations. This leads later to an increase in simulated permeability with loss of image contrast.

While looking for the limits in image contrast, we implemented algorithmic methods to downsample original 3D images from their high-resolution versions, $r≈60$⁠, to medium and low resolutions, $3<r<25$⁠. The goal behind implementing downsampling was to preserve the spatial distribution of gray level while using a lower number of voxels. (The level of gray in each voxel of CT image is determined by the density and atomic number of the scanned material.) We employed two simple downsampling methods referred to as “via center” and “via volume.” The downsampling via center method merged f³ voxels of high-resolution image into one new voxel of the resulting lower resolution image. The factor f was chosen as an odd integer number $3,5,…,19,21$⁠. The gray value for each downsampled voxel was equal to the gray value of the central voxel of the cube formed by f³ voxels of the high-resolution image [Fig. 23(e)]. In other words, information from a single center high-resolution voxel was used, while the rest of $f3−1$ voxels were discarded. In downsampling via volume method, the gray value of each downsampled voxel was the average value of f³ high-resolution voxels [Fig. 23(f)].

Figure 23(f) shows that downsampling via volume reduces image noise because this downsampling acts, in principle, as a mean 3D filter. However, volume downsampling also blurs contacts between beads, which introduces additional bridges in the segmented image. Conversely, the downsampling via center provides noisy but higher-contrast image, also after segmentation. Figures 23(d)–23(g) show the identical regions of a sample from high and medium-resolution scans and the downsampling procedures. Obviously, the highest contrast and less bridging are observed in panel (d), and the contrast decreases toward (g) together with increasing bridging. Figures 23(a) and 23(b) provide the corresponding normalized gray-level histograms for the full 3D images, where the contrast loss manifests itself as an increasing central plateau. The high-resolution and downsampled via center histograms are identical. The quantitative impact of contrast loss and contact bridging on permeability is discussed below.

Acquired CT images at high resolutions are subject to so-called ring artifacts, Fig. 23(c), top left, which are mostly pronounced in high-resolution scans. As mentioned above, we did not apply image filters to remove ring artifacts because they can introduce other image distortions. Oscillations in gray levels due to ring artifact or minimal image noise can lead to solid voxels “hanging” in the middle of a pore, away from the solid boundaries [Fig. 23(c), bottom left]. These voxels break the parabolic-like flow profile inside a pore and slowdown flow in general. In order to remove such voxels, we applied a 3D median filter with (1,1,1) radius using the Fiji software,⁶⁷ Fig. 23(c). However, similar to downsampling via volume, 3D median filtering introduced bridging. We found an efficient solution to recover or “revert” some fluid voxels, which were turned into solid by 3D median filtering. Namely, we compared the original and 3D median filtered segmented images and identified clusters of new solid voxels introduced by 3D median filtering. The clusters were found using the bwconncomp MATLAB function with 26 nearest neighbor connectivity in three dimensions. If the number of solid voxels inside a given cluster exceeded 1, this cluster was reverted back to fluid. The locations of such reverted fluid voxels were recorded into a binary “revert mask,” Fig. 23(c), black. Then, this mask was applied to the original gray image, and the corresponding gray voxels were changed to very dark gray value to ensure that they would be fluid after segmentation. Finally, the 3D median filtered and reverted gray image was segmented again using a porosity-based global segmentation threshold. The result is shown in Fig. 23(c), bottom right, where the hanging voxel was removed while bridging was avoided. The implemented procedure also removed some minimal roughness at the caps of beads [the green dashed circles in Fig. 23(c)]. The 3D median filtering and revert mask were applied only to the four highest resolutions, r > 50. Some statistics for the numbers of fluid and solid voxels before and after filtering and reverting is available in supporting files. Namely, the fraction of affected solid voxels was $<0.05%$⁠, while the difference in permeability before and after this procedure was about 0.2%, i.e., mostly cosmetic. However, the presence of occasional hanging voxels was rather annoying, and the 3D median filtering with reverting was our way to do something about it. The binary revert mask as well as the original, 3D median filtered, and masked images in gray and segmented versions are available as supplementary material.

Comparing dimensions of the object scanned either with CT or with calibrated optical microscope (which is more precise equipment and has ∼12 times smaller pixel size compared to CT voxel size), we found the difference between the reported and actual voxel sizes of CT scans. The estimated difference was about + 0.25%, and, therefore, we applied the correction factor of 1.0025 to all voxel sizes. For this correction, we added a subjective uncertainty estimation of ±0.05%, assuming that the correction factor is between 1.002 and 1.003. The estimated porosities reported above then will have the following uncertainties: $εP3=35.414¯ ± 0.1$%, $εP4=35.689¯ ± 0.1%$⁠, or $εP3=35.324$ − 35.518%, $εP4=35.600$–35.792%.

To simulate Stokes flow of an incompressible fluid on the segmented 3D CT images, we employed the lattice Boltzmann method (LBM). LBM simulates evolution of fictitious particles, which collide at lattice nodes and propagate along lattice links at discrete time steps. We used a two-relaxation-times (TRT) version of the collision operator,⁶⁸ which is a truncated version of the multiple-relaxation-times collision,⁶⁹ and which is different in formulation from the popular Bhatnagar–Gross–Krook (BGK) collision operator.⁷⁰ However, all the obtained results can be reproduced with the BGK LBM, although at higher computational costs. A detailed description of the TRT LBM implementation can be found in our previous study,³⁶ together with a comparison with the established semi-analytical results for Stokes flow in complex geometries.^71,72

The no-slip boundary condition at solid–fluid interface was implemented using the “bounce-back” rule. For all simulations, we used the body force value of $10−9$ and the viscosity value of $νLBM=0.05$⁠, both in lattice units. The particular choice of $νLBM=0.05$ was motivated by a faster solution convergence, see, e.g., Fig. 7 in Ref. 37 and Fig. 6 in Ref. 73. The value of the so-called dimensionless “magic parameter” Λ (controlling truncation errors⁷⁴) was taken as 0.25. Such a choice of Λ makes the converged TRT permeability values equal within machine accuracy to the BGK values obtained with a commonly accepted relaxation rate of τ = 1. To demonstrate this point, we also ran a few simulations with 1/6 of the lattice viscosity, which together with $Λ=0.25$ is equivalent to BGK with τ = 1. Please note that the simultaneous choice of $νLBM=0.05$ and $Λ=0.25$ is not possible with τ parameter in the BGK collision operator.

At boundaries of the computational domain, we used periodic boundary conditions, including inlet and outlet, and applied a body force to drive the fluid. The periodic boundary conditions used together with the body force are equivalent to the pressure boundary conditions if the simulated geometry is periodic (e.g., see Table I in Talon et al.⁷³). CT scans of real porous samples are not periodic, and below we describe two approaches, with and without tubing, to deal with this case.

Without tubing, the simulations were performed on 80% of the length of each sample, Fig. 24, middle and bottom. This was motivated by the porosity distribution, where the porosity oscillations near both sample ends were observed. We also tested the 90%, 70%, and 60% length fractions, and, as expected, the particular choice had no impact on the following simulation steps. After extracting 80% length from each sample, the geometry was mirrored in the longitudinal direction, and flow was simulated. Because of the proper sample alignment with the axis of rotation during CT scans [Figs. 21(g) and 21(h)], the body force was applied only along the longitudinal direction, accurately mimicking experimental conditions.

The simulations with tubing were performed using hollow cylinders directly attached to 100% length of the P3 or P4 samples, including meshes. To mimic flexible tubing and plug hole, we used two hollow cylinders with diameters of 1.73 and 1.89 mm, respectively. These cylinders were discretized on a uniform cubic mesh using resolution equal to the resolution of each sample. Similarly to the case without tubing, proper sample alignment with the axis of rotation during CT scans [Figs. 21(g) and 21(h)] allowed the resulting geometry to accurately mimic the experimental geometry.

The lattice Boltzmann method updates the lattice links (and the resulting permeability) iteratively. After a sufficient number of iterations, the flow field for a given geometry does not change anymore and can be considered as converged. The iterative evolution of the calculated permeability is shown in Figs. 25(a) and 25(b), together with the relative change in permeability between two consequent iterations, Figs. 25(c) and 25(d). According to Figs. 25(a) and 25(b), the simulations on 80% of the sample length require about 10³ iterations to converge, while those on the full system comprised of the sample and tubes—about 10⁶. The difference originates from the heterogeneities of these two geometries: the 80%-sample-length geometry consists of similar pores between beads, while the full system additionally contains the long highly permeable tube. Because of the relative homogeneity of the 80%-sample-length, less strict convergence criterion can be used, such as $10−7$⁠, for the relative change in permeability between two iterations. (We used about $10−9$⁠). For the full system, we used a stricter criterion of $10−11$ to terminate simulations. According to Fig. 25(d), at the relative change of $∼10−11$⁠, the solution converges exponentially, which supports the assumption that the obtained permeability values for the full system are converged.

In Fig. 25, we demonstrate the solution convergence for both (i) TRT with $Λ=0.25$ and viscosity of 0.05, and (ii) BGK with τ = 1 and the corresponding viscosity of 1/6. Clearly, the TRT collision operator provides converged permeability values identical to BGK [insets in Figs. 25(a) and 25(b)] but with fewer iterations [Figs. 25(c) and 25(d)]. We also mention that taking the viscosity value of 1/6 and $Λ=0.25$ for TRT collision will produce results identical to BGK after each iteration.

After establishing criteria for the solution convergence, we use two methods to estimate accurate permeability values. The first method is based on our previous work,³⁶ and uses the linear extrapolation of permeability value for $r→∞$⁠. The second method is based on the similarity of the relative errors in permeability for different geometries. Both methods rely on a resolution study. Due to prohibitive computational cost [growing as $O(r5)$], we perform the resolution study only on 80% length of either P3 or P4 samples.

The first method suggests using $Λ=0.05$⁠, calculate permeability values at several higher resolutions (⁠$r≈60$ or higher) and estimate the “true” permeability value from linear extrapolation. The main idea is that the discrete truncation error always vanishes with resolution increase, but at a different convergence rate controlled by Λ. This error includes two contributions: the second-order convergence of the LBM away from solid boundaries and the first-order convergence of the boundary condition (Fig. 3 in Ref. 36). Both contributions can have opposite signs, which can lead to very slow transient error convergence, also slower than first-order. With $r→∞$⁠, the second-order contribution should decay faster leaving only the first-order contribution. However, for real complex geometries and $Λ=0.25$ (τ = 1 in BGK), the first-order-only convergence is never observed for any reasonable resolution range. Taking $Λ≈0.05$ increases the truncation error magnitude but also greatly accelerates transition to the first-order-only convergence with increasing r (Fig. 4 in Ref. 36). Hereafter, the true permeability value for $r→∞$ can be determined within ±0.1% using linear extrapolation (Figs. 5 and 6 in Ref. 36).

The second method is based on the observation that the relative error in permeability is similar for different geometries. This was observed at least for regular and irregular computer-generated sphere packings of similar porosity but different permeabilities (Fig. 9 in Ref. 37). Following the methodology described in Subsection 3.4.2 in Khirevich et al.,³⁷ we calculate the relative difference between the extrapolated and simulated permeabilities.

Below, we employ computer-generated packings, discretizing them at resolutions of the physical CT scans or higher. We find the relative difference between $Λ=0.05$-extrapolated permeability values (the first method) and permeabilities obtained with $Λ=0.25$ (the second method). Then, for CT-scanned geometries we correct their Λ = 0.25-permeabilities for the relative difference from the second method, and compare the result with the Λ = 0.05-extrapolated values.

We computer-generated two packings using about 12 000 diameters of beads1 and beads2, packing the virtual beads in the periodic containers with dimensions of $10×10×100$ bead diameters. Both generated packings have porosities of about 36.2% and are shown in Fig. 24, top. Hereafter, each generated packing is discretized on a uniform cubic mesh at various resolutions, and flow simulations are performed using LBM. The discretization resolutions are chosen to match the original CT scans of samples P3 and P4, or to exceed them.

The permeability values at higher resolutions, r > 50, are shown in Figs. 26(a) and 26(c). We performed two sets of simulations, with $Λ=0.25$ and $Λ=0.05$⁠, to apply both methods for accurate permeability estimation. As mentioned, $Λ=0.05$ results in larger error magnitude but close-to-linear error convergence behavior. The linear extrapolation for $Λ=0.05$ is shown in Figs. 26(b) and 26(d), for several sets of resolutions, $r=55,…,64,r=80,…,100, and r=112,…,140$⁠. The extrapolated values are close, with a slight improvement with the resolution increase, and the most accurate values are obtained for $r=112,…,140$⁠. The relative difference in extrapolated permeabilities between $r=80,…,100$ and $r=112,…,140$ for both geometries is about 0.1%, which is consistent with our previous results, Fig. 6 in Ref. 36. Therefore, we take the reference values for computer-generated packings as

and assume them to be ±0.1% accurate, as indicated with the gray rectangles in Fig. 26. Using the second method for permeability estimation with $Λ=0.25$⁠, we get the relative difference of about $−2.8,…,−2.6$% between the simulated permeabilities at $r=55,…,64$ and the corresponding reference value, both for P3 and P4 samples.

Before proceeding to CT-scanned images, we need to comment that due to discreteness in positioning the first and the last slice, resolution variation on real images introduces additional scatter in porosity values. For example, extracting 80% from the image with r = 55 and r = 64 will result in the slightly different porosity values, despite that both 100%-long geometries were segmented with the same porosity threshold for r = 55 and for r = 64. This is also noticeable in the porosity oscillations for the computer-generated and CT-scanned geometries, where the former overlap at two resolutions (gray solid vs black dashed in Fig. 24, middle and bottom), while the latter are slightly shifted (solid blue and solid green in Fig. 24, middle and bottom). To compensate for this porosity difference, we used $ε3/(1−ε)2$ correction present in, for example, the Kozeny–Carman equation. That is each permeability value is corrected by $[ε3/(1−ε)2]/[εref3/(1−εref)2]$⁠, where $εref$ is the porosity value for the 80%-length averaged over four highest resolutions. Such a correction is small (∼0.1% for most cases) but may have a visible impact on the extrapolated values.

Next, we apply two accurate permeability estimation methods to CT images scanned at high resolutions, r > 50. The extrapolation results for $Λ=0.05$ are shown in Figs. 26(f) and 26(h), with and without porosity correction. Hereafter, we apply the second method and adjust the $Λ=0.25$ permeability values for $−2.77,…,−2.59$% for P3, and $−2.77,…,−2.60$% for the P4 samples, averaging the result. The second method provides slightly higher estimations for the reference permeability, which is totally consistent with the computer-generated case, subject to porosity oscillations [blue dashed vs red solid in (e) and (g) and in (a) and (c) of Fig. 26]. Therefore, we used the permeability values from the second method as the reference, resulting in

These reference values are for the 80% length case.

We note that two accurate permeability estimation methods can be performed for different Λ-values, but following the described logic will lead to the same result.

After establishing the reference ±0.1% accurate permeability values, we can construct the full permeability error–resolution curves, Figs. 27(a) and 27(b), obtained for $Λ=0.25$ and the 80% length case. Here, we include not only the original scans and computer-generated packings but also the geometries downsampled from the high-resolution scans. The rationale here is to address not only the resolution but also image contrast impact on the simulated permeability values. Figures 27(a) and 27(b) demonstrate that (i) the permeability error overlaps between the computer-generated and CT-scanned geometries above r = 45, (ii) with a decrease in resolution, the lower-contrast images exhibit higher permeability error, and (iii) the P3 and P4 samples and their computer-generated twins demonstrate identical behavior. From the reference permeability estimation methods and from Figs. 27(a) and 27(b), it is clear that at low resolutions, the error is positive, it crosses zero at $r≈20$⁠, and with the resolution increase, the error starts very slow convergence from below to the reference.

For completeness, we note that this whole study was originally designed to obtain permeability error–resolution curves [Figs. 27(a) and 27(b)] for real geometries and to analyze their difference from the computer-generated geometries, if any.

Figures 27(a) and 27(b) suggest that the most accurate permeability values are not for $r≈60$⁠, but for $r≈20$⁠, if we use the images downsampled via center and $Λ=0.25$⁠. The employed permeability error–resolution analysis is performed for 80% length of either P3 or P4. For comparison with the experiment, we need to use the 100% sample length and add hollow cylinders mimicking tubing and plugs. The cylinders also introduce the error in flow simulations, which is straightforward to estimate using the Hagen–Poiseuille equation. The inner diameters of the hollow cylinders are between 63 and 81 lattice nodes. The relative error in the LBM flow rate relative to the analytical value is about −0.5% for $Λ=0.25$⁠. Here, both the negative sign and the error magnitude are consistent with the error analysis for porous samples. In order to alleviate the negative error from the hollow cylinders, and as the final flourish, we intentionally chose the r-value resulting in a minor positive error of about + 0.3%. Such geometries were obtained by downsampling via center with the factor of three the high-resolution images indicated with the black “ii” for P3 and red “ii” for P4 samples, Fig. 27(b). Each P3 and P4 downsampled geometry was appended with the hollow cylinders, resulting in the S3 and S4 geometries, and used for comparison with the experiment in Fig. 4.

According to Fig. 27(b), the minor positive error in permeability could also be obtained from the original (not downsampled) scans at $r≈28$ and $Λ=0.25$⁠, or at different values of Λ and accordingly chosen resolutions. However, using the approach presented here, one will end up with the similar final permeability values. In summary, we have eliminated the arbitrariness of choice of the Λ (or τ in BGK) with the corresponding choice of the discretization resolution.

To give more insights into the three-dimensional model system used to obtain results in Fig. 5, we provide the two-dimensional velocity and pressure profiles obtained from the LBM simulations. Figures 28(a)–28(c) and 28(d)–28(f) illustrate the results for a cylindrical holder with and without sample installed. Please notice the pressure values from both sides of the tube–sample connection and compare them for the case of installed sample and empty holder.

In order to compare our model with the well-established results, we plot the flow streamlines for the case of $Dsmpl/Dtube=4$ with the empty holder. Previous simulation and experimental studies suggest formation of a corner vortex,^75–77 which is reproduced well by the LBM simulations.

Figure 29 presents experimental setups from previous studies used for comparison in Fig. 3. The red and green color drawings in panel “Ma2011” sketch mesh and tube inside the geometry according to a private communication with the authors (Paolo Macini).

The results of Wyllie and Gregory were taken as the Kozeny–Carman constant “k” from their Table II and estimating permeability according to their equation (1) with $S0=6$⁠. The results from Khirevich's Ph.D. thesis were taken from Table 1.1. For the results of Loudon, we used their Table 4 (“mean spherical diameter” and “present experiments” columns), translating their hydraulic conductance values into permeability with g = 9.81 $m/s2$ and viscosity of water $μ=1.31$ mPa $·$ s at 10 °C.

The permeability values of James and Mclaren were calculated after digitizing their friction factor (f_JM) and Reynolds number (⁠$ReJM$⁠) from their Fig. 1, considering only 0.011 cm and 0.022 cm-diameter beads. The corresponding porosities were taken as 37.2% and 37.7%. The product of f_JM and $ReJM$ was used to estimate the constant (⁠$≈180$⁠) in the Kozeny–Carman equation (see Fig. 3 caption).

For the Fand et al. permeability values, we used the values of porosity and the Kozeny–Carman constant “k” from their Table I, “simple media” rows, together with their equation (10). The minor correction to their data was not to use their corrected “d_h” bead diameters from Table I (2.098, 3.072, and 4.029 mm), but “nominal diameters” from the paragraph below Table I (2, 3, and 4 mm).

The permeability values from the Macini et al. study were taken from their Table I, “glass beads” rows, “$kgas$” and “$kliq$” values. The average bead diameters were calculated as the average between the minimum and maximum values in the first column of their Table I.

Table II summarizes individual contributions and total values of the uncertainties for the experimental and simulation parts of this work. Each total value was calculated as the square root of sums of squares of the individual contributions.

In addition to the uncertainties discussed above, the uncertainty in the sample area, $ΔA/A$⁠, was estimated using the square of the uncertainty in the voxel size correction factor (±0.05%). According to Table II, the major experimental uncertainty originates from the measured viscosity values.

The uncertainty in permeability estimation due to discretization resolution, $Δk(r)/k(r)$⁠, is taken as its upper bound ±0.3% due to the intentional choice of that error for the sample region of the simulation geometry.

The main source of uncertainty in the simulated values we attribute to the uncertainty in voxel size, which defines the internal glass tube diameter and the resulting porosity of the sample region. We calculated $k(Δε/ε)$ for the aforementioned porosity ranges, $εP3=35.324$–35.518% and $εP4=35.600$–35.792%, using the relation $k∼0.8ε3/(1−ε)2$⁠. The additional estimated factor of 0.8 is due to the presence of hollow cylinders [tube and plug regions in Fig. 1(d)], which reduce the impact of sample permeability variation on the full system permeability.