Insights Into Estimation of Sand Permeability: From Empirical Relations to Microstructure-based Methods

The selected soil samples represent three distinctive types of sandy soils. Sample 1, identified as fine sand, exhibits a uniform grain size distribution, sample 2 was characterized as fine sand with traces of lignite, while sample 3 was classified as medium sand, demonstrating a less uniform grain size distribution. The samples were first subjected to macroscopic analysis, measurements of grain size distribution, bulk density, and specific density and falling-head permeability test in oedometer; the latter performed in accordance with the procedure described in EN ISO 17892-11:2019-05. These preliminary laboratory measurements were treated as reference data for this study. All soils were air-dried and in loose state of compaction. The following table summarizes the measured properties of each sample.

The grain size distribution analysis was conducted according to polish standard PN-EN 933-1:2012. Its results are provided in Figure 1.

Figure 1:

In order to acquire microstructure data for analysis, samples were scanned in Bruker Skyscan 1172 system, which consisted of a conical X-ray source with a maximum power of 10 W and 11 MPx CCD image sensor. The first general scan of the whole sample was conducted to accurately measure the length of the soil specimen for permeameter test, utilizing 1 MPx sensor configuration and oversize scanning function. The second scan, utilizing 4 MPx sensor, was to provide microstructure geometry data for further analyses.

3D reconstructions were performed with the NRecon software, which uses Feldkamp algorithm (Feldkamp et al., 1984). Smoothing, misalignment compensation, beam hardening, and ring artifact corrections were further performed to enhance the quality of the reconstructed images. Acquisition and reconstruction parameters were as follows: –

Resolution: 7.8492 μm/px

When compared between testing of soil permeability in a falling-head test in oedometer and micro-CT-based imaging methods, a significant challenge arises due to the difference in size of the samples. The inner diameter of the oedometer fixture used in the measurement (50 mm) is significantly larger than the diameter of the sample reconstructed in microtomography (10 mm). The scale effect may lead to problematic interpretation of the resulting permeability values, as the volume of the scanned sample may not be representative in dimensions of the sample in oedometer, and thus, the comparison of both approaches is of questionable validity. To address this challenge and overcome scale discrepancy, a custom unsteady-state permeameter was developed, which is capable of measuring the conductivity of the exact same sample as the one reconstructed in terms of micro-CT.

The core component of this system is a permeameter fixture fabricated using a 3D printing technology. It features a sample holder 10 mm in diameter, as well as inlet and outlet filters. Corrugated wall is incorporated to minimize the influence of boundary flow during testing. Accurate water-level measurement is achieved with the help of custom-made laboratory setup and the use of the computer vision (OpenCV 4.6.0., n.d.).

During testing, multiple soil samples are measured simultaneously. The water level as a function of time is later used to fit a theoretical flow curve as in Eq. (3): (3) Ht=H0⋅exp−K⋅Aa⋅L⋅t−t0 H\left( t \right) = {H_0} \cdot exp\left( { - K \cdot {A \over {a \cdot L}} \cdot \left( {t - {t_0}} \right)} \right) where H(t) is a hydraulic height as a function of time [m], H₀ is the hydraulic height during the start of the test [m], A is a sample cross-sectional area [m²], a is an inlet pipe cross-sectional area [m²], L is the sample length [m], and t₀ is a parameter accounting for starting time of the test [s].

Hydraulic height at the start of the test may be adjusted using a regulated outlet overflow. To account for the resistance induced by the apparatus itself, resulting hydraulic conductivity is corrected by the conductivity of permeameter without the soil sample. Measuring setup and permeameter fixture schematic are presented in Figures 2 and 3.

Figure 2:

Figure 3:

In the study, despite their widely known limitations, results of several empirical methods were compared to bridge the gap between commonly used laboratory approaches and micro-CT-based methods. As put together by Vukovic & Soro (1992) and reported by Řiha et al. (2018), those equations may be presented in standardized forms (Eq. 4, 5): (4) K=gv⋅C⋅fφ⋅(de)m K = {g \over v} \cdot C \cdot f\left( \varphi \right) \cdot {({d_e})^m} (5) K=C′⋅fφ⋅(de)m K = C' \cdot f\left( \varphi \right) \cdot {({d_e})^m} where K is the hydraulic conductivity [m/s], g is the gravitational acceleration [m/s²]; C, C′, and exponent m are unitless coefficients, f(φ) is a porosity function, and d_e is an effective grain diameter. Values of C, f(φ), m and percentile of grain size used as an effective diameter differ between methods. The equations used, values of coefficients, and applicability conditions are presented in Table 2.

Samples were scanned in a FDM-printed fixture later used in laboratory permeability measurements. To mitigate the possibility of soil particles displacing due to high angular acceleration of micro-tomograph holder while scanning, the top cover accommodated a piston providing small compressive force to the sample. Figure 6 shows rendered volumetric data of samples as processed by Bruker CTvox software. Figure 7 shows the process of extracting volumes of interests (VOI) from the retrieved data and their resulting binarized form. VOIs of arbitrarily chosen dimensions of 400³, 600³, and 800³ voxels have been selected to systematically gauge each method's sensitivity to varying VOI sizes.

Figure 6:

Figure 7:

The conductivities of the specimens were concurrently measured under consistent environmental conditions. Prior to the final measurement, samples underwent immersion in water for a duration of 12 h. Subsequently, a series of test runs was conducted under elevated pressure difference to ensure the complete saturation of the soil within the fixture. Raw data acquired with the help of computer-vision-augmented measuring system were subjected to preprocessing to eliminate readings distorted by the pipe holders and erroneous readings. Outcomes of these tests as well as the optimal fits of the theoretical flow curves are depicted in Figure 8.

Figure 8:

The measured hydraulic conductivity values have been corrected to account for the internal flow resistance imposed by the apparatus (i.e. filters) using the following equation (Eq. 8): (8) Kcorr=LsampLsampKmeas−LapKap {K_{corr}} = {{{L_{samp}}} \over {{{{L_{samp}}} \over {{K_{meas}}}} - {{{L_{ap}}} \over {{K_{ap}}}}}} where L_samp is the length of the soil sample measured from the general CT scan [m], K_meas is the measured series conductivity of the sample and system [m/s], L_ap is the length of the fictional sample, reflecting flow resistance of the measurement setup with no sample attached [m], and K_ap is the measured conductivity of the system in the reference test (i.e. without any sample) [m/s]. The results of calculations are depicted in Table 3.

The corrected soil conductivity values obtained through the employed apparatus demonstrated fair coherence with the results of a conventional falling-head test in oedometer. However, they present larger variabilities among individual tests despite similar environmental conditions and constant water temperature. Additionally, observations unveiled that the highly decreasing trend in permeability of sample 2 across successive trials (brown trace in Fig. 8 b) represents the first test, while the green trace denotes the last measured run through the sample, which may indicate influence of swelling, compaction-induced changes in the microstructure of this sample, or filter clogging in the duration of the trial.

Moreover, it is also worth noting that the theoretical curves employed as best-fitting models to the experimental measurements of sample 2 exhibited the worst alignment with the outcomes of the tests. This inconsistency between measured results and theoretical predictions may suggest the need of deeper inquiry into the testing methodology for samples containing larger amounts of silt and organic compounds (i.e., lignite).

Hydraulic conductivity of samples was assessed with the use of empirical formulae described in Section 2.4. Effective diameters used for calculations were derived from interpolated values of grain size distribution curves obtained from results of granulometric analysis.

The application of empirically derived formulae to estimate soil conductivity yielded a highly varied range of results as anticipated. Among them, the Seelheim formula stood out by correctly showing that sample 3 had the highest conductivity, while sample 2 had the lowest. The USBR formula in comparison generally provided results that were closest to the actual experimentally acquired values (8.94E-5 m/s, 1.15E-5 m/s, and 8.53E-5 m/s for samples 1–3 accordingly). Beyer, Hazen, Chapuis, and Seelheim formulae have consistently provided overestimated values of conductivity for samples 1 and 3. Interestingly, those formulae (apart from Chapuis formula) have shown smaller error for finely grained sample 2.It is also important to note that there was a significant variation in the results for each sample when different methods were used. This highlights the need for utilizing several distinct formulations when estimating permeability with the use of empirical relations, as the outcomes may vary significantly based on the chosen approach. The resulting conductivity values were summarized in Table 4.

Figure 9 presents grain size distribution curves acquired through the simulated sifting method, contrasting them with the actual grain size distribution curves derived from laboratory granulometric analysis. Figure 10 shows relative differences between hydraulic conductivity values calculated with the use of data from simulated sifting and those from granulometric analysis.

Figure 9:

Figure 10:

The simulated sifting method demonstrated satisfactory outcomes across all three samples. For soils with finer particle compositions (samples 1 and 2), slight differences appeared along the lower part of the grain size distribution curves. This behavior may be attributed to the possible aggregation of finer particles into larger grains, likely arising from intricacies in micro-tomography scanning, reconstruction, and further binarization of the datasets. Despite these differences, the simulated curves provided a satisfactory alignment with the measured values, which was reflected by the high determination coefficients for these samples (R²≥0.945). The magnitude of error can be potentially reduced by fine-tuning the scanning resolution and parameters of the segmentation process for such samples. Notably, expanding the volume of interest did not significantly alter the quality of the results.

For sample 3, characterized by the broader distribution of grain sizes, the method yielded slightly worse results for the overall fitting; however, this reduction was more prominent for the larger sizes of grains. Additionally, expansion of the VOI led to an improvement in a curve alignment, suggesting that the accuracy is impeded by the high size of representative elementary volume or lack of larger grains in scanned sample (given that the physical dimension of VOI at 800³ voxels roughly corresponds to 6.3 mm in the sample size).

However, it is important to note that for hydraulic conductivity estimation with empirical formulae, greater significance lies in the precise capture of the lower part of the grain size distribution curves, which notably yields best results for sample 3.

Semi-theoretical permeability estimation in this study was performed using the Kozeny–Carman equation (Eq. 7) and methodology detailed in Section 2.5.2. In contrast to simplified approaches commonly presented (Carrier, 2003; Kaviany, 1995), both tortuosity and specific surface area (SSA) were numerically computed using 3D-image data. The selected number of random walkers in the tortuosity simulation proved to be adequate, as the outcome remained fairly constant with its further increase. An illustrative representation of walkers' trails, acquired with the help of ParaView software (Ahrens et al., 2005), is provided in Figure 11.

Figure 11:

In the case of samples 1 and 2, changes in tortuosity and SSA were not notably influenced by the expansion of the VOI extent, leading to relatively stable estimations of conductivity for these samples. However, for sample 3, there was a consistent increase of conductivity as a result of altering the VOI. The summarized outcomes of the described procedure were presented in Table 5.

Computations were conducted according to the procedure described in the documentation of the OpenPNM framework (Gostick et al., 2016). Figure 12 (a)–(c) presents rendered pore-network models of samples 1–3 at 400³ VOI sizes. Each sphere represents a pore-type element, while each tube represents a throat-type element. Color denotes the physical diameter of each element. The flow was simulated along the Z-axis of the samples presented here, which was coherent with the flow direction in permeameter. Results of simulations are shown in Table 6.

Figure 12:

Extracted pore networks reflected the observed characteristics of each specimen's pore space. However, their topological simplification may be a major factor contributing to the observed overestimation of permeability values, especially for the second sample. This possible reason is supported by comparison of the tortuosity (1.479), specific surface area (62941 1/m), and porosity (0.584) of the network in question with values obtained from image analysis, which reveals substantial discrepancies.

In this study, lattice-Boltzmann simulations were conducted using Taichi-lang (Hu et al., 2019) based implementation (Yang et al., 2022). To facilitate high computational demands, a different machine had to be used as LBM is demanding in terms of both CPU/GPU and RAM. Despite its well-recognized suitability for massively parallel computations, the practical applicability of that advantage was hindered by the substantial dimensions of the selected VOI and resolution of the scans. These factors prevented the full harnessing of the parallelism, as the sizes of the samples exceeded what even high-end highly VRAM-equipped GPUs could accommodate – such as the Nvidia RTX3090 with 24GB of VRAM, where the practically feasible domain size was found to be constrained to less than 500³ voxels. Consequently, for the purpose of benchmarking, all simulations were conducted with CPU-based implementation for accurate comparison and evaluation. Even though analyses were made only for 400³ and 600³ VOIs, larger sizes would have made the computation impractically long.

Furthermore, it was also observed that the rate of convergence differed not only based on the domain size but also between samples themselves. Sample 2 exhibited the fastest convergence of flow rate toward its asymptotic limit, while sample 1 displayed the slowest. Thus, for accurate results, the flow rate approach to its limit should be closely monitored. Finally, it is worth noting that for samples 1 and 3, the computed permeability values decreased with the increase of the VOI extent, contrary to the behavior of sample 2. Table 7 presents the permeability and conductivity values, as derived from the LBM flow simulation. Figure 13 (a)–(c) depicts streamlines of the flow through the samples.

Figure 13:

Computational efficiency is a matter of great importance in investigation of soil properties from reconstructed images, especially considering large datasets to process. The minimal and maximal times of computation observed for different methods have been provided in Table S1 in supplementary materials, as well as the overview of the computational platform specifications.

The computation time varied among individual samples, with sample 3 generally requiring the shortest calculation time, while sample 2 the longest. That can be attributed to the fact that these specimens represent opposite ends of the porosity and feature number spectrum. The pore-network-based simulation exhibited the shortest computation time for all assessed VOI sizes, but it also showed the lowest overall accuracy. Interestingly, the computationally expensive lattice-Boltzmann method did not yield the most accurate results for the analyzed samples.