Effective hydraulic properties of 3D virtual stony soils identified by inverse modeling

Stony soils that have a considerable amount of rock fragments (RFs) are widespread around the world. However, experiments to determine the effective soil hydraulic properties (SHPs) of stony soils, i.e., the water retention curve (WRC) and hydraulic conductivity curve (HCC), are challenging. Installation of measurement devices and sensors in these soils is difficult, and the data are less reliable because of their high local heterogeneity. Therefore, effective properties of stony soils especially under unsaturated hydraulic conditions are still not well understood. An alternative approach to evaluate the SHPs of these systems with internal structural heterogeneity is numerical simulation. We used the Hydrus 2D/3D software to create virtual stony soils in 3D and simulate water flow for different volumetric fractions of RFs, f. Stony soils with different values of f from 11 % to 37 % were created by placing impermeable spheres as RFs in a sandy loam soil. Time series of local pressure heads at various depths, mean water contents, and fluxes across the upper boundary were generated in a virtual evaporation experiment. Additionally, a multistep unit-gradient simulation was applied to determine effective values of hydraulic conductivity near saturation up to pF=2. The generated data were evaluated by inverse modeling, assuming a homogeneous system, and the effective hydraulic properties were identified. The effective properties were compared with predictions from available scaling models of SHPs for different values of f. Our results showed that scaling the WRC of the background soil based on only the value of f gives acceptable results in the case of impermeable RFs. However, the reduction in conductivity could not be simply scaled by the value of f. Predictions were highly improved by applying the Novák, Maxwell, and GEM models to scale the HCC. The Maxwell model matched the numerically identified HCC best.

1 Introduction

Stony soils are soils with a considerable amount of rock fragments (RFs) and are widespread in mountainous and forested watersheds around the world (Ballabio et al., 2016; Novák and Hlaváčiková, 2019). RFs in soil are particles with an effective diameter larger than 2 mm (Tetegan et al., 2015; Zhang et al., 2016). Their existence in soil influences the two constitutive soil water relationships known as soil hydraulic properties (SHPs), i.e., the water retention curve (WRC) and the hydraulic conductivity curve (HCC) (Russo, 1988; Durner and Flühler, 2006). The accurate identification of SHPs is a prerequisite for the adequate prediction of water flow in soil with the Richards equation (Farthing and Ogden, 2017; Haghverdi et al., 2018). The SHPs depend on soil texture and structure (Kutilek, 2004; Lehmann et al., 2020) and are influenced by the presence of RFs in soil. It is generally accepted that RFs decrease the water storage capacity of soils as well as their effective unsaturated hydraulic conductivity. In contrast, the formation of macropores in the vicinity of embedded RFs may lead to an increase in saturated hydraulic conductivity. While experimental evidence and theoretical analyses show that the volumetric fraction of RFs, f ($v/v$), has the highest influence on the effective SHPs of stony soil, the effect of other characteristics of RFs such as their porosity, shape, size, arrangement, and orientation towards flow is less clear (Hlaváčiková and Novák, 2014; Hlaváčiková et al., 2016; Naseri et al., 2020). To date, two approaches have been dominant in identifying the hydraulic behavior of stony soils: (i) experimental setups with the aim of measuring the SHPs of stony soils in the field or in controlled systems in the laboratory (Cousin et al., 2003; Dann et al., 2009; Grath et al., 2015; Beckers et al., 2016; Naseri et al., 2019) and (ii) the development of empirical, physical, or physico-empirical approaches to scale hydraulic properties of background soil based on the value of f and characteristics of RFs (Novák et al., 2011; Naseri et al., 2020). These two approaches have some systematic limitations that restrict their application in investigating the hydraulic behavior of stony soils. Installation of sensors and measurement instruments in stony soils are technically demanding (Cousin et al., 2003; Verbist et al., 2013; Coppola et al., 2013; Stevenson et al., 2021), undisturbed sampling is laborious (Ponder and Alley, 1997), relatively larger samples are required (Germer and Braun, 2015), and the measured data might be more inconsistent due to the higher local heterogeneity of such soils (Baetens et al., 2009; Corwin and Lesch, 2005). Furthermore, some of the available scaling models to obtain effective SHPs are conceptually oversimplified, and they exclusively consider the value of f as the only input parameter (Bouwer and Rice, 1984; Ravina and Magier, 1984). Additionally, they assume impermeable RFs and are proposed mainly for saturated flow conditions. These scaling models require systematic verification under variably saturated conditions using experimental data or 3D simulations. Some reviews of these models and their evaluation are available in the literature (Brakensiek et al., 1986; Novák et al., 2011; Beckers et al., 2016; Naseri et al., 2019).

Hlaváčiková and Novák (2014) proposed a model to scale the HCC of the background soil, parameterized with the van Genuchten–Mualem (van Genuchten, 1980) model, using the model of Bouwer and Rice (1984). Hlaváčiková et al. (2018) used the water content of RFs as an input parameter to scale the WRC of the background soil. Naseri et al. (2019) used the simplified evaporation method (Peters et al., 2015) to experimentally determine the effective SHPs of small soil samples containing various amounts of RFs. Their study criticizes the application of the scaling models developed for saturated stony soils to unsaturated conditions and emphasizes the need to develop approaches that consider more RF characteristics to calculate the SHPs of stony soils.

Recent advancements in computational hydrology and computing power suggest the numerical simulation of soil water dynamics as a promising alternative to the measurement of the effective SHPs of heterogeneous soils (Durner et al., 2008; Lai and Ren, 2016; Radcliffe and Šimůnek, 2018). Numerical simulations have several advantages: they do not demand strict experimental setups; they are repeatable under a variety of initial and boundary conditions; and, in contrast to laboratory experiments, spatial and temporal scales are not restrictive factors in the simulations. These assets have made them a favorable tool in water and solute transport modeling in heterogeneous soils (Abbasi et al., 2003; Šimůnek et al., 2016). However, with few exceptions, heterogeneous soils like stony soils have been simulated only for simplified cases, i.e., either under fully saturated conditions or with reduced dimensionality, i.e., simulations of stony soils in two spatial dimensions (2D). Novák et al. (2011) calculated the effective saturated hydraulic conductivity (K_s) of soils containing impermeable RFs using steady-state simulations with the Hydrus-2D software, which solves the Richards equation in two spatial dimensions. They derived a linear relationship between the K_s of stony soils and f. Hlaváčiková et al. (2016) simulated different shapes and orientations of RFs in Hydrus-2D to obtain the effective K_s of the virtual stony soils. Beckers et al. (2016) used Hydrus-2D simulations to extend the investigations towards the impact of f, shape, and size of RFs on the HCC. They also identified the effective SHPs of a clay stony soil using laboratory evaporation experiments for f values up to 20 % ($v/v$).

The inverse modeling approach has been applied to identify effective hydraulic properties of soils in laboratory experiments (Ciollaro and Romano, 1995; Hopmans et al., 2002; Nasta et al., 2011), using lysimeters and field studies (Abbaspour et al., 1999, 2000), utilizing virtual lysimeters with internal textural heterogeneity (Durner et al., 2008; Schelle et al., 2013), and for the WRC of stony soils through field infiltration experiments (Baetens et al., 2009). Although theoretical studies and laboratory investigations on packed samples are insufficient to fully understand the hydraulic processes in stony soils, they do lead the way to the improvement and validation of effective models and their application at the field and even larger scales. Inverse modeling is arguably the best approach to achieve these aims because it allows one to validate effective models using process modeling. Our aim in this study was to investigate the application of inverse modeling to identify the effective SHPs of 3D virtual stony soils and to explore its applicability to these soil systems as an example of internal structural heterogeneity. We were interested in answering the following questions:

• i.

  Is it possible to describe the dynamics in the heterogeneous 3D system with the 1D Richards equation assuming a homogeneous soil?

• ii.

  If so, what are the effective SHPs of stony soils and how do they relate to the SHPs of the background soil?

To answer these questions we conducted forward simulations of water movement in 3D using the Richards equation as a variably saturated flow model. We created stony soils by embedding voids representing impermeable spherical RFs as inclusions into a homogeneous background soil. We then simulated transient evaporation experiments and stepwise, steady-state, unit-gradient infiltration experiments in 3D. The generated data were used as an input to a 1D inverse model to obtain the effective SHPs of stony soils, and these properties were used to evaluate and compare the available scaling models of SHPs for stony soils.

2 Materials and methods

2.1 Simulation model

The Hydrus 2D/3D software was used to generate virtual stony soils and simulate the water flow in the created 3D geometries. Water flow in Hydrus 2D/3D is modeled by the Richards equation (Šimůnek et al., 2006, 2008), which is the standard model for variably saturated water flow in porous media. The Hydrus 2D/3D software solves the mixed form of the Richards equation numerically using the finite-element method and an implicit scheme in time (Celia et al., 1990; Šimůnek et al., 2008, 2016; Radcliffe and Šimůnek, 2018). The 3D form of the Richards equation under isothermal conditions, without sinks/sources, and assuming an isotropic hydraulic conductivity is as follows:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}& ={\displaystyle \frac{\partial }{\partial x}}\left[K\left(h\right)\left({\displaystyle \frac{\partial h}{\partial x}}\right)\right]+{\displaystyle \frac{\partial }{\partial y}}\left[K\left(h\right)\left({\displaystyle \frac{\partial h}{\partial y}}\right)\right]\\ & +{\displaystyle \frac{\partial }{\partial z}}\left[K\left(h\right)\left({\displaystyle \frac{\partial h}{\partial z}}+\mathrm{1}\right)\right],\end{array}\end{array}$$

where θ is the volumetric water content (cm³ cm⁻³); t is time (s); h is the pressure head (cm); K(h) is the hydraulic conductivity function (cm d⁻¹); x and y (cm) are the horizontal Cartesian coordinates; and z (cm) is the vertical coordinate, which is positive upwards. We used the van Genuchten–Mualem model to parametrize the WRC and HCC (van Genuchten, 1980):

$$\begin{array}{}\text{(2)}& S_{\mathrm{e}}\left(h\right)={\displaystyle \frac{\mathit{\theta }\left(h\right)-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}={\left[\mathrm{1}+{\left(\mathit{\alpha }h\right)}^n\right]}^{-m}\end{array}$$

and

$$\begin{array}{}\text{(3)}& K\left(h\right)=K_{\mathrm{s}}{S}_{\mathrm{e}}^{\mathit{\tau }}{\left[\mathrm{1}-{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^m\right]}^{\mathrm{2}},\end{array}$$

where θ_s and θ_r are the saturated and residual water contents (cm³ cm⁻³), respectively; S_e(h) is the effective saturation (–); α (cm⁻¹) is a shape parameter; n is an empirical parameter related to the pore size distribution (–); $m=\mathrm{1}-\mathrm{1}/n$; K_s is the saturated hydraulic conductivity; and τ is a tortuosity/connectivity parameter (–).

2.2 The 3D geometries representing stony soils

The virtual stony soils in 3D were created by placing spherical inclusions in a background soil. In accordance with real laboratory experiments (not reported here), we generated virtual soil columns as cylinders with a diameter of 16 cm, a height of 10 cm, and a total volume of ≈2011 cm³. The inclusions were considered as voids representing impermeable RFs embedded in the background soil. Configurations and characteristics of the created 3D geometries of stony soils are illustrated in Fig. 1. Each spherical inclusion had a diameter of 3.04 cm and a volume of ≈ 14.7 cm³. Stony soils with different values of f were created by including different numbers of spherical inclusions in the soil column. A total number of 15, 27, 39, and 51 spherical inclusions in each column led to four volumetric RF contents of 11.0 %, 19.8 %, 28.5 %, and 37.3 % ($v/v$). Spheres were arranged in the column in three layers. The spheres' centers were at depths of 2.5, 5.0, and 7.5 in the column, and each layer was packed with one-third of the total number of intended spheres. Furthermore, observation points at selected nodes of the numerical grid were inserted at each of the three depths of the column (i.e., 2.5, 5.0, and 7.5 cm) in the background soil and not in close vicinity to the inclusions to provide time series of the soil water pressure head for the inverse simulations. For the background soil, a homogenous sandy loam soil was considered with the van Genuchten–Mualem model parameters θ_s=0.410 (cm³ cm⁻³), θ_r=0.065 (cm³ cm⁻³), α=0.01 (cm⁻¹), n=2.0 (–), τ=0.5 (–), and K_s=100 (cm d⁻¹). The targeted mesh size for the different simulations was set to 0.25 cm. The dependency of the numerical solution on the mesh size was tested with some refined meshes, and negligible differences in the results were obtained for different mesh sizes.

Figure 1 Visualization of the generated stony soils in 3D, including the dimension of the RFs and the soil cylinder as well as the location of the observation points (a). The bottom row shows the RF contents of 11.0 % (b), 19.8 % (c), 28.5 % (d), and 37.3 % (e), from left to right.

2.3 Forward simulations

We simulated evaporation (EVA) (Peters and Durner, 2008) and multistep unit-gradient (MSUG) experiments (Sarkar et al., 2019). For EVA, a linear distribution of the pressure head (−2.5 cm top, +7.5 cm bottom) was used as the initial condition. The boundary conditions were no-flux at the bottom, and they were atmospheric with a constant potential evaporation rate of 0.6 (cm d⁻¹) and zero precipitation at the top. The EVA experiments were simulated for 10 d, and the time series of pressure heads at each observation point, the initial volumetric water content, and the cumulative evaporation and evaporation rate were collected for later use in the inverse simulations.

In the virtual MSUG experiment, the soil column was initially fully saturated with a constant pressure head of 0 cm. A sequence of stepwise decreasing constant pressure heads was assigned to the upper and lower boundaries of the column. The duration of the virtual MSUG experiment was 100 d, and the pressure head in the upper and lower boundaries was simultaneously decreased stepwise to a pressure head of −100 cm. The applied pressure heads h_i were 0, −1, −3, −10, −20, −30, −60, and −100 cm, respectively. Time steps were chosen such that a steady-state flow condition was reached for each pressure step, indicated by identical water fluxes at the top (inflow) and bottom (outflow) boundaries and constant pressure heads at the observation points. The hydraulic conductivities at the respective pressure heads h_i were calculated by dividing the steady-state water flux rates (cm³ d⁻¹) by the total surficial area of the soil column (≈202 cm²).

The converging and diverging flow fields around obstacles produce spatially different pressure heads, even under unit-gradient conditions; as opposed to saturated conditions, these different pressure heads are associated with different water saturations and different local hydraulic conductivities under unsaturated conditions. We were interested in the extent to which this could lead to nonlinear effects in the derivation of the effective hydraulic properties, in particular the effective HCC. Furthermore, as the flow field for a given volume fraction of obstacles depends on dimensionality, i.e., is different in a 2D simulation compared with a 3D simulation, studying the effects in the unsaturated region was one of the main motivations for performing this numerical analysis in 3D.

2.4 Inverse modeling of evaporation in 1D

A 10 d EVA experiment in 1D was simulated with the Hydrus-1D software package (Šimůnek et al., 2006, 2008) to obtain the SHP parameters using inverse modeling. The generated data from the EVA and MSUG forward simulations in 3D were used as input to the 1D inverse simulations. Time series of the pressure heads at three observation depths, mean volumetric water contents in the column during the virtual EVA experiment, and the data points of the effective HCC from the virtual MSUG experiment were used in the objective function. The time series of the mean volumetric water content was calculated from the initial water content, cumulative evaporation, and soil volume. The measurement range for pressure heads used in the objective function was from saturation down to −2000 cm. This reflects a setup with laboratory tensiometers with boiling delay (Schindler et al., 2010). The time series of the simulated evaporation rates from the 3D simulations were used as the time-variable atmospheric boundary condition for the 1D inverse simulations. The 1D soil profile was 10 cm long and was discretized into 100 equally sized finite elements. Similar to the 3D simulations, three observation points were defined at depths of 2.5, 5.0, and 7.5 cm. A no-flux boundary condition was used at the bottom. The six parameters of the van Genuchten model occurring in Eqs. (2) and (3) were all simultaneously estimated by inverse modeling. The weighted least squares objective function was minimized by the Shuffled Complex Evolution algorithm (SCE-UA, Duan et al., 1992). The data obtained from the EVA experiment allow one to identify the WRC from saturation to the pressure where the tensiometers fail as well as the HCC in the middle to dry range of the SHPs (pressure heads of roughly between −100 and −2000 cm), whereas the MSUG provides a precise determination of the HCC in the wet range (Sarkar et al., 2019; Durner and Iden, 2011). The EVA experiments do not provide information on hydraulic conductivity near water saturation (Peters et al., 2015). Therefore, we included the data obtained from the virtual MSUG experiment in the objective function for the inverse simulation of the virtual EVA experiments in order to improve the uniqueness of the inverse solution and the precision of the identified HCC near saturation (see Schelle et al., 2010, for another example).

2.5 Predicting the SHPs of virtual stony soils by scaling models

The SHPs of stony soils obtained by inverse modeling were compared to SHPs that were predicted by available scaling models and used for their evaluation. Considering that f has the dominant influence on the WRC of a stony soil, a common approach is partitioning the WRC and HCC of stony soil based on the volume of each component in the soil–rock mixture and calculating the effective SHPs of stony soil using the volume averaging or the composite-porosity model. The general form of the WRC model considers the moisture contents of the background soil θ_soil (h) (cm³ cm⁻³) and embedded rock fragments θ_rock (h) (cm³ cm⁻³) to calculate the effective WRC of stony soils θ_m (h) (cm³ cm⁻³) (Flint and Childs, 1984; Peters and Klavetter, 1988) with the following form in the full moisture range (Naseri et al., 2019):

$$\begin{array}{}\text{(4)}& {\mathit{\theta }}_{\mathrm{m}}\left(h\right)=f{\mathit{\theta }}_{\mathrm{rock}}+\left(\mathrm{1}-f\right){\mathit{\theta }}_{\mathrm{soil}}.\end{array}$$

A typical assumption in stony soil hydrology is that the porosity of RFs is negligible. In this case, following Bouwer and Rice (1984), Eq. (4) reduces to

$$\begin{array}{}\text{(5)}& {\mathit{\theta }}_{\mathrm{m}}\left(h\right)=\left(\mathrm{1}-f\right){\mathit{\theta }}_{\mathrm{soil}}.\end{array}$$

For the effective hydraulic conductivity of stony soils, some scaling models are developed for saturated conditions that might apply to the hydraulic conductivity at any pressure heads. The simplest scaling model accounts only for the reduction in the cross-sectional area available for the flow of water. This leads to the following equation (Ravina and Magier, 1984):

$$\begin{array}{}\text{(6)}& K_{\mathrm{r}}=\mathrm{1}-f,\end{array}$$

where K_r (–) is the relative hydraulic conductivity of stony soil, i.e., $K_{\mathrm{r}}=K_{\mathrm{m}}/K_{\mathrm{soil}}$, where K_m is the effective hydraulic conductivity of the stony soil (cm d⁻¹), and K_soil is the conductivity of the background soil (cm d⁻¹).

In a more recent approach, Novák et al. (2011) developed a linear relationship based on the 2D numerical simulation results as a first approximation to scale the saturated hydraulic conductivity of stony soils:

$$\begin{array}{}\text{(7)}& K_{\mathrm{r}}=\mathrm{1}-\mathit{\alpha }f.\end{array}$$

The parameter α was reported to depend on the texture of the background soil, with a range from 1.1 for sandy clay to 1.32 for clay. This model is easy to apply, but it requires the estimation of the parameter α to calculate K_r. For our calculations, we assumed α=1.2 for the sandy loam background soil used in our study.

Another model that has been developed for mixtures with spherical inclusions is the Maxwell model (Maxwell, 1873; Corring and Churchill, 1961; Peck and Watson, 1979; Zimmermann and Bodvarsson, 1995). It takes the value of f, the hydraulic conductivity of the background soil, and the hydraulic conductivity of inclusions into account to calculate the hydraulic conductivity of stony soil. In the special case of impermeable inclusions, the Maxwell model reduces to

$$\begin{array}{}\text{(8)}& K_{\mathrm{r}}={\displaystyle \frac{\mathrm{2}(\mathrm{1}-f)}{\mathrm{2}+f}}.\end{array}$$

A recently developed model by Naseri et al. (2020), which is based on the general effective medium theory (GEM), allows one to consider the effects of the permeability, shape, and orientation of RFs on the effective HCC. For impermeable RFs, the GEM model reduces to the following form:

$$\begin{array}{}\text{(9)}& K_{\mathrm{r}}={\left(\mathrm{1}-{\displaystyle \frac{f}{f_{\mathrm{c}}}}\right)}^t,\end{array}$$

where f_c is the critical f with values between 0.84 and nearly 1, and t is a shape parameter with values between 1.26 and nearly 1.5 for spherical RFs. In this study, we set f_c=0.982 ($v/v$) according to the size ratio of the RFs to the background soil, and we set t=1.473 for spherical RFs (for details, see the appendix in Naseri et al., 2020).

It should be noted that all approaches apply at any pressure head h_i, i.e., the scaling that is originally developed for saturated conditions with locally constant hydraulic conductivity in the background soil is equally applied to unsaturated conditions.

3 Results and discussion

3.1 Flow field and the variability of state variables in the virtual MSUG experiment

Figure 2 visualizes the pressure head (cm), water content (cm³ cm⁻³), and velocity (cm d⁻¹) in a 2D cross section in the center of the soil column through the forward simulation of the MSUG experiment. The profile is shown for the steady-state flux situation with a pressure head of −100 cm and stony soil with f=28.5 %. The figure shows a considerable change in the flow velocities, even at the upper boundary. Moreover, as Fig. 2 illustrates, the conditions above an obstacle might be slightly wetter than below an obstacle, but the variations in the pressure head and the water content fields is very small. We note that this general finding was equally applicable for all other pressure head steps in the virtual MSUG experiment.

Figure 2 Visualization of the pressure head (cm), water content (cm³ cm⁻³), and velocity (cm d⁻¹) in a 2D profile in the center of the soil column during the forward simulation of the MSUG experiment.

3.2 Comparison of the relative K_s of the scaling models and the 3D simulations under saturated conditions

The dependency of the relative saturated hydraulic conductivity (K_r) on the percentage of RFs, calculated by different scaling models, and the obtained values from the first pressure in the virtual MSUG experiment are presented in Fig. 3. The results of the models are shown up to f=37.3 %, which was the highest value of f simulated in 3D. However, some of the evaluated models are theoretically valid for higher or lower values of f, e.g., 40 % for the Novák et al. (2011) model and higher values for the GEM model (Naseri et al., 2020).

Figure 3 Comparison of the values of K_r (–) from the virtual MSUG experiment in 3D (circles) as well as those calculated by different scaling models (solid lines) for f values up to 37.3 %. The dashed line connects the simulated data points of K_r shown by circles.

Obviously, the results of our simulations confirm a linear reduction in K_r with an increasing volume of RFs up to f=37.3 % ($v/v$) in the soil. The numerically obtained values of K_r are shown using circles and are connected by the dashed line in Fig. 3. The dashed line has a slope of −1.29, representing a higher reduction rate of K_r compared with the scaling of K_r that would be proportional to f, which is expressed by Eq. (6) and predicted by the model of Ravina and Magier (1984) (solid green line). This result supports the fact that, even in a stony soil with spherical impermeable RFs, the reduction in the hydraulic conductivity is higher than the reduction in the average cross-sectional area (which is statistically equivalent to the value of f). Hlaváčiková et al. (2016) found an even higher value of −1.45 for spherical RFs with a diameter of 10 cm. The model of Novák et al. (2011) performs better but also leads to a slight underprediction of the reduction in the effective saturated conductivity. The performance of this model could be improved by adjusting the parameter α to match data of the 3D simulation, but doing this would lead to an unfair comparison with the other models. The two models predicting a nonlinear relationship between K_r and f, GEM and Maxwell, show similar results at low RF contents of up to 10 %, with minor differences in the outputs of the models. Among all of the evaluated models for scaling K_s, the Maxwell model yields the closest match to the numerically identified values of K_r.

We note that these results may differ in natural soils, where an increase in the saturated hydraulic conductivity might be expected because of macropore flow in lacunar pores at the interface between the background soil and RFs (Beckers et al., 2016; Hlaváčiková et al., 2016; Arias et al., 2019). We have not included such a process in our 3D simulations.

3.3 Inverse modeling results for effective hydraulic properties

The observed and fitted time series of the pressure heads at the three representative observation points are shown in Fig. 4 for the simulated experiments of the four cases with different values of f. In each case, the fitted pressure heads at the three depths of the column (2.5, 5.0, and 7.5 cm) match well with the time series of the corresponding data from the 3D virtual experiments. Specifically, the match of the pressure heads at 7.5 and 5.0 cm is excellent, whereas there are slight systematic deviations at the uppermost level at the later stage of the virtual EVA experiment.

Figure 4 The time series of the 3D-simulated (circles) and 1D-fitted (solid lines) pressure heads at the observation points at three depths of the stony soil columns with the f values of (a) 11.0 %, (b) 19.8 %, (c) 28.5 %, and (d) 37.3 %. The observation depths are indicated using different colors: orange (upper, 2.5 cm), green (middle, 5.0 cm), and blue (lower, 7.5 cm from the top).

Table 1 shows the values of the root-mean-square error (RMSE) and mean absolute error (MAE) between the observed and fitted time series of the pressure heads at three observation points for values of f. According to Table 1, the fit is best for the lower f and in the middle of the column. The highest deviations occur for the highest f, but there is no clear trend. Overall, the RMSE and MAE values are in an acceptable range with respect to the observed values of pressure heads up to −2000 cm. This indicates that the time series of the pressure heads at multiple depths generated by the 3D simulations of the EVA experiments can be described successfully by the 1D Richards equation assuming a homogenous system with effective SHPs.

Table 1 The RMSE and MAE values between the observed and fitted pressure heads at three observation points for different values of f.

Figure 5 The WRCs and HCCs of the background soil (solid black line) and the identified effective WRCs (a) and HCCs (b) of the virtual stony soils (solid gray lines) with different values of f. The HCCs are also presented on a logarithmic scale. The dashed lines show the effective WRCs and HCCs calculated by the models of Bouwer and Rice (1984) and Ravina and Magier (1984). The circles on the HCCs present the data points of hydraulic conductivity obtained by the virtual MSUG experiment under near-saturated conditions up to pF≈2.

The identified SHPs are presented in Fig. 5. The solid lines in the figure show the WRCs and HCCs of the virtual stony soils obtained by inverse simulation (except the solid black lines, which are the WRCs and HCCs of the background soil); the dashed lines represent the WRCs and HCCs scaled using Eq. (5) and Eq. (6), respectively; and the circles on the HCC plots represent the discrete data points of hydraulic conductivity obtained by the virtual MSUG experiment. The WRCs and HCCs are presented on a pF scale, which is defined as pF=log 10(|h|), where h is the pressure head in centimeters (Schofield, 1935). The van Genuchten model parameters of the background soil and stony soils are shown in Table 2.

Table 2 The van Genuchten model parameters of the SHPs of the background soil and of the inversely determined effective SHPs of the virtual stony soils with different values of f.

According to Fig. 5 and Table 2, the value of the shape parameter α is independent of the value of f, and the change in the value of n is negligible (but might be systematic). The inversely identified WRC and the predictions from the Bouwer and Rice (1984) scaling model match almost perfectly for all RF contents. In agreement with this, there is also an excellent agreement of the values of the saturated (θ_s) and residual water contents (θ_r) (cm³ cm⁻³) between scaled and identified WRCs. The values of θ_s and θ_r in the WRCs are directly related to the values of f, and the respective values of the background soil and the WRC are scaled by this factor over the whole range of the soil water pressure head.

Similar to the WRC, an increase in f reduces the hydraulic conductivity over the whole pressure head range covered by the virtual experiments. However, in contrast to the WRC, the simple scaling model based on Eq. (6) cannot describe the reduction in HCC. Figure 5 shows that the model of Ravina and Magier (1984; dashed lines) underestimates the reduction in the effective HCC for all RF contents. The reason for this might be related to the local variations in the flow velocity in the soil column. It is shown in Fig. 2 that the variations in the water flow velocity might be considerable. The nonlinearities in the flow field and changes in the local conductivities, as well as an increased average flow path length, force a stronger overall conductivity reduction. Thus, the arrangement of RFs might affect the reduction in hydraulic conductivities, leading to different conductivities at the same value of f (Naseri et al., 2020). The degree depends on how the flow area is altered in the soil column due to the presence of RFs (Figs. 1, 2). This result is in agreement with Novák et al. (2011), who reported a higher reduction in conductivity compared with a reduction that is proportional to the RF content. Furthermore, it may also differ depending on the characteristics of RFs, such as their size, shape, and orientation towards flow (Novák et al., 2011).

We had to include the data points of hydraulic conductivity from the virtual MSUG experiment in the inverse objective function to get a precise identification of the HCC obtained by inverse modeling near saturation. The information content from the virtual EVA experiment gives a unique identification only when the flux rate in the system reaches the magnitude of the unsaturated hydraulic conductivity, which is around pF=1.5 to pF=2 for many soils (Peters and Durner, 2008). Although there are some discrepancies visible near saturation for the case with a high value of f, the resulting values of hydraulic conductivity from the virtual MSUG experiment and the inversely identified HCC using the virtual EVA experiment join well around pF=2 for all of the values of f. Therefore, the HCC could be described successfully from the saturation up to pF=3 using the inverse modeling of the virtual EVA experiment with added K support points from the virtual MSUG experiment. The overall results suggest that the effective hydraulic parameters of stony soils could be obtained by the corresponding real experiments, and the result is robust for both WRC and HCC, even if the uncertainty in the identified HCC is higher than that of the WRC (Singh et al., 2020, 2021).

3.4 Evaluation of the Novák, Maxwell, and GEM models using the identified HCC

As stated above, the model of Ravina and Magier (1984), which is a linear scaling approach of the hydraulic conductivity (Eq. 6), underestimates the reduction in conductivity in stony soils. We used the identified HCC as a benchmark to evaluate and compare more advanced models for scaling HCCs, namely the Novák, Maxwell, and GEM models (Eqs. 7, 8, and 9, respectively). Figure 6 illustrates the calculated HCCs of stony soils with different values of f using these models to scale the HCCs and then compares them to the HCC identified by the inverse modeling.

Figure 6 Evaluation of the Novák, Maxwell, and GEM models for scaling HCC of stony soils using the identified HCC as a benchmark. In each case, the HCCs were obtained for f = 11.0 %, 19.8 %, 28.5 %, and 37.3 % ($v/v$). The inverse identified curves are shown using solid lines, and the model results are shown using dashed lines. The value of f_c in the GEM model was set as 0.982 with the corresponding shape parameter t = 1.473, and the parameter α=1.2 was selected for the Novák model.

The HCCs calculated by the three models are generally in good agreement with the identified HCC in the observed range of pressure heads. All three models result in a more realistic estimate of the HCC compared with the simple linear scaling approach. While the Novák model slightly underestimates the identified HCC for all four RF contents, the GEM model, in contrast, overestimates the reduction in the hydraulic conductivity. The Maxwell model shows the same results as the GEM model except that it underestimates the HCC for the stony soil with f=37.3 %.

In order to compare the performance of the three models, the average deviation (d_avg) between the calculated and identified HCCs (logarithmic scale) was calculated to quantify the error of each model in the pF range from 0 to 3 (Table 3). The signs of numbers in Table 3 represent the tendency of the model to over- or underestimate the identified hydraulic conductivities. A negative number means that the model underestimates the reduction in hydraulic conductivity.

Table 3 Performance of the Novák, Maxwell, and GEM models quantified by the average deviation hydraulic conductivities (d_avg) on a logarithmic scale for different values of f.

Table 3 confirms the qualitative tendency toward an underestimation of the conductivity reduction by the Novák model and toward overestimation by the GEM and Maxwell models; however, Table 3 also shows that the difference between the three models is not large and is probably not of relevance in practice (the GEM model at the high value of f, which has the highest deviation, corresponds to a relative mismatch of hydraulic conductivities of about 6 %). Nevertheless, despite the potential of the three models with respect to predicting the HCC of stony soils, we think they require further evaluation using field-measured data of hydraulic conductivity under different experimental conditions.

4 Conclusions

In this study, we created virtual stony soils with different volumetric fractions of RFs in 3D and identified their effective SHPs from saturation up to pF=3 by inverse modeling of virtual EVA and MSUG experiments in 1D. We used the identified SHPs to investigate the performance of the available scaling models in stony soils, namely the Bouwer and Rice (1984) model for scaling WRCs, and the Ravina and Magier (1984), Maxwell (Peck and Watson, 1979), Novák (Novák et al., 2011), and GEM (Naseri et al., 2020) models for scaling HCCs of stony soils.

Our results show that the boundary fluxes and the internal system states in the virtual 3D EVA experiments, represented by the observed time series of pressure heads at multiple depths, could be matched well by 1D simulations, and the effective WRCs and HCCs of the considered stony soils were determined accurately. Comparison with the scaling models showed that, by assuming a homogeneous background soil and impermeable RFs, the effective WRC can be calculated from the WRC of the background soil using a simple correction factor equal to the volume fraction of background soil, (1−f). This is a result with practical implications for obtaining the WRC of stony soils. In addition, the scaling results for the HCC were promising. Our results confirmed that the reduction in K_r was stronger than that calculated by a simple proportionality to (1−f). The three models of Novák, Maxwell, and GEM consider that, and they performed adequately in predicting the effective HCCs of stony soils. The Maxwell model matched the numerical results best.

Care must be taken before generalizing these results to arbitrary conditions, e.g., highly dynamic boundary conditions with sequences of precipitation and higher and lower evaporation rates, which might yield different results due to the occurrence of nonequilibrium water dynamics and hysteresis. For real stony soils, changes in the pore size distribution of the background soil may result from the presence of RFs (Sekucia et al., 2020), with corresponding consequences for the effective SHPs. This influence was reported to be more common in compactible soils with a shrinkage–swelling potential (Fiès et al., 2002). In highly stony soils, where RFs are not embedded completely in the background soil, the existence of effective SHPs is still an open question. Finally, the impact of arrangement and size of RFs on evaporation dynamics and the effective SHPs needs to be understood. Tackling these problems requires a combination of experimental and modeling approaches.