Estimating scale dependence of saturated hydraulic conductivity in soils

Understanding the effect of scale on hydraulic and physical properties of soils has broad applications to scaling problems in hydrogeology, soil physics, and environmental engineering. The scale dependence of flow and transport is attributed to spatial heterogeneities, such as pore-size distribution and pore connectivity at small scales (e.g., core), fracture orientation and long-range correlations at large scales (e.g., field). In this study, we apply concepts from percolation theory to estimate the scale dependence of saturated hydraulic conductivity, 𝐾 𝑠𝑎𝑡 . For this purpose, we use a database including undisturbed soil samples from four Danish sites (Jyndevad, Tylstrup, Estrup, and Silstrup). The value of 𝐾 𝑠𝑎𝑡 was measured at small (100 cm 3 ) and large (6280 cm 3 ) scales. First, we apply a classification approach, widely used in petroleum engineering, to group soils based on their similarities in hydraulic properties using porosity and 𝐾 𝑠𝑎𝑡 measurements at the small scale. We detect nine different soil classes with the average flow zone indicator (FZI) from 0.05 𝜇𝑚 in class 1 to 9 𝜇𝑚 in class 9. Next, using percolation theory, we characterize the scale dependence of critical pore-throat radius. We use the critical path analysis to link the critical pore-throat radius to 𝐾 𝑠𝑎𝑡 and, consequently, determine the scale dependence of 𝐾 𝑠𝑎𝑡 . Comparing the theoretical estimations with the experimental measurements show that the percolation theoretic model reasonably well estimates the 𝐾 𝑠𝑎𝑡 at the large scale from the soil water retention curve and 𝐾 𝑠𝑎𝑡 measured at the small scale. We find the root mean square log-transformed error (RMSLE) values 0.45, 0.77, 1.9, and 2.05 (cm/day) for sites Jyndevad, Tylstrup, Silstrup, and Estrup, respectively. Results show that the theory tends to provide more accurate estimations in coarser textures and unstructured soils as well as soil classes with FZI values greater than 0.7 𝜇𝑚 .


Chapter 1 -Introduction
Understanding flow and transport in porous media has been an active subject of research due to broad applications in science and industry, particularly contaminated site remediation and wastewater disposal. Fluid flow and solute migration in field, such as groundwater discharge, transmissivity and dispersion, may be studied at different scales e.g., from a few meters to several kilometers. However, field measurements are typically time consuming and expensive.
Thus, scaling analyses and more specifically upscaling techniques are required to estimate one property at a larger scale from characteristics available at a smaller scale (Hopmans et al., 2002).
For this purpose, one needs to understand the effect of scale on flow and transport. For instance, the scale dependence of dispersion (De Smedt et al., 1986;Pachepsky and Hill, 2017;Gao et al., 2012;Younes et al., 2020), tortuosity (Ghanbarian et al., 2013;Ghanbarian, 2022a;Matyka et al., 2008), and hydraulic conductivity (Lauren et al., 1988;Guimerà et al., 1995;Schulze-Makuch et al., 1999;Hunt, 2006;Fallico et al., 2012;Ghanbarian, 2022b) have been investigated at various scales. However, accurate estimations of soil physical and hydraulic parameters at a larger scale from soil properties available at a smaller scale is still a great challenge.
The scale dependence of hydraulic conductivity has traditionally been investigated in subsurface hydrology at large scales using pumping tests (Garbesi et al., 1996;Javaux and Vanclooster, 2006;Sánchez-Vila et al., 1996;Schulze-Makuch, 1996;Schulze-Makuch et al., 1999). Such studies show that hydraulic conductivity increases with scale. For instance, Schulze-Makuch et al. (1999) proposed the following power-law relationship where Vs is the sample volume, Ksat is the saturated hydraulic conductivity, and a and b are empirical parameters. Although widely used by soil scientists and hydrologists to link to sample volume or scale (Zota et al., 2009;Fallico et al., 2010;Pachepsky et al., 2014;Brunetti et al., 2022), Eq. (1) is purely empirical (Fallico et al., 2012).
Numerical simulations have also been used to detect the representative elementary volume (REV) of a domain and to study the effect of scale (Aslannejad et al., 2017;Esmaeilpour, 2021;Mostaghimi et al., 2013;Sahimi et al., 1986). REV is the smallest domain size above which physical or hydraulic properties do not change with scale (Bear and Braester, 1972). Mostaghimi et al. (2013) investigated the REV for consolidated and unconsolidated porous media by numerically calculating specific surface area, , and porosity using micro-CT images of different sizes. They showed that the REV value for porosity and specific surface area was different than . Mostaghimi et al. (2013) demonstrated that the REV could be up to two times greater than the porosity and specific surface area REV. Recently, Esmaeilpour et al. (2021) investigated the effect of scale on both hydraulic and electrical conductivities using pore-scale simulations. They found that both hydraulic and electrical conductivities increased as domain size increased. Ghanbarian et al. (2021), however, demonstrated that depending on pore coordination number, hydraulic conductivity may increase or decrease with scale using pore network simulations. Ghanbarian et al. (2015) applied a machine learning method called contrast pattern aided regression and developed scale-dependent pedotransfer functions to estimate soil water retention curve (SWRC) and using soil samples from the UNSODA database. By including sample diameter and length as input variables, they showed that the accuracy of scale-dependent pedotransfer functions was higher than those developed without those that incorporate sample dimensions. Using a large database of nearly 20000 samples, Ghanbarian et al. (2017) demonstrated that scale-dependent pedotransfer functions developed by Ghanbarian et al. (2015) estimated Ksat substantially more accurately than several other models widely used in the literature.
Within the critical path analysis framework, depends on some critical pore-throat radius and electrical conductivity (Katz and Thompson, 1986;Ghanbarian, 2020). By applying the critical path analysis approach and assuming that Ksat is dominantly controlled by critical pore-throat radius (Hunt, 2006),  proposed Eq.
(2) to explain the scale dependence of hydraulic conductivity as follows: where ( ) is the scale-dependent saturated hydraulic conductivity, ( ) is the at the smallest scale, ( ) is critical pore-throat radius at the large scale, and ( ) is critical pore-throat radius at the smallest scale. Eq. (2) ignoring the effect of electrical conductivity on may be applied when electrical conductivity data are not available.
By comparing with pore-network simulations,  showed that Eq.
(2) estimated the scale dependence of permeability with high accuracy. They reported the relative error ranged from -3.7% to 0.74%.
Although the performance of Eq.
(2) was evaluated by  using pore-scale simulations, it has not yet been assessed with experimental measurements and agricultural soils, which are complex systems (Minasny et al., 2008). Therefore, the main objective of this study is to evaluate the accuracy of Eq.
(2) using a database of agricultural soil samples from four sites and different horizons.

Chapter 2 -Theory
Percolation theory addresses flow and transport in a network of pores through scaling relationships by taking the effect of interconnectivity into account (Sahimi, 2011;Hunt et al., 2014). Although initial percolation-based models were developed based on bond and site percolation classes and regular networks, more realistic and representative models were later proposed using irregular and disordered systems (Tyč and Halperin, 1989;Kogut and Straley, 1979).  generalized the methodology originally proposed by Hunt (2006) to link the scale dependence of to the scale dependence of critical pore-throat radius.
Using concepts of percolation theory, where is the critical volume fraction of pores, ( ) is the pore-throat radius distribution, and are respectively the minimum and maximum pore-throat radii, is the pore-throat length, and is the critical pore-throat radius.  related the critical pore-throat radius to the pore-throat radius distribution, ( ), typical pore-throat length, 0, and the system size, , as follows which implicitly explains the scale dependence of . As L increases, the ratio ( 0⁄( + 0)) decreases, if 0 remains constant. Accordingly, , as the lower limit of the integral, has to increase, if ( ) does not vary with the scale. in Eq. (4) is a universal scaling exponent whose value is 0.88 in three dimensions (Hunt et al., 2014). Based on Eq. (4), the larger the system size, the greater the , if ( ) and 0 remain scale-invariant. This agrees with the results of Koestel et al. (2020) who showed that as system size increased, the value of increased following a nonlinear trend.

-Experiments
The data used in this study are from a database of 196 soil samples from four sites (i.e., Jyndevad, Tylstrup, Estrup, and Silstrup) in Denmark published by Iversen et al. (2001). The samples, collected at two to three soil profiles by coring, are agricultural soils that range from sand to loam in texture according to Soil Survey Division Staff (1993). Table 3.1 summarizes the salient properties of the samples at each site. The large-scale samples were collected using a 6280 cm 3 (inner diameter 20 cm) hydraulic press core. The small-scale samples were collected using a 100 cm 3 (inner diameter 6.1 cm) core that was driven into the soil. All samples were protected from evaporation and physical disruption and stored at 2-5˚ C until experimental measurements commenced (Iversen et al., 2001).
At the small scale (100 cm 3 ), soil samples were back saturated and then placed on top of a sandbox to measure the SWRC at tension heads h = 10, 16, 50, and 100 cm H2O. The pressure plate method was then used at h = 160 and 1000 cm H2O. The SWRC at h = 15850 cm H2O was measured on disturbed samples after grinding and sieving through a 2-mm sieve. After that, the samples were re-saturated for at least 24 hours, and then the was measured using the constant head method and the Darcy equation (Klute and Dirksen,1986). At the large scale (6280 cm 3 ), was measured using the constant head method suggested by Klute and Dirksen (1986). Table 3.1 Some properties of the soil samples used in this study (Iversen et al., 2001).

-Classifying Soils
In this section, we present the theory of soils classification based on porosity and data measured at the small scale. For this purpose, we first invoke the approach proposed by Amaefule et al. (1993), widely applied in petroleum engineering to group rocks in sandstone (Bhattacharya et al., 2008) and carbonate (Martin et al., 1997) reservoirs, and then adopt it for soils. Amaefule et al. (1993) used the concept of mean hydraulic radius to modify the Kozeny-Carman model and generalized it as follows: where k is the permeability ( 2 ), is the porosity (cm 3 /cm 3 ), Fs is a shape factor (= 2 for a circular cylinder), is the tortuosity (> 1), and Sgv is the surface area per unit grain volume (1/ ). Amaefule et al. (1993) divided both sides of Eq. (5) by porosity and took the square root to have 1− √ Those authors then defined √ k as the quality index (hereafter SQI, soil quality index). The concept of SQI is similar to that of average linear velocity widely used in groundwater literature (Freeze and Cherry, 1979). Since normalized porosity , the pore volume-to-grain volume ratio (also known as void ratio), is , Eq. (6) can be rewritten as in which = 1 √ is called the flow zone indicator. Amaefule et al. (1993) proposed plotting the SQI against the to identify different types or classes. Based on their terminology, samples that lie on a straight line with a unit slope would correspond to the same hydraulic flow unit and group. Each hydraulic flow unit is characterized with a unique FZI value, which is the intercept of the unit-slope line with = 1.
To classify soils from four different sites, we determined soil permeability from dividing smallscale by fluidity factor. Then, the SQI value was calculated from soil porosity and permeability both measured at the small scale. The value of was also determined from porosity measurement at the small scale. By plotting the SQI against the , we draw lines with unit slopes using different FZI values and captured soil samples belonged to the same trend.

-Estimating Large-Scale Ksat
Since the soil water retention data were not available at the large scale, we assumed that the pore-throat radius distribution, ( ), does not change with the scale. To determine the scale dependence of via Eq. (4), we first calculated the value of at the smaller scale using the SWRC by setting equal to (ℎ = 15450 )/ and solving Eq. (3) numerically. For this purpose, the ( ) was determined from the SWRC measurements and by fitting the van Genuchten (1980) model to the data where and are the saturated and residual volumetric water contents, respectively, h is the tension head, and , n, and m are shape parameters. We fitted Eq. (8) directly to the measured SWRCs using the nonlinear least square approach in MATLAB and optimized the values of , , , and simultaneously. For this purpose, we considered and as two independent fitting parameters. For the sake of higher flexibility, we also let the parameter n to be less than 1 and m greater than 1.
After optimizing the van Genuchten (1980) shape parameters, the pore-throat radius distribution was calculated as (Dexter, 2004) 0.149 where all variables have previously been defined. The typical pore-throat length, 0, was not known in our samples. We, therefore, estimated its value for each sample via Eq. (4) from and ( ) determined at the small scale. We also assumed that 0 does not change with the scale, and the optimized value of 0 was used to determine the at the large scale via Eq. (4). Porethroat lengths and their distribution are not typically measured in soil surveys. Since such 1 ∑ information were not known for the soil samples studied here, we treated in Eq. (4) as a constant and canceled it out.
After calculating the values of at the small and large scales, we used Eq.
(2) to estimate the large scale from the small scale one. To evaluate the predictability of Eq.
(2), we calculated the root mean square log-transformed error as follows: where is the number of samples, est is the estimated , and meas is the measured value.

Chapter 4 -Results and Discussion
The average values of the fitted van Genuchten (1980) model parameters and their standard deviations are presented in Table 4.1. We found the average and n values for the Silstrup and Estrup sites smaller than the Jyndevad and Tylstrup sites. More specifically, the average value of n for the Estrup site was 0.55. Although n > 1 is typically reported in literature, our results are consistent with those estimated by Zakizadeh Abkenar and Rasoulzadeh (2019).
They found n to be 0.73 and 0.63 in a loamy and a clay loamy soil sample, respectively. As discussed by van Genuchten and Nielsen (1986), one may expect n values less than 1 when m and n are considered to be independent variables, particularly in structured and/or coarsetextured soils. In our study, we also found n < 1 for sandy loam and clay loam samples in the Estrup site (Table 3.1). Since empirical parameters n and m are independent fitting parameters in our study and the purpose of fitting the van Genuchten SWRC model to the data was deriving a continuous form of the ( ), we do not expect n values less than 1 impact the obtained results. The average values of were 0.121, 0.039, 0.172, and 2.06 × 10 −6 cm -1 for the Jyndevad, Tylstrup, Silstrup, and Estrup sites, respectively (Table 3.1). We found similar average values for the Jyndevad and Silstrup sites, although the two sites have different soil textures (sand versus sandy clay loam). Generally speaking, the samples from the Silstrup site are finer in texture compared to Jyndevad. We should note that is inversely proportional to the air entry tension head, a parameter that is affected by soil structure more than its texture, particularly in undisturbed samples (Ghanbarian-Alavijeh et al., 2010). Recall that the SWRC data were only available at the small scale. As can be seen, the volumetric water contents measured at 15540 cm H2O are greater in the Silstrup (Fig. 4.1c) and Estrup (Fig. 4.1d) samples compared to Jyndevad (Fig. 4.1a) and Tylstrup (Fig. 4.1b). This indicates the soils from Silstrup and Estrup are finer in texture compared to those from Jyndevad and Tylstrup. This is consistent with the soil textures reported from the different horizons in Table 3   which indicates the effect of scale (Bear and Braester, 1972). The histograms from one site are also different than those from another site, which is due to spatial heterogeneities, particularly differences in soil textural and structural characteristics. Recall that soil samples from the Jyndevad and Tylstrup sites are sandy, while those from Estrup and Silstrup are loamy and more structured (Iversen et al., 2001)

-Soil Classes
Soil classification results based on the Amaefule et al. (1993) approach are presented in Fig. 4.3a. We detected nine soil classes each of which is characterized by a different average FZI value. We found that the FZI values spanned over two orders of magnitude in variation from 0.05 to 9 (Table 4.2), which indicates substantially different soil hydraulic properties within the soil database studied here. As reported in Table 4.2, classes 4 and 9 are the smallest and largest groups with respectively 5 and 48 soil samples. We found that soils with higher FZIs (≥ 2.35 ) were more frequent in the database analyzed in this study. This can also be inferred from Fig. 4.3a showing that the majority of samples have SQI > 0.6 .
The FZI is a parameter that differentiates pore geometrical facies based on soil texture and mineralogy (Amaefule et al., 1993). Soil samples in classes with greater FZIs (e.g., classes 8 and 9) incorporate larger pore throats in general. Fig. 4.3b displays the distribution of soil classes in each site. We found that the Jyndevad and Tylstrup sites included samples from four and five classes, respectively, were the least diverse sites. In contrast, the Silstrup and Estrup sites were the most diverse ones contained soil samples from all nine classes detected. This is well in accord with the fact that the loamy soils from Silstrup and Estrup are more structured than those from Jyndevad and Tylstrup.  We should note that using the Amaefule et al. (1993) approach one should expect unstructured/coarse-textured and structured/fine-textured soils with similar permeability and porosity values to be grouped into the same soil class. However, even if their permeabilities are similar, an unstructured/coarse-textured soil typically has a lower porosity compared to a structured/fine-textured soil. This means their values might be different and so their soil classes. Table 4.3 summarizes values used to estimate the large-scale . The average values of minimum and maximum pore-throat radii (rtmin and rtmax) were determined respectively from the maximum and minimum tension heads using the Young-Laplace equation (assuming zero contact angle) and the measured SWRC. The rtc values were calculated at the small and large scales. As can be seen, the average rtc at the large scale is greater than that at the small scale (Table 4.3). Table 4.3 The average values of minimum and maximum pore-throat radii (rtmin and rtmax) determined from maximum and minimum tension heads and the measured SWRC, critical porethroat radius (rtc) at the small scale computed by setting equal to ( = 15450 )/ and solving Eq. (3) numerically, and critical pore-throat radii (rtc) at the large scale computed by solving Eq. (4). Numbers in parentheses represent the standard deviations.

Site
The estimated values at the large scale against the measured ones are presented in reasonably well for all samples except three (all from class 7) for which the was underestimated by more than one order of magnitude. In structured soil samples, such a great difference between theory and experiment may be attributed to the presence of macropores (Jarvis, 2008;Iversen et al., 2012) that are not captured by the SWRC. However, by means of computed tomography imaging one should be able to characterize soil structures and quantify macropores (Heijs et al., 1995;Elliot et al., 2010;Bölscher et al., 2021). Most samples from the Jyndevad site are unstructured and sandy with low organic contents (Iversen et al., 2001).
However, those three samples from class 7 are from the Ap horizon (depth 0-32 cm) with high organic matter (Table 3.1). This means they might be weakly structured with some macropores present.   we found that the average ratio of large-scale to small-scale was 171.6. This means that the large-scale was, on average, more than two orders of magnitude greater than the small-scale . Interestingly, the average ratio for other samples from the Tylstrup site was 1.12. Given that those three samples are from the Ap horizon (depth 0-32 cm) with high organic matter (  Although for some soil samples the proposed model overestimated the large-scale , results indicate that its value was mainly underestimated in this site. This is most probably because the soil samples from the Silstrup site are loamy and structured (Iversen et al., 2001), and, thus, one may expect macropores to be present, as discussed above. We found the RMSLE value of 1.9 cm/day, which is nearly four times greater than that in Jyndevad and three times greater than that in Tylstrup. Similar results were obtained for the Estrup site ( Fig. 4.7). The Silstrup and Estrup

-Limitations
The underestimation is probably because the lowest tension head on the SWRC was h = 10 cm H2O. This means some large pores corresponding to tension heads less than 10 cm H2O were not captured, although they effectively contributed to and its value. We also assumed that the pore-throat radius distribution and typical pore-throat length do not significantly vary with the scale. Such an assumption was required because unfortunately the SWRC measurements as well as pore-throat length data were not available at the large scale in this study. Results of Tinni et al. (2012), Chen et al. (2015), and Shu et al. (2020)  estimates the at the large scale from the and SWRC measured at the small scale in coarse-textured soils more reliably than that in fine-textured soils. Recall that the soils from Silstrup and Estrup are finer in texture and more structured compared to those from Jyndevad and Tylstrup (Table 3.1). As stated earlier, Eq.
(2) only includes the effect of critical pore-throat radius. However, theoretically, depends not only on the value of critical pore-throat radius, but also on electrical conductivity or formation factor (Katz and Thompson, 1986;Ghanbarian, 2020). Therefore, one should expect more accurate estimations, if soil water retention and electrical conductivity data both are available across scales.
Since Eq.
(2) does not take into account the effect of electrical conductivity,  argued that the relationship between and may not perfectly follow the quadratic relationship. They accordingly generalized such a relationship by replacing the power 2 in Eq.
(2) with an optimized exponent. For their simulated pore networks,  found an exponent of 3.03 instead of 2. By comparing with pore-scale simulations, they showed that Eq.
(2) with the generalized exponent yielded more accurate estimations than the exponent 2. The exponent 2 is a reasonable choice for the Jyndevad and Tylstrup sites with coarser soil textures. It is also close to 2.22, the experimental exponent reported by (Ghanbarian and Skaggs, 2022) for coarse-textured soil samples from the GRIZZLY database. One may, however, need an exponent greater than 2 to improve the large-scale estimations for the Silstrup and Estrup sites with finer soil textures.
In this study, we estimated the at the large scale (6280 cm 3 ) from the and SWRC measured at the small scale (100 cm 3 ). This was required because the Iversen et al.
(2001) dataset did not include porosity or SWRC at the larger scale. Since within the critical path analysis framework depends on both critical pore-throat radius and electrical conductivity, further investigations are required to evaluate the accuracy and reliability of Eqs.
(2)-(4) using a large database including samples of various textural properties with SWRC, electrical conductivity, and measured at difference scales.

Chapter 5 -Conclusion
Understanding the effect of scale on hydraulic and physical properties of soils has been challenging in soil science and hydrology. In this study, we used concepts from percolation theory and critical path analysis to estimate the at the large scale (6280 cm 3 ) from soil water retention curve and saturated hydraulic conductivity measured at the small scale (100 cm 3 ). To evaluate the proposed approach, we used 196 soil samples collected from different horizons and sites in Denmark. We first adapted a classification technique widely used in petroleum engineering and grouped soils using the porosity and permeability measurements available at the small scale. We detected nine soil classes with average flow zone indicator spanned from 0.05 in class 1 to 9 in class 9. Using the proposed model, we then estimated the at the large scale from the and SWRC measured at the small scale. By comparison with experiments, we found RMSLE = 0.45 (cm/day) for the Jyndevad site with 73 samples, 0.77 (cm/day) for the Tylstrup site with 25 samples, 1.90 (cm/day) for the Silstrup site with 45 samples, and 2.05 (cm/day) for the Estrup with 53 samples. We discussed that although the estimations were reasonable, for higher accuracy one needs to collect electrical conductivity and soil water retention data consistently at all scales. This is because based on critical path analysis not only depends on critical pore-throat radius but also electrical conductivity. Further investigations are required to evaluate the proposed model using a broader type of soils with various levels of heterogeneity.