Effect of Stoniness on the Hydraulic Properties of a Soil from an Evaporation Experiment Using the Wind and Inverse Estimation Methods

: Stony soils are distributed all over the world. The study of their characteristics has gained importance lately due to their increasing use as agricultural soils. The effect that rock fragments exert on the soil hydraulic properties is difﬁcult to measure in situ, and is usually derived from the ﬁne earth properties. However, the corrections used so far do not seem accurate for all types of stony soils. Our objective was to assess the adequacy of estimating the hydraulic properties of a stony soil from the ﬁne earth ones by correcting the latter by the volume occupied by rock fragments. To do that, we ﬁrst assessed the validity of different approaches for estimating the hydraulic properties of a stone-free and a stony (40% rock fragments) cylinder prepared with samples from the same silt loam soil. The functions relating to the soil hydraulic properties ( θ - h , K - h - θ ) were estimated by the Wind method and by inverse estimation, using data from an evaporation experiment where the soil water content and pressure head were measured at different soil depths over time. Results from the evaporation experiment were compared to those obtained by applying the equation that corrects ﬁne earth properties by the rock fragments volume. Wind and the Inverse Estimation methods were successfully applied to estimate soil water content and hydraulic conductivity from the stony soil experiment, except for some uncertainties caused by the limited range of suction in which the experiment was conducted. The application of an equation for adjusting the soil water content at different pressure heads (allowing for deﬁning the soil water retention curve, SWRC), and the unsaturated hydraulic conductivity ( K ) directly from the stone content was not satisfactory. K values obtained from the measured data were higher than those inferred by the correcting equation in the wet range, but decreased much faster with a decreasing pressure head. The use of this equation did therefore not take into account the effect that the creation of lacunar pores by the presence of rock fragments likely exerts on water ﬂow processes. The use of such correction needs therefore to be revised and new approaches are needed for estimating the hydraulic conductivity in stony soils. In relation to SWRC, a new equation to calculate the water content of a stony soil accounting for the inﬂuence of possible lacunar pores is proposed.


Introduction
Soils that contain rock fragments (particles > 2 mm in diameter) are usually denoted as stony soils. Over the last few years, the interest on studying the effect that rock fragments exert on soil physical properties has grown because soils containing these fragments exist worldwide. Only in the Mediterranean region, stony soils already occupy more than 60% of the land [1]. In Africa, the Sahara, Sahel, and most soils in the geological mantle in the west of the African continent are gravelly [2]. Areas with stoniness from 15% to 80% are predominant in Central America and many areas of Chile, Peru, Venezuela and Brazil [3]. In some regions, with an increasing demand for food, fiber and renewable energies, these soils are being converted from less intensive land use to agricultural crop production.
Rock fragments play an important role in soil water dynamics because they have an influence on the soil water content and hydraulic conductivity [4,5]. This influence is complex and depends on the volume fraction that they occupy in the bulk soil, their geological origin and grade of weathering, their size and also their position in the soil [4]. Intuitively, the presence of rock fragments reduces the cross-sectional area of water flow, leading to lower hydraulic conductivity [6][7][8][9][10]. However, their existence can create lacunar pores at the rock fragment-fine earth interface that could become paths for preferential water flow and thus increase the saturated hydraulic conductivity [11][12][13][14][15]. Additionally, their presence reduces the total porosity and the spatial pore network architecture depending on their volume, size, shape and the type of soil studied [16]. Some rock fragments are even able to hold water depending on their origin [4,17]. Under unsaturated conditions, a volume of water that enters a soil with rock fragments will cause a larger matric potential or soil water pressure head (h) increase than the same volume of water entering the same soil without stones. With the larger h-increase, a more rapid increase in hydraulic conductivity is expected for the soil with rock fragments than for the soil not containing them.
Our knowledge of the hydraulic properties of soils with rock fragments is, however, still limited [18]. One reason for this is that the common sampling methods for obtaining undisturbed cores used in non-stony soils are impracticable when large volumes of stones are present. Moreover, larger volumes of soil need to be sampled in order to obtain a representative elementary volume (REV) [9] in stony soils than in non-stony soils. Inserting probes, access tubes or installing lysimeters in-situ, which turns out to be complicated without altering the soil structure [4][5][6]8] in non-stony soils, becomes nearly impossible in soils containing stones.
Due to these challenges, some authors have tried to derive stony soils´hydraulic properties from those observed in the fine fraction of the soil [5,6,9,10] and corrected for the volume of rock fragments [7,19]. However, for certain types of soils, the existence of rock fragments can lead to complications and unrealistic derivations of hydraulic properties from the fine earth fraction [6,8,18]. Due to discrepancies in the results the question remains open whether a simple correction can be applied for different soil types, shapes and sizes of rock fragments and whether this correction would be valid over a wide range of soil water content (θ), from wet to relatively dry conditions [20].
Some studies on stony soils have mostly focused on measuring infiltration rates and hydraulic conductivity [12][13][14][15] or, on estimating the soil hydraulic properties by using soil water transport models [21][22][23]. The results found in the literature show a different influence of stoniness on soil hydraulic properties depending on the specific characteristics of the soil (e.g., texture, organic matter) and stones used (e.g., geological origin, shape, size) and also the volume which the latter occupy [4,8,9,24]. The discrepant results could be also attributed to the lack of a reference method for measuring hydraulic properties in stony soils and the different assumptions made by the authors during calculations. All in all, there is little information on the behavior of stony soils in terms not only of water retention, but also in relation to hydraulic conductivity. In particular, not many studies have looked at hydraulic conductivity and its relation to the pressure head allowing for the identification of a continuous function.
Repacked soil cores have also been used to study the hydraulic properties of stony soils [8,11,[25][26][27]. Although they could not reproduce soil structure, experiments based on repacked soils with a controlled volume of rock fragments helped to derive a basic understanding on the impact of different degrees of stoniness on the hydraulic properties, to test the accuracy of the estimation methods used to derive the soil hydraulic parameters, and also to evaluate the adequacy of estimating these soil properties by simply correcting fine earth hydraulic properties by the volume the rock fragments occupy.
The evaporation method is a reliable and simple method to determine soil hydraulic properties, i.e., the soil water retention curve (θ(h)) and the unsaturated hydraulic conductivity (K) as a function of h. It was developed by Reference [28], modified and validated by Reference [29] and used in various ways, e.g., by Reference [30].
Another approach to analyzing experimental data obtained in an evaporation experiment is the inverse parameter estimation technique. Inverse methods are based on a water transport simulation algorithm, such as the finite element solution of Richards' equation, HYDRUS [31]. For simulating water transport, soil hydraulic properties can be parameterized according to Reference [32]. The simulation is coupled with a parameter optimization algorithm. Knowing from the evaporation experiment the upper and lower flux boundary conditions, h at different depths and the total water storage S in the sample over time t, HYDRUS can be used to estimate soil hydraulic function parameters. If soils with and without rock fragments are measured and analyzed, the impact of the volume fraction of these rock fragments on the hydraulic function parameters can be derived.
In summary, despite its relevance for modeling and management decisions, the way in which the presence of stones affects unsaturated soil hydraulic parameters compared to the same soil without stones is still not well understood and quantified. Thus, for a given soil without rock fragments, it is not yet known how to adjust soil hydraulic property parameters if the same soil contains a defined volume of rock fragments. In this context, the objectives of this study were to: (1) estimate the water retention characteristics and hydraulic conductivity of a silt loam from evaporation experiments using the Wind method and the inverse parameter estimation, and (2) test the validity of corrections applied to stony soils for estimating the soil water content and the hydraulic conductivity from the fine earth properties for this soil.

Experimental Design
Evaporation experiments were conducted in the soil physics laboratory of the Kentucky Agricultural Experiment Station, University of Kentucky, Lexington, to derive soil hydraulic property functions (θ-h, K-h-θ) from pressure head and water content observations over time and at three soil depths. Two soil columns of 18 cm in length and an internal diameter of 10 cm were packed with air-dry silt loam soil (upper horizon from Bluegrass-Maury silt loam, classified as mixed, semiactive, mesic Typic Paleudalf, [33]) containing 11.8% sand, 70.3% silt, 17.9% clay, and 3.97% organic matter that had been passed through a 2-mm sieve. This soil was selected because it is commonly found in Kentucky and the southeastern United States, and does not exhibit volume changes with different water contents. The air-dry soil water content was 0.025 cm 3 cm −3 . One of the two soil columns was packed entirely with sieved soil while the other one was repacked with sieved soil and rock fragments. The cylinders were repacked under the premise of reaching a bulk density of 1.3 g cm −3 for the fine soil material. Rock fragments were distributed randomly and they occupied 40% of the total volume. The 40% volume of rock fragments was selected because this percentage falls inside the United States Bureau of Reclamation (USBR) arable land class number three. From this point on, subsequent land classes (four, five and six) show restrictions for irrigation purposes [34]. The rock fragments were sub-angular and sub-rounded shaped fragments of quartzite, with an equivalent diameter of 2-6 cm and with a negligible water holding capacity. The particle density of the rock fragments was measured by water displacement of their mass and it was of 2.65 g cm −3 . The rock fragments used had a mean mass of 15 grams, which means that the total soil mass in the cylinder was >100 times larger than the averaged mass of the rock fragments used. The soil column with rock fragments is further denoted as stony.
Three electronic pressure transducer tensiometers were inserted horizontally in the soil columns at 3, 8 and 13 cm from the surface. The insertion length of the ceramic cup was 4 cm and the diameter 0.6 cm. The holes for the ceramic cups in the soil were drilled using a screw auger with a diameter slightly smaller than the ceramic cup's diameter to ensure good cup-soil contact. The pressure transducers (26PCCFA6D) had been calibrated before. In order to indirectly measure the soil water content, CS640 TDR probes, with a length of 7.5 cm, were installed in the soil columns at the same depths as the ceramic cups. As shown in Figure 1, the cylinder columns were conceived as having 3 different compartments, based on the depth center between two measurement depths, 0-5.5 cm depth, 5.5-10.5 cm depth and 10.5-18 cm depth. content, CS640 TDR probes, with a length of 7.5 cm, were installed in the soil columns at the same depths as the ceramic cups. As shown in Figure 1, the cylinder columns were conceived as having 3 different compartments, based on the depth center between two measurement depths, 0-5.5 cm depth, 5.5-10.5 cm depth and 10.5-18 cm depth. Before the onset of evaporation, the soil column was placed on a porous plate that was connected to a hanging water column 15 cm in length. Soil columns were subjected to two wetting-drying cycles prior to the onset of the experiment in order to assure structural settlement.
The cylinder was covered with a lid to avoid evaporation and hydraulic equilibrium was established at a pressure head of −15 cm at the lower boundary. For the evaporation experiment, the soil column was removed from the suction plate and placed on a Plexiglas plate that was deposited on a balance (10 −2 g of precision). The balance was connected to a data logger and total mass was recorded once every hour for quantifying the amount of water evaporated. Next, the lid was removed from the top of the column. Evaporation continued until the upper tensiometer showed a pressure head reading close to −450 cm.
Pressure head and water content pairs were logged every minute and further processed to be used as inputs for quantifying the hydraulic properties of both soil columns.

Theory
The soil water content (θ) is a non-linear function of the soil water pressure head (h) that can be described by the closed-form equation proposed by van Genuchten [32]: Before the onset of evaporation, the soil column was placed on a porous plate that was connected to a hanging water column 15 cm in length. Soil columns were subjected to two wetting-drying cycles prior to the onset of the experiment in order to assure structural settlement.
The cylinder was covered with a lid to avoid evaporation and hydraulic equilibrium was established at a pressure head of −15 cm at the lower boundary. For the evaporation experiment, the soil column was removed from the suction plate and placed on a Plexiglas plate that was deposited on a balance (10 −2 g of precision). The balance was connected to a data logger and total mass was recorded once every hour for quantifying the amount of water evaporated. Next, the lid was removed from the top of the column. Evaporation continued until the upper tensiometer showed a pressure head reading close to −450 cm.
Pressure head and water content pairs were logged every minute and further processed to be used as inputs for quantifying the hydraulic properties of both soil columns.

Theory
The soil water content (θ) is a non-linear function of the soil water pressure head (h) that can be described by the closed-form equation proposed by van Genuchten [32]: where, θ s is the saturated water content (L 3 L −3 ), θ r is the residual water content (L 3 L −3 ), and α (L −1 ), n and m (m = 1 − 1/n) are the empirical fitting parameters describing the curve shape. The desorptive water retention curve obtained with the evaporation experiment was iteratively derived based on Wind´s method [28], using a fourth order polynomial equation ] + e), as described in Reference [30]. The hydraulic conductivity (K) was estimated for each time interval following the procedure described by Reference [29]. The hydraulic conductivity between the 0-8 cm depth and 8-18 cm depth compartments was calculated from the water flux across the 8 cm depth (q, Figure 1). Hydraulic gradients below a threshold value of 0.3 cm cm −1 were not considered for calculating K because of the precision of pressure transducer measurements [29]. Data obtained from the experiment were fitted to the unsaturated hydraulic conductivity function expressed in terms of the pressure head resulting from combining the Mualem [35] model for predicting unsaturated hydraulic conductivity with van Genuchten's model [30]: where K s is the saturated hydraulic conductivity and is a dimensionless pore tortuosity parameter that was assumed to be 0.5 as in many previous soil studies with and without stones [5,13,21], and to keep the number of parameters low for the inverse estimation solution.
It has been proposed that some hydraulic properties of stony soils can be derived from the fine earth fraction's properties if the water retention capacity of rock fragments is assumed to be negligible. Using the equation proposed by Reference [6] it is possible to obtain the bulk relationship of the soil water retention curve (SWRC, θ(h)) for a stony soil if the relation between the water content and pressure head of the same soil not containing stones is already known, using: where θ b (L 3 L −3 ) is the bulk volumetric water content of the stony soil, i.e., the volume of water per total volume of soil including the volume occupied by stones, V r is the volume fraction of rock fragments and θ fe (L 3 L −3 ) is the water content of the fine earth, i.e., the volume of water per volume of fine soil, excluding the volume occupied by stones. The equation to estimate the saturated hydraulic conductivity of a soil containing rock fragments based on K of the same soil not containing stones was derived by [10]: where K s,b denotes the saturated hydraulic conductivity of the stony soil and K s , fe the saturated hydraulic conductivity of the soil without stones. Hydraulic parameters defined in Equation (2) (α, n, m and ) for a stony soil are here assumed to be identical to those obtained for its fine earth as explained by Reference [5]. Thus, as pointed out by Reference [6], the unsaturated hydraulic conductivity function of the stony soil is assumed to have the same shape as the corresponding function of the fine earth but with a lower saturated hydraulic conductivity coefficient K s,b , calculated by Equation (4). This assumption has not been validated so far. A modification of Equation (4) was proposed by Reference [9] by introducing a parameter that corrects the non-linear relationship between stoniness and saturated hydraulic conductivity resulting from rock fragment resistance to water flow.
where a is a dimensionless parameter that depends on size, shape and rock fragments distribution, and soil texture.

Approach
In order to obtain the desorption branch of the soil water retention curve θ(h), as well as the unsaturated hydraulic conductivity as a function of pressure head (K(h)) for the soil without stones and the soil with rock fragments, different approaches were used in this study. The inputs used in every approach are summarized in Table 1. Table 1. Inputs used in every methodological approach employed to obtain the soil water retention curve and soil hydraulic conductivity function from the soil columns.

Soil Water Storage Evaporation Rate
Wind Method The results obtained from the numerical inversion in HYDRUS were parameter sets according to Equations (1) and (2).

Wind's Method
Wind's procedure [28] for determining the soil water retention curve (SWRC and K as a function of h is a transient method based on an evaporation experiment. An iterative procedure is carried out to estimate the SWRC based on h and the total water storage measured over time. No depth-specific measurements of soil water content over time are used as those were not available when the original method was developed. The iterative procedure is briefly described here. For further details see References [29,30,36]. The initial set of polynomial coefficients that are used for describing the θ(h) relationship was arbitrarily chosen. Given the pressure head measurements, water content values were calculated for each time and depth for their respective depth increments and converted to layer-specific water storage (L) by multiplying the computed θ values by their respective compartment thickness. The total soil water storage at a given time is the sum of soil water storage heights in the three different compartments. This calculated water storage was then compared to the soil water storage determined by weighing at the same time. Water contents estimated based on the polynomial coefficients were updated for each time by a correction factor which was the ratio between the measured and calculated total storage, assuming that the estimation error was equally distributed between compartments. Pressure heads were then plotted against the new set of water contents and new polynomial coefficients for the θ(h) relationship were obtained through fitting. This process was repeated until the maximum change in water content for all depths and times between two iterations was <10 −3 cm 3 cm −3 . Once the solution converged, the final water content values were adjusted by the correction factor one last time and used for calculating water fluxes between depth compartments ( Figure 1).
Accordingly, the fluxes and hydraulic gradients were applied to calculate the hydraulic conductivity for each time interval and depth compartment. The resulting θ-h-K data sets were fitted using the RETC code of Reference [37].

Inverse Estimation
The inverse solution of Richards´equation is a non-linear parameter optimization method. It minimizes an objective function Φ which describes the difference between the simulated and measured values of h and θ obtained at different depths over time during the evaporation experiment, using the Levenberg-Marquardt algorithm [38]. The inverse mode of the finite-element program HYDRUS-1D was used to estimate the van Genuchten parameters θ r , θ s , α, n and K s from the evaporation experiment data. The objective function was characterized either to pressure head and water content recorded at three different depths and at different times, or to pressure head and total soil water storage at different times. In the first case, the inverse solution parameter estimation was accomplished using tensiometer readings and water content measurements from the time-domain reflectometry (TDR) probes at the different depths. Results of this scenario were denoted as Estimation h-θ. In the second inverse estimation scenario, tensiometer readings and total water storage in the soil core measured by weighing at every time interval were used. Hence, depth-specific soil water content data measured with TDR probes were not used in this second inverse estimation scenario. Results obtained in this second scenario were defined as Estimation h-SWS.
The initial estimates of the soil hydraulic parameters were those associated to the average parameters for a silt loam soil selected from the HYDRUS-1D textural class catalog ( Table 2). The upper boundary condition for numerical simulations was the average evaporation rate measured during the evaporation experiment and the bottom boundary condition was set to zero flux. The magnitude of the evaporation rates measured experimentally was small, around 0.22 ± 0.05 cm d −1 (stony cylinder experiment) and around 0.16 ± 0.05 cm d −1 (stone-free cylinder). The fluctuation in the evaporation rates during the experiment was associated with small fluctuations in the temperature and the humidity in the laboratory. As these fluctuations did not change much in magnitude and were cyclic, the average evaporation rate was used as representative of the experimental conditions.

Conversion of Data from the Stone-Free Soil to Soil with 40% Rock Fragments
All equations presented so far are based on the assumption that the water retention characteristic and the hydraulic conductivity function of a stony soil are the same as those of its fine earth fraction but being uniformly corrected by the volume of rock fragments present. Based on this assumption, a new set of data was simulated to study if the influence of a rock fragment content of 40% on the soil hydraulic parameters was linearly related to the stone-free SWRC and K(h) functions by a correction factor of 0.4, or whether this relationship was non-linear due to rock fragment interference with water flow.
Equation (3) was used to convert the SWRC data from the stone-free soil to pairs of θ-h data representing a soil with 40% rock fragments. Hence, Vr of 0.4 cm 3 cm −3 represents a soil of 40% rock fragments. New water retention curves were obtained using this approximation. The resulting SWRC were then compared to those obtained using Wind's and inverse estimation methods from the experiment on the stony cylinder. A similar procedure was applied for K(h) functions obtained from the different methodologies to estimate the functions for a stony soil based on the results from the stone-free soil. The unsaturated hydraulic conductivity functions resulting for the stone-free soil were fitted according to Equation (2), and then converted for functions in the soil containing stones based on their saturated hydraulic conductivities using Equation (4).

Data Treatment
During the evaporation experiment, h and θ were measured every minute and total soil water storage every hour. In order to reduce noise from the TDR and tensiometric measurements, the data at a particular time of soil water storage measurement were smoothed by averaging four readings before and four readings after that time.

Soil Water Content Adjustment
Total water storage measured through weighing every hour was over or under estimated when calculated from TDR measurements. It is known that the default equation given by Reference [39] that relates both parameters is not exactly valid for any particular soil, and that a specific soil calibration may become necessary. It was assumed that the error of the TDR-measured soil water content was systematic and not constant but changing over the range of moisture conditions. Therefore, in order to minimize the error, water content values measured by TDR probes were adjusted based on the relationship between soil water content and h. The ratio between total soil water storage calculated from TDR measurements (S TDR ) and soil water storage obtained by weighing (S mass ) was correlated with the average log(h) of the soil column. In Figure 2a the relationship between this ratio (S TDR /S mass ) and the average log of h observed at the same time during the evaporation experiment are displayed for the stony and stone-free cylinders.
The ratio (S TDR /S mass ) differed depending on the soil column analyzed. For the column without rock fragments, the ratio slightly decreased as the soil became drier. However, the ratio for the stony cylinder showed the largest difference between water storage at the beginning of the experiment when the soil was wettest. This difference decreased as the soil became drier but towards the end of the experiment it increased again. Therefore, the ratio between the soil water storage calculated and weighed did not vary linearly with h for the stony soil column. This result led to the conclusion that soil water content measurements using TDR for the two soil columns had to be adjusted independently.
For this purpose, the (S TDR /S mass ) − log(h) relationship observations were fitted using a polynomial equation of second order: where S TDR /S mass (h) i,j was the adjustment factor for every time (i) and depth (j), and a, b and c were the polynomial coefficients. This approach led to a set of correction factors, one for each time interval and depth, used to obtain the adjusted θ values with Equation (7): where θ i,j is the water content measured by TDR for every time and depth, and θc i,j was the adjusted water content. The effect of the adjustment of TDR measurements on error minimization for the stone-free cylinder is shown in Figure 2b. With the adjustment of water content values, soil water storage based on sample mass and based on adjusted TDR measurements became very similar. The error, before and after the adjustment, is shown as the RMSE in Table 3. Using raw TDR measurements in the 40% rock fragments cylinder would have overestimated the total soil water storage. After the adjustment, the difference between mass and adjusted TDR based soil water storage became negligible ( Figure 2c, Table 3).

Pressure Transducer Tensiometer Adjustment
Before the evaporation experiment, the cylinders were left on a Plexiglas plate to reach hydrostatic equilibrium. After several hours of stability of the pressure reading, it was assumed that equilibrium was reached. At that point, the readings obtained from tensiometers installed at the different depths were expected to differ from each other by a height that reflected the 5-cm elevation difference between them. Deviations from this number were adjusted as described in references [29] and [36].
All h and θ results presented in the following sections were processed as described in 2.4.2 and 2.4.3.

Parametrization of the Stone-Free Cylinder
The results from assessing the validity of the Wind method and the HYDRUS inverse estimation for estimating the SWRC in a silty loam repacked stone-free soil column are shown in Figure 3. For

Pressure Transducer Tensiometer Adjustment
Before the evaporation experiment, the cylinders were left on a Plexiglas plate to reach hydrostatic equilibrium. After several hours of stability of the pressure reading, it was assumed that equilibrium was reached. At that point, the readings obtained from tensiometers installed at the different depths were expected to differ from each other by a height that reflected the 5-cm elevation difference between them. Deviations from this number were adjusted as described in References [29,36].
All h and θ results presented in the following sections were processed as described in Sections 2.4.2 and 2.4.3.

Parametrization of the Stone-Free Cylinder
The results from assessing the validity of the Wind method and the HYDRUS inverse estimation for estimating the SWRC in a silty loam repacked stone-free soil column are shown in Figure 3. For the three soil depths, θ(h) pairs as well as the average θ(h) set are depicted in Figure 3a. At the 3 cm depth, the soil water content was lower at the same h value than at the 8 and 13 cm depths. The θ(h) points resulting from averaging measurements from the three soil depths were closer to points representing the lower depth compartments.  An excellent agreement was found between the averaged θ(h) data measured with tensiometers and TDR probes during the experiment and the data points determined with the Wind [28] method. Their analytical description is shown in Figure 3b, and resulting parameters obtained with RETC (R 2 = 0.979) are listed in Table 4. Hydraulic parameters were also estimated using the inverse estimation method under the two different scenarios described in 2.3.2, i.e., h-θ and h-SWS, and the results are shown in Figure 3c and 3d, respectively. Soil water pressure heads and soil water contents measured at the three depths were used in the objective function under the inverse estimation h-θ scenario. The fitting between measured and predicted values resulted in R 2 equal to 0.989 (Figure 3c).
In order to analyze the quality of these estimated values, the correspondence between measured and simulated values was analyzed. Simulated water content values deviated from those measured mainly for the 3 cm depth and soil water pressure heads best matched the predicted values but An excellent agreement was found between the averaged θ(h) data measured with tensiometers and TDR probes during the experiment and the data points determined with the Wind [28] method. Their analytical description is shown in Figure 3b, and resulting parameters obtained with RETC (R 2 = 0.979) are listed in Table 4. Hydraulic parameters were also estimated using the inverse estimation method under the two different scenarios described in Section 2.3.2, i.e., h-θ and h-SWS, and the results are shown in Figure 3c,d, respectively. Soil water pressure heads and soil water contents measured at the three depths were used in the objective function under the inverse estimation h-θ scenario. The fitting between measured and predicted values resulted in R 2 equal to 0.989 (Figure 3c).
In order to analyze the quality of these estimated values, the correspondence between measured and simulated values was analyzed. Simulated water content values deviated from those measured mainly for the 3 cm depth and soil water pressure heads best matched the predicted values but deviated once pressure head reached −200 cm.
A more flexible approach for optimizing the hydraulic parameters by the inverse estimation h-θ was to divide the soil column into three depth compartments: 0-5.5 cm, 5.5-10.5 cm and 10.5-18 cm. As a result, three sets of soil hydraulic properties were obtained, one for each compartment. The agreement between simulated and measured results manifested in R 2 = 0.99, shown in Table 5, and was improved compared with the results obtained in the scenario assuming a homogeneous soil column (Table 4). In addition, the root-mean-square-error (RMSE = 0.025) was also improved implying that the goodness of fit between the measured and the estimated data was greater. Thus, water content and pressure head observations were better simulated by the inverse estimation h-θ approach, which allowed for vertical heterogeneity and dividing the soil cylinder into three compartments.    Figure 4a depicts the agreement between the soil water retention curve obtained by inverse estimation (estimation h-θ scenario) and the θ(h) points measured for every soil depth when the soil compartments were analyzed independently. Averaging the θ s values of the three different compartments resulted in 0.48, similar to the value obtained for the homogeneous soil scenario (Table 4). However, averaging parameter n resulted in a higher value than the one obtained for the unique parameter set scenario in which the cylinder was not divided. Parameter alpha was similar for each scenario studied.
In the inverse estimation scenario, tensiometer readings were used as inputs along with the total soil water storage determined by weighing the total mass at each measurement time (estimation h-SWS). The retention curve obtained by this optimization scenario corresponded well with the averaged θ(h) data points measured experimentally, as shown in Figure 3d. The value of R 2 = 0.996 indicated a better fit than for the inverse estimation scenario h-θ, which characterized the objective function using both h and θ measurements. This result could be explained by the fact that soil water measurements by TDR probes were not included in the h-SWS scenario. In this way, the optimization accuracy reduction associated to soil water content estimation when the cylinder was not divided into three compartments was avoided. The agreement between the measured and simulated soil water storage values using the inverse estimation h-SWS method was acceptable. Also, measured and simulated soil water pressure head values did not show severe discrepancies. The agreement between all the retention curves derived from the different estimation methods and the average θ(h) points measured during the evaporation experiment was good despite the fact that for the three depth compartments, different hydraulic parameter values were obtained. The numerical inversion resulted in lower predicted n values and slightly higher values of θs than those obtained by the Wind method.
The differences between various estimates of the parameter n could be caused by the fact that the range of measurements in our experiment was limited to the wet and medium range. Therefore, measured data did not support the estimation of this parameter with the accuracy that would have been achieved if our measurements had covered the range between water saturation and drier soil conditions more completely [30].
The unsaturated hydraulic conductivities determined by the different methodological approaches, considering the soil as a uniform column, are shown in Figure 5a. For the Wind method, no results for K could be obtained in the wet range as long as the threshold value of 0.3 cm cm −1 for the hydraulic gradient was not passed. Therefore, K(h) are only represented for the range of h < -50 cm. Unsaturated hydraulic conductivity functions obtained via parameter optimization with both inverse estimation scenarios (h-θ and h-SWS) were similar to the K(h) function obtained by the analytical description of the Wind method using RETC. The agreement between all the retention curves derived from the different estimation methods and the average θ(h) points measured during the evaporation experiment was good despite the fact that for the three depth compartments, different hydraulic parameter values were obtained. The numerical inversion resulted in lower predicted n values and slightly higher values of θ s than those obtained by the Wind method.
The differences between various estimates of the parameter n could be caused by the fact that the range of measurements in our experiment was limited to the wet and medium range. Therefore, measured data did not support the estimation of this parameter with the accuracy that would have been achieved if our measurements had covered the range between water saturation and drier soil conditions more completely [30].
The unsaturated hydraulic conductivities determined by the different methodological approaches, considering the soil as a uniform column, are shown in Figure 5a. For the Wind method, no results for K could be obtained in the wet range as long as the threshold value of 0.3 cm cm −1 for the hydraulic gradient was not passed. Therefore, K(h) are only represented for the range of h < -50 cm. Unsaturated hydraulic conductivity functions obtained via parameter optimization with both inverse estimation scenarios (h-θ and h-SWS) were similar to the K(h) function obtained by the analytical description of the Wind method using RETC.

Parametrization of the Stony-Soil Cylinder
Results for demonstrating the applicability of the different estimation approaches in the stony soil to obtain hydraulic parameters are presented in this section. Depth-specific and average θ(h) pairs obtained during the evaporation experiment for the stony soil are presented in Figure 6a. Averaging the measurement results of the three depths resulted in a similar water retention relationship as the one observed at the central position (8 cm depth).
The Wind method was applied to the three tensiometer and total water storage readings obtained over time for the stony soil. This iterative estimation gave a good water retention curve fit to a 4th degree polynomial equation with a R 2 value of 0.998. A promising result was also obtained after fitting data points obtained with the Wind method to van Genuchten´s model by RETC. The shape of the soil water retention curve derived from the fitted Wind results agreed well with the soil water retention points measured (Figure 6b). The results obtained from the different estimation approaches for the stony soil column are presented in Table 4.

Parametrization of the Stony-Soil Cylinder
Results for demonstrating the applicability of the different estimation approaches in the stony soil to obtain hydraulic parameters are presented in this section. Depth-specific and average θ(h) pairs obtained during the evaporation experiment for the stony soil are presented in Figure 6a. Averaging the measurement results of the three depths resulted in a similar water retention relationship as the one observed at the central position (8 cm depth).
The Wind method was applied to the three tensiometer and total water storage readings obtained over time for the stony soil. This iterative estimation gave a good water retention curve fit to a 4th degree polynomial equation with a R 2 value of 0.998. A promising result was also obtained after fitting data points obtained with the Wind method to van Genuchten´s model by RETC. The shape of the soil water retention curve derived from the fitted Wind results agreed well with the soil water retention points measured (Figure 6b). The results obtained from the different estimation approaches for the stony soil column are presented in Table 4.
The inverse estimation h-θ scenario was not as good as for the stone-free cylinder. Despite the fact the R 2 value between measured and simulated data was 0.997, the water retention curve shape derived from this inverse parameter estimation did not effectively represent the θ(h) measured points in the wet range (Figure 6c). Its RMSE value was the lowest for the different inverse estimation scenarios. The assessment of the quality showed that simulated water content values diverted from the 3-cm and 13-cm depth TDR measurements. These simulated values were similar to those measured at the 8-cm depth. The inverse estimation h-θ scenario was not as good as for the stone-free cylinder. Despite the fact the R 2 value between measured and simulated data was 0.997, the water retention curve shape derived from this inverse parameter estimation did not effectively represent the θ(h) measured points in the wet range (Figure 6c). Its RMSE value was the lowest for the different inverse estimation scenarios. The assessment of the quality showed that simulated water content values diverted from the 3-cm and 13-cm depth TDR measurements. These simulated values were similar to those measured at the 8-cm depth.
In the next step, the inverse estimation was performed for the three depth compartments 0-5.5 cm, 5.5-10.5 cm and 10.5-18 cm independently, in the same fashion as for the stone-free soil. The results of this scenario are provided in Table 5. Assuming depth-specific soil hydraulic property functions improved the estimation results, and the differences between predicted and measured water content were reduced. Figure 4b shows the agreement between the θ(h) points measured for every soil depth and the water retention curves shape obtained by inverse parameter estimation (estimation h-θ) when soil compartments were analyzed separately.
The inverse estimation h-SWS scenario (Figure 6d) gave the best R 2 result with a value of 0.999 and the lowest RMSE (Table 4). The measured and the simulated water content were in good agreement. In addition, measured and simulated soil water pressure head values were very similar except at the end of the experiment when small discrepancies were noted.
Resulting retention curves from the different estimation methods showed a close match with the average θ(h) points measured during the evaporation experiment except for the inverse estimation h-θ solution, where the curves deviated from each other near water saturation. The θs obtained was slightly lower compared to other solution values.
Differences existed between the estimation methods studied in terms of θr and n. Similar as for the stone-free soil, the lack of agreement for the n parameter value could be caused by the range of In the next step, the inverse estimation was performed for the three depth compartments 0-5.5 cm, 5.5-10.5 cm and 10.5-18 cm independently, in the same fashion as for the stone-free soil. The results of this scenario are provided in Table 5. Assuming depth-specific soil hydraulic property functions improved the estimation results, and the differences between predicted and measured water content were reduced. Figure 4b shows the agreement between the θ(h) points measured for every soil depth and the water retention curves shape obtained by inverse parameter estimation (estimation h-θ) when soil compartments were analyzed separately.
The inverse estimation h-SWS scenario (Figure 6d) gave the best R 2 result with a value of 0.999 and the lowest RMSE (Table 4). The measured and the simulated water content were in good agreement. In addition, measured and simulated soil water pressure head values were very similar except at the end of the experiment when small discrepancies were noted.
Resulting retention curves from the different estimation methods showed a close match with the average θ(h) points measured during the evaporation experiment except for the inverse estimation h-θ solution, where the curves deviated from each other near water saturation. The θ s obtained was slightly lower compared to other solution values.
Differences existed between the estimation methods studied in terms of θ r and n. Similar as for the stone-free soil, the lack of agreement for the n parameter value could be caused by the range of measurements and accentuated by the presence of rock fragments in the soil. As pointed out by Reference [30], θ r is associated with a high level of uncertainty if available measurements belong to the wet range. Similarly, Reference [18] reported larger errors in a modelling approach for stony soils at higher tensions attributing it to soil processes that modelling does not take into account.
Hydraulic conductivity versus pressure head results for the stony soil cylinder are presented in Figure 5b. Discrepancies became apparent when comparing the resulting unsaturated conductivity functions obtained from the stony cylinder data considering the soil as a whole compartment. Figure 5b reveals larger differences between the hydraulic conductivity functions obtained by both inverse parameter estimation scenarios and the Wind method compared with the stone-free cylinder results. Despite the fact that both methodological approaches use similar inputs, the difference between results obtained by the Wind method and by inverse estimation could be explained by the fact that the latter is a more restrictive approach for optimizing the hydraulic parameters. Thus, a more flexible approach was used to derive K(h) functions from both inverse estimation scenarios (h-θ, h-SWS) by dividing the soil column into three depth compartments. Three resulting sets of soil hydraulic parameters, one for each compartment, were geometrically averaged in order to obtain a unique set. The resulting K(h) functions were compared to those of the stone-free cylinder inverse estimation approach (Section 3.1).

Assessment of the Estimation of Hydraulic Properties in the Stony Soil
In this section the validity of estimating the stony hydraulic properties of soil from those of its fine earth by simply correcting the volume of rock fragments is analyzed.
The difference between the water content of the stone-free and the stony soil columns was not manifested in a constant shift as might be expected from Equation (3), but it depended on h. Under the conditions of this experiment it seems that Equation (3) (water content of the fine earth cylinder corrected by the volume of rock fragments) underestimated around 5% to 7% the bulk water content of the stony soil column mainly in the wetter range for all scenarios compared (Figure 7). In addition, this overestimation depended on the range of h and seemed to be weaker at more negative values.
One hypothesis for explaining this result could be the presence of rock fragments changing the pore-size distribution by creating voids and pores of big size that Equation (3) was not accounting for [16,18,40]. Different authors' findings support this hypothesis. Reference [22] estimated water retention curves from an alluvial gravel soil showing that retention curves derived from equations seem to underestimate the water content close to saturation, and thus the porosity. A faster decrease of water content in a stony soil with decreasing matric potentials than in a soil without rock fragments was reported by Reference [21], which attributed this result to the increase of pores larger than 0.25 cm. Also, References [11,27] found that soils with rock fragments resulted in a different fine pore structure than soils without.
Considering these results, the relationship between the water content of the stony soil and the water content of its fine earth fraction seems non-linear and dependent on h. In addition, the water content at saturation for a silt-loam soil with 40% of rock fragments was not exactly the water content at saturation of its fine earth corrected by the volume of rock fragments, but higher. Therefore, deriving the water content of a stony soil by using Equation (3) seems not correct as Reference [10] already acknowledged by suggesting that correcting the soil water content of a soil by simply using the rock fragment volume it contains could imply an over-simplification to represent the effect that these rock fragments exert, at least for determined types of soils like aggregated clays at the lower suction range of the water retention curve. Besides, Reference [5] applied Equation (3) to determine the bulk water content of a stony soil from the water content of its fine earth. They stated that this equation could only be used under the assumptions that rock fragments do not hold water, macropores are not present in the stony soil, and van Genuchten water retention parameters are identical to those of the fine fraction. In addition, Reference [8] pointed out that the air-entry value and the water retention curve shape coefficients of some repacked stony soils were not similar to those of the stone-free soil. Considering these results, the relationship between the water content of the stony soil and the water content of its fine earth fraction seems non-linear and dependent on h. In addition, the water content at saturation for a silt-loam soil with 40% of rock fragments was not exactly the water content at saturation of its fine earth corrected by the volume of rock fragments, but higher. Therefore, In our experiment, one solution to improve the accuracy of the estimation of the water content in the stony soil using data from the stone-free soils could be to complement Equation (3) as proposed by Equation (8): (8) where θ r,fe (L 3 L −3 ) is the residual water of the fine earth, θ s,fe (L 3 L −3 ) is the saturated water of the fine earth and b is a dimensionless parameter dependent on size, shape, rock fragment distribution and soil texture. Under the conditions of this evaporative experiment using rock fragments with negligible water retention capacity, a mean diameter between 2-6 cm, and sub angularly shaped stones distributed randomly in a silt loam soil, parameter b is equal to 0.875. The new equation proposed above (Equation (8)) includes the fitting parameters of van Genuchten´s model obtained for the water retention curve of fine earth but with θ s,fe and α adjusted by parameter b. This adjustment reduced to a minimum the differences between the SWRCs of the stony soil column and the SWRCs derived from the stone-free column using Equation (3) (Figure 7). Thus, Equation (8) considers the effects of the porosity change created by the rock fragments along the h range of measurements that Equation (3) does not account for. Figure 7 shows how the SWRC derived from the new Equation (8) was more similar to the SWRC obtained from the stony soil column (40% stones) with the Wind and HYDRUS estimation methods. There were some discrepancies at lower pressure head values, explained by some uncertainties in the adjustment of parameter n in the intermediate and wetter range of measurements. The agreement was also lower between the SWRC of the stony soil (40%) and the SWRC obtained using Equation (8) in the inverse estimation h-θ scenario due to less accurate results obtained by this approach when using a homogeneous soil compartment.
In relation to the hydraulic conductivity, Figure 8 shows the comparison between the hydraulic conductivities of the different data sets studied. The K(h) curves of the stony soil belonging to both inverse estimation scenarios (h-θ, h-SWS) were obtained dividing the cylinder into three depth compartments. For all the other cases shown in Figure 8, K(h) curves were derived considering the soil as a whole compartment.
Regardless of the methodological differences (see above in Section 3.2), two common features of the K(h) curve of the stony cylinder can be observed in Figure 8. First, the values of K derived from the stony (40% stones) and stone-free (0% stones) cylinders measurements in the higher range of pressure head values (h) were systematically higher than those derived from Equation (4) for the three different approaches. Second, the shape of the K(h) curves obtained for the stony cylinder denoted a faster decrease of K values with decreasing soil water pressure head than both K(h) curves from the non-stony cylinder and the curve determined using Equation (4).
These results support the hypothesis of rock fragments creating large pores that change the soil porosity distribution that Equation (3) and also Equation (4) are not able to account for. Large pores drain faster than smaller ones at lower h. However, when the water flow is reduced, their connectivity changes, decreasing K more rapidly. Similarly, Reference [22] measured K(h) in-situ by an infiltration experiment and found that K(h) values close to saturation were higher than those obtained using Equation (4). They stated that models used to estimate K(h) generated differences because they did not take into account structural features like macropores. Also, Reference [13] reported that rock fragments could have influenced the flux and the saturated conductivity of the soil by creating lacunar pores. In this infiltration experiment they observed that these pores were concentrated in the wetter range and, because of their size, were not able to transport water at low water pressure heads. This is also in agreement with Reference [12], who observed, during an infiltration experiment, that saturated hydraulic conductivities on stony soils were higher compared to stone-free ones, and Reference [18] observed a similar result for soils with 20 to 60% of rock fragments. However, that was not the case for the unsaturated hydraulic conductivity when h values became more negative. Both experiments had in common a resulting steeper decrease of K as h decreased in the presence of stones. lacunar pores. In this infiltration experiment they observed that these pores were concentrated in the wetter range and, because of their size, were not able to transport water at low water pressure heads. This is also in agreement with reference [12], who observed, during an infiltration experiment, that saturated hydraulic conductivities on stony soils were higher compared to stone-free ones, and reference [18] observed a similar result for soils with 20 to 60% of rock fragments. However, that was not the case for the unsaturated hydraulic conductivity when h values became more negative. Both experiments had in common a resulting steeper decrease of K as h decreased in the presence of stones. Therefore, the presence of larger pores able to conduct water mainly at higher h values (less negative) could be an explanation for the steeper decrease of K(h) in the stony soils column under all estimation scenarios compared in the present study. Furthermore, the effect of rock fragments on the wetting behavior and resulting K(h) could be a function of whether it is an upward capillarity-driven flow process or a gravity-driven infiltration process.
The percentage of volume occupied by rock fragments could be also a factor determining the effect that the presence of rock fragments exerts on a soil. Reference [25] tested how hydraulic Therefore, the presence of larger pores able to conduct water mainly at higher h values (less negative) could be an explanation for the steeper decrease of K(h) in the stony soils column under all estimation scenarios compared in the present study. Furthermore, the effect of rock fragments on the wetting behavior and resulting K(h) could be a function of whether it is an upward capillarity-driven flow process or a gravity-driven infiltration process.
The percentage of volume occupied by rock fragments could be also a factor determining the effect that the presence of rock fragments exerts on a soil. Reference [25] tested how hydraulic conductivity was affected by different stone contents. The findings revealed that contents higher than 30% of stones in soils could show higher saturated hydraulic conductivities than soils without stones. However, the steeper decrease of the hydraulic conductivity function was not so obvious. Reference [9] ran an inverse estimation simulation for different sizes and content of stones while using rounded stones. They identified that saturated hydraulic conductivity decreased with the rock fragment content depending on its size and soil type. However, Reference [8] found stony soil hydraulic conductivities higher than or similar to the stone-free soil depending on the volume occupied by the stones. On the other hand, Reference [41] found that rock fragments contents higher than 30% of volume contributed positively to the macroporosity. More similar to the conditions of the present study, Reference [26] reported that in their experiment a 40% rock fragment volume in a soil was the fraction at which the water movement in the soil was most benefited by the created macropores than restricted by tortuosity. These results agreed with the results obtained for the stony soil column at the wet range.
Therefore, applying Equation (4) for estimating the stony soil unsaturated hydraulic conductivity from the fine earth seems to be not sufficient for a silt loam soil with 40% of rock fragments with sizes from 2-6 cm. This is due to the effect of rock fragments on a soil not depending only on their volume, size, shape, origin or disposition as mentioned in previous studies, but also on the pressure head in which water flows are taking place.
The correction factor applied to the stoniness in Equation (5) suggested by Reference [9] could help to estimate the saturated hydraulic conductivity of a stony soil from its fine fraction by incorporating a parameter that takes into account the stones characteristics. However, the effect of rock fragments on the K(h) shape should be also addressed in order to infer the behavior of K through the whole pressure head range.
The results from the present study could be applicable to field measurements because even though the effect of macropores in repacked soils could be smaller than the effects in the field, the presence of rock fragments in repacked soil columns could also result in the creation of macropores if the structure of the cylinder is allowed to stabilize [22].

Conclusions
Nowadays, in a scenario of water scarcity and increasing food demand, the optimization of water use in agriculture is a major challenge. Stony soils are spread worldwide and conducting in situ experiments on them is complicated or near impossible. In these cases, inferring the parameters of water retention from the fine earth properties is the answer and should be accomplished as accurately as possible.
Under the experimental conditions used in this work, both the Wind method and the inverse estimation technique were successfully applied to estimate the soil hydraulic properties from an evaporation experiment on a silt loam soil even with a rock fragment volume of 40%. The range of suction in which the experiment was conducted, restricted to medium-wet, limited the accuracy of the estimation of the parameter n. Water retention curves derived from the inverse estimation scenario θ-h were improved when a more flexible approach, dividing the soil cylinder into layers, was used.
Results from this work enable us to understand that using a simple correction of the fine earth properties by the volume occupied by rock fragments to estimate the SWRC of a stony soil is a too simplified approach that does not represent reality. The effect that a volume of 40% of rock fragments of 2-6 cm exerted on a silt loam soil hydraulic properties was not a linear function of those of the same soil without rock fragments, and varied along the range of pressure head measurements.
A volume of 40% of rock fragments resulted in a positive influence for K in the wet range but this value decreased faster with decreasing pressure head, confirming the nonlinear effect of rock fragments on the behavior of K. Therefore, when the unsaturated hydraulic conductivity of a stony soil is derived from applying a correction to the fine earth values based only on the volume the rock fragments, the accuracy of the estimation will be low and will vary depending on the value of the pressure head considered.
Based on this research, a new equation to calculate the water content of a stony soil was proposed (Equation (8)) in which the key element for calculations was not only the volume of stones but the effect that the presence of lacunar pores and the pressure head has on water retention. Further research should be conducted to validate parameter b values for different types of soils, rock fragments volume, size, shape and origin in order to infer its validity. In addition, more studies should be conducted along the same lines that Reference [9] started by inferring K s values for stony soils with different textures and rock fragments properties but taking into account the different behavior of K values depending on h as we have observed in this study.