Geoderma 429 (2023) 116240

Contents lists available at ScienceDirect

Geodermajournal homepage: www.elsevier.com/locate/geoderma

Inverse estimation of hydraulic parameters of soils with rock fragments
Paulo Ivonir Gubiani a, *, Suélen Matiasso Fachi a, Quirijn De Jong Van Lier b,
Rodrigo Pivoto Mulazzani a, Fabrício de Araujo Pedron a, Jirka Šimůnek caSoil Science Department, Federal University of Santa Maria (UFSM), Av. Roraima 1000, Santa Maria, Rio Grande do Sul 97105-900, Brazil
Center for Nuclear Energy in Agriculture (CENA/USP), University of São Paulo, Piracicaba, São Paulo 13416-903, BrazilcDepartment of Environmental Sciences, University of California Riverside, Riverside, CA 92521, USAb

A R T I C L E I N F O

A B S T R A C T

Handling Editor: Morgan Cristine L.S.

The application of inverse modeling to determine the hydraulic properties of layered soils with rock fragments(RF > 2 mm) is associated with complex challenges such as selecting suitable hydraulic functions and defining astrategy to optimize many parameters. These two issues were addressed in this study by performing fielddrainage experiments in layered soil profiles (three layers) with different particle size distributions and different
percentages of RF. Water contents (θ) and pressure heads (h) were monitored at different soil depths duringdrainage experiments (h ≥ − 200 cm). Data of θ and h for a drier range (h < − 5000 cm) were determined with adew point potentiometer on disturbed soil samples. The unimodal and bimodal van Genuchten-Mualem functions
were evaluated, and their parameters were optimized using the Hydrus-1D “Inverse solution”. The RETC curve
fitting software was used to estimate the second set of parameters of the bimodal function. The parameters for theuni- and bimodal models (15 and 27 parameters, respectively) were repeatedly optimized by fixing some pa­
rameters while estimating others. The unimodal van Genuchten-Mualem function provided a satisfactory fit of θand h measurements only when either drainage or dew point potentiometer measurements were used. On theother hand, the bimodal van Genuchten-Mualem function provided a satisfactory fit of θ and h measurements
when both sets of data (drainage and dew point potentiometer measurements) were used simultaneously.
Furthermore, the optimized parameters were physically consistent except for a few high saturated hydraulicconductivity values. Therefore, bimodal functions should be considered to represent the hydraulic properties ofsoils with RF. The problem of many calibrated parameters when using bimodal functions can be managed using
sequential optimization, which allowed us to estimate 27 parameters for layered soils with RF successfully.

Keywords:
Drainage
Soil Hydraulic Functions
Heterogeneous Soils
Hydrus-1D

1. Introduction
Soils with rock fragments (RF > 2 mm) impose additional challenges
in assessing their hydraulic properties and hydrological processes
compared with soils without RF, mainly because of their high hetero­geneity and the impossibility of collecting undisturbed soil samples.
While worldwide research on several characteristics of these soils has
increased in the last decades, a consensual understanding of the effectsof RF on hydrological processes still demands further investigation(Zhang et al., 2016). One of the biggest challenges addressed in the lastyears is determining soil hydraulic parameters needed to model waterflow and contaminant transport through the vadose zone of these soils.
Direct laboratory measurements on undisturbed soil samples are com­plex and impractical depending on the amount of RF, which requires theuse of altered porous media created by mixing RF with finer soil particles(Ma and Shao, 2008; Beibei et al., 2009; Ma et al., 2010; Parajuli et al.,2017; Naseri et al., 2019). This strategy allows mimicking the effects of
RF on hydraulic functions, but replacing the natural porous system withan artificial one introduces bias in the soil hydraulic functions. Alter­
natively, an indirect determination of soil hydraulic parameters throughinverse modeling has been increasingly applied over the last decades (e.
g., Šimůnek et al., 1998; Lazarovitch et al., 2007; Jacques et al., 2012;
Rashid et al., 2015; Pinheiro et al., 2019; Brunetti et al., 2019) and isconsidered a feasible technique to estimate hydraulic parameters fromdata collected during non-destructive in-situ experiments in soils with
RF. However, challenges to field applications of inverse modeling toexperiments in complex porous media like those of soils with RF (Weiler
and Flühler, 2004) still remain a major research gap.
For any inverse modeling problem, optimized parameters depend on
the suitability and quality of (i) the experimental design, which* Corresponding author.
E-mail address: paulo.gubiani@ufsm.br (P. Ivonir Gubiani).
https://doi.org/10.1016/j.geoderma.2022.116240
Received 1 June 2022; Received in revised form 12 October 2022; Accepted 25 October 2022
0016-7061/© 2022 Elsevier B.V. All rights reserved.

P. Ivonir Gubiani et al.

Geoderma 429 (2023) 116240

comprises the boundary conditions, locations and time resolution of
measurement sensors, and the degree of accuracy of the experimental
data; (ii) the suitability of the transient variably-saturated water flow
model and hydraulic functions; and (iii) the robustness of the optimi­
zation algorithm as determined by its convergence to a global minimum(Hopmans et al., 2002). To prevent issues with numerical convergence,
the number of parameters being simultaneously fitted can be reduced(Ritter et al., 2003; Abbasi et al., 2003). A suitable model option is theHydrus-1D software (Šimůnek et al., 2013), which constitutes a robust,transient, variably-saturated flow model often applied on inverse
modeling problems (e.g., Lazarovitch et al., 2007; Jacques et al., 2012;
Rashid et al., 2015; Pinheiro et al., 2019; Brunetti et al., 2019). Thus, themain challenge, especially for soils with RF, is defining suitable soil
hydraulic functions and performing well-designed experiments.
A suitable method for determining parameters of hydraulic functions
for layered soils is the field (in-situ) drainage experiment, known as theinternal drainage or instantaneous profile method (Hillel et al., 1972;
Dikinya, 2005; Mavimbela and Van Rensburg, 2013). The method
consists of quasi-saturating a bare soil in which sensors (probes tomeasure water contents and pressure heads) have been installed atselected depths. The soil surface is then covered with a plastic sheet to
avoid any evaporation, and the soil is allowed to drain by gravity while
the transient drainage process is monitored by simultaneous measure­ments of soil water contents and pressure heads. Efforts to apply inverse
modeling in the field drainage experiment have focused on several issues
as, for example, the definitions of the hydraulic head gradient at thebottom boundary (Romano 1993), the likelihood of non-uniqueness andinstability of the inverse solution due to an increase in the number of
optimized parameters for layered soils (e.g., Zijlstra and Dane, 1996;
Ritter et al., 2003), and the combination of water content and pressure
head measurements (Dikinya, 2005). All these issues demand further
investigations, but evaluating the most suitable hydraulic functions forsoil with RF is an important and still unresolved investigation gap.
The choice of the hydraulic function shapes could be guided by priorknowledge about the soil pore size distribution. However, almost no
knowledge of pore size distribution is available for in-situ soils with RF.
Multimodal (soil hydraulic functions, having more than one maximum
in the soil water capacity curve and usually corresponding to a poresystem composed of more than one kind of pores, e.g., structural andintra-aggregate pores) rather than unimodal soil hydraulic functions,with only one maximum in the soil water capacity curve, may be moresuitable for soils with RF due to the presence of large pores and cracks
around RF and small textural pores in the aggregated fine material be­
tween them (Durner, 1994; Ma and Shao, 2008; Ma et al., 2010; Zhang
et al., 2016). Yet, multimodal functions increase the number of param­
eters to be optimized in the inverse problem, not only increasing the
likelihood of non-uniqueness of the inverse solution, but also requiring
measured water contents and pressure heads at a finer temporal reso­
lution so that the flow model can capture the effects of different poredomains on hydraulic functions.
However, water contents and pressure heads in field drainage ex­periments without evapotranspiration vary only over the narrow wet
range of soil hydraulic functions. The soil hydraulic parameters opti­
mized considering only such data may thus not represent complete poredomains and adequately define the behavior of soil hydraulic functions
over the dry range beyond the drainage phase. This limitation can be
avoided and/or minimized by defining the inverse solution using addi­
tional data involving water retention in the dry range (Hopmans et al.,2018; Pinheiro et al., 2019). Such data can be measured in the labora­
tory, for example, using a dew point potentiometer, an instrument that is
also suitable for soils with RF (Parajuli et al., 2017; Gubiani et al.,2021b). Thus, the shape of soil hydraulic functions can be better
determined by combining both data sets (in the wet and dry range) in the
optimization algorithm.
Avoiding ill-posed problems due to an increase in the number of soil
hydraulic parameters may represent the most significant challenge

when combining multimodal soil hydraulic functions with layered soil.
For example, combining the van Genuchten bimodal water retention
curve (Durner, 1994) with a soil profile with three soil horizons results
in a vector of 21 parameters (24 when including Ks) to be optimized. A
solution to this inverse problem with so many soil hydraulic parameters
may become ill-posed (e.g., Ritter et al., 2003; Abbasi et al., 2003),
mainly when local gradient-based optimization algorithms, such as theLevenberg-Marquardt scheme, are used (Šimůnek et al., 2012). To
address this problem, some parameters may be independently estimatedand their values fixed in the inverse solution, reducing the dimension­ality of the inverse problem (Abbasi et al., 2003). Another strategy is toperform sequential optimization (Ritter et al., 2003; Abbasi et al., 2003),
which repeatedly runs the optimization algorithm for selected parame­ters at each step.
Therefore, the main objective of this study is to determine the hy­
draulic parameters of soils with RF without disturbing the soil’s struc­ture while addressing the main investigation gaps discussed above. The
additional objectives of this study are to determine suitable hydraulicfunctions for soils with RF and find a sequential optimization strategy toestimate their parameters using inverse modeling and data from fielddrainage experiments and dewpoint potentiometer measurements.
2. Material and methods
2.1. Experimental site and soil characterization
Field drainage experiments were performed at four locations insouthern Brazil (Fig. 1). Locations were selected to include soil profiles
with different particle size distributions and different percentages of RF,varying from a soil profile with only a few particles > 2 mm (Location 2)
to a soil profile with a few particles < 2 mm (Location 3). Soils wereclassified as Typic Udorthent (Location 1 and 4), Humic Dystrudept(Location 2), and Lithic Udorthent (Location 3), according to Soil Tax­
onomy (Soil Survey Staff, 2014), formed from igneous rocks (rhyolite or
rhyodacite). Based on typical soil morphological properties such as
texture and color, three layers were defined in each profile, which also
defined the material depths in the discretized soil profile in Hydrus-1D(see Section 2.4.1). Location 1: 0–20, 20–40, and 40–65 cm; Location
2: 0–10, 10–35, and 35–55 cm; Location 3: 0–10, 10–25, and 25–40 cm;
and Location 4: 0–15, 15–35, and 35–55 cm.
Disturbed soil samples were collected from each layer to characterize
particle size distribution (Fig. 1), determined according to Schoene­
berger et al. (2012). Due to the presence of relatively fragile sand par­ticles in these soils, the samples used for the determination of sand, silt,
and clay (20 g of the dry soil fraction passed through a 2 mm mesh sieve)
were shaken for 2 h in a solution of 1 mol/L NaOH (Gubiani et al.,2021a). The clay fraction was determined with the pipette method (Geeand Bauder, 1986), and the silt fraction was calculated as the remaining
part by subtracting sand and clay from the total amount.
2.2. Instrumentation
Capacitance sensors with a factory calibration (ECH20 10HS Large
Volume Soil Moisture Sensors by Meter, Inc.) to measure water contents
were horizontally installed in the soil profiles at depths of 5, 20, 40, and65 cm at Location 1, depths of 5, 20, 35, and 55 cm at Location 2, depths
5, 17, and 40 cm at Location 3, and 5, 20, 35, and 55 cm at Location 4.
Additionally, soil water potential sensors with a factory calibration(Teros 21 by Meter Inc.) were installed at Locations 2 and 4 (at depths of
5 and 20 cm) and Location 3 (5 cm). Note that the sensors were notinstalled at the center of each layer (see Section 2.1), and some layers
had more than one sensor. This is because the sensors’ depths were
defined according to other studies’ needs on water balance. A detailed
description of the use of water content and pressure head measurements
is provided in Section 2.4.1. After the sensor installation on November 7,
2020, soil trenches were backfilled with the original soil materials of2

P. Ivonir Gubiani et al.

Geoderma 429 (2023) 116240

Fig. 1. Locations of experimental sites (top) and particle size distributions (bottom) of evaluated soil profiles (middle). Coarse gravel (CCG, 250 – 20 mm), medium
and fine gravel (MFG, 20 – 2 mm), coarse sand (CS, 2–0.25 mm), fine sand (FS, 0.25–0.05 mm), silt (0.05 – 0.002 mm), and clay (<0.002 mm).

each layer, followed by compaction to minimize preferential flow
(Fig. 2).

The experiments were performed shortly (two or three days) after
significant rainfall events when the initial water content was close to
field capacity. This guaranteed that the wetted soil volume was sur­rounded by wet soil, avoiding considerable lateral flow during thedrainage experiments. No rainfall occurred shortly after the onset of theexperiments, allowing enough time for drainage rates to decrease to
negligible values. The duration of observations was 50 h for the exper­
iments on November 30, 2020 (a rainfall occurred after 50 h), and 180 h
for the experiments on December 3, 2020, April 25, 2021, and May 22,2021.

2.3. Drainage experiments
Four drainage experiments were performed at each location on
November 30, 2020, December 3, 2020, April 25, 2021, and May 22,2021. The area was planted with soybean until March 15 (during thefirst two experiments) and with oat after April 10 (during the last two
experiments). During all experiments, plants were at their initial growthstage, and they did not interfere in the experimental determinations.
A circular area of approximately 2 m2 (sensors at its center) was
delimited with a 20-cm high plastic sheet, 10 cm above the soil surface
and 10 cm buried. Water was ponded on the surface of the boarded soil,
and an increase in the soil water content was monitored with 10 HS
sensors every 2 min. Water supply was ceased after five similar readingsof the highest water content were observed at all depths. When all
ponded water had infiltrated, the soil surface was covered with a plastic
sheet to prevent evaporation. Water contents and pressure heads werethen measured over time while redistribution and drainage took place(Fig. 2).

2.4. Definition of hydraulic functions and parameter optimization by
inverse modeling
Parameter optimization was performed using the inverse routine of
the Hydrus-1D model (Šimůnek et al., 2013). The one-dimensional
Richards equation solved by Hydrus-1D for vertical drainage flow isgiven as follows:
[ (
)]
∂h ∂
∂h
C =
K
+1
(1)
∂t ∂z
∂z
3

P. Ivonir Gubiani et al.

Geoderma 429 (2023) 116240

Fig. 2. Schematic of the timeline of the setup of drainage experiments (steps 1, 2, and 3) and the data analysis using parameter optimization (steps 4, 5, and 6).

where C is the water capacity (dθ/dh, cm− 1), θ is the volumetric watercontent (cm3 cm− 3), h is the soil pressure head (cm), t is time (min), z is
the vertical spatial coordinate (cm), and K is the hydraulic conductivity(cm min− 1).
From several hydraulic functions available in Hydrus-1D to describethe θ-h-K relations in equation (1) we selected one of the most common,a unimodal water retention function (van Genuchten, 1980), and its
parameters were optimized considering the measured and estimated θand h in the wet range of the drainage experiment. Subsequently, weused an expansion of this unimodal function to a bimodal one (Durner,
1994) to reduce inconsistencies that would prevent the extrapolation of
the unimodal function to drier domain. These strategies are described indetail in the next sections.

Se1 =

1− 1
θ − θr1= [1 + (α1 h)n1 ]n1
θs1 − θr1[l
K = Ks1 Se1

1−

(

n1
n1 − 1

1 − Se1)1− n1 ]2
1

(2)

(3)

where θs1 and θr1 are the saturated and residual water contents (cm3
cm− 3), α1 (cm− 1) and n1 (–) are fitting shape parameters, K and Ks1 (cmmin− 1) are the hydraulic conductivities for unsaturated and saturated
conditions, Se1 is the effective saturation (–), and l (–) is the poreconnectivity parameter. Subscript “1” distinguishes the parametersoptimized for a unimodal curve from those optimized for a bimodal
curve. Hysteresis was disregarded in all simulations.
During the inverse modeling, the parameters θs1, θr1, α1, n1, and Ks1
were sequentially optimized (see below), while the pore-connectivity
parameter l was kept at 0.5 according to Mualem (1976). All soil pro­files (Fig. 1) were represented in Hydrus-1D by a discretized schemewith 101 equally spaced (1 cm) nodes. Each soil profile consisted of

2.4.1. Parameter optimization of the unimodal water retention curve
The van Genuchten-Mualem functions (van Genuchten, 1980) wereused to describe the θ-h-K relation in equation (1).

4

P. Ivonir Gubiani et al.

Geoderma 429 (2023) 116240

three materials to represent the sequence of soil layers (see Section 2.1).
The upper boundary condition was set as a zero flux (a covered soil
profile), and the lower boundary was set as free drainage, in agreementwith field conditions where no water table or impeding soil layers were
present within several meters of depth. A pressure head of − 10 cm wasassumed as the initial condition along the entire domain. Although the
Teros 21 sensor cannot accurately measure pressure heads higher than
− 100 cm, a value close to saturation was considered a reasonableapproximation of the field-saturated condition, close to the air-entry
value for coarse-textured soils (Rawls et al., 1982). The suggested
default values by the Hydrus-1D model were used for the water content
error tolerance (0.001 cm3 cm− 3), pressure head error tolerance (-1 cm),and the maximum number of iterations (1 0 0) for the numerical solutionof the Richards equation. Hydrus-1D automatically selects an optimal
time step, allowed to vary between the minimum and maximum limits
set as 0.0144 min and 7200 min, respectively.
Parameters θs1, θr1, α1, n1, and Ks1 were optimized by minimizing theobjective function (Šimůnek et al., 1998):
nqj
∑

mq
∑

Φ(p, b) =

vj
j=1

[
]2
wi,j q*j (x, ti ) − qj (x, ti , b)

At the locations where no measurements of the pressure heads were
performed, estimates of h were used. At Location 1, the pressure heads h
at depths 5 and 20 cm were considered equal to the mean h value atcorresponding depths at Locations 2 and 4. The pressure heads h atlarger depths without measurements were considered the same or about− 3 to − 10 cm less than h at the overlying position. This strategy facil­
itated the setting of initial values for hydraulic parameters before
running the automatic Hydrus-1D “Inverse solution” routine (see step 2
in the next section). Observations points were defined in the discretizedprofile scheme of Hydrus-1D to obtain estimated qj(x, ti, b) values cor­
responding to measured qj*(x, ti) values.
2.4.2. Sequential optimization
A sequential calibration was implemented to reduce nonconvergence of the optimization algorithm and the likelihood of nonuniqueness when a large number of parameters had to be optimized(e.g., Hopmans et al., 2018; Ritter et al., 2003; Abbasi et al., 2003). The
entire vector of 15 parameters (θs1, θr1, α1, n1, and Ks1 for each of thethree soil layers) was reduced to a few selected parameters to be cali­brated in different steps as described below.

(4)

i=1

1. The initial estimate of θs1, θr1, α1, n1, and Ks1 for each material (soillayer) was taken from the Hydrus-1D catalog for a loam soil (θr1 =0.078, θs = 0.4, α1 = 0.036 cm− 1, n1 = 1.56, Ks1 = 1.7⋅10-2 cmmin− 1).
2. Before running the Hydrus-1D “Inverse solution” routine, most pa­
rameters were several times updated in Hydrus-1D by a trial-anderror procedure. Parameters were approximately refined by consid­
ering measured data of θ and visual inspection of the similarity be­
tween the measured and simulated time series of θ in the Hydrus-1D
graphical interface. A good improvement in the fit between
measured and simulated data was obtained by trial-and-error when
considering θs1 and θr1 close to the initial and final values, respec­tively, of measured θ in the drainage experiments. Subsequently, α1and n1 were also modified to improve the similarity between simu­
lated and measured water content values. All these improvements
were best achieved with Ks1 fixed at 0.2 cm min− 1.
3. Parameter optimization with the Hydrus-1D inverse solution routine
was accomplished in sequential steps to optimize only a few opti­
mized parameters in each step, thus minimizing ill-posed problems
and avoiding non-uniqueness and unidentifiability. First, only α1 andKs1 were optimized, and they were forced to have the same values forall soil layers. This was done by manually setting Kod = 2 (a code
indicating that the same parameter value should be used for all soil
layers) in the input file Fit.in. Then, the initial estimates of Ks1 and α1
were updated with their respective optimized values. Next, only θs1and θr1 for all layers were optimized (with Kod = 1, allowing
different parameter values in different layers). After that, the initial
estimates of θs1 and θr1 for all layers were updated with optimized
values. In the next step, α1, n1, and Ks1 for all layers (with Kod = 1)
were optimized and updated. Finally, a last minor refinement wasdone for θs1 and θr1, optimizing them for all layers. Thus, the numberof optimized parameters increased from two to six and then to nine,respectively, for the first, second, and third optimization steps. The
number of optimized parameters in the third step (nine) was higher
than usually recommended to avoid ill-posed problems (Ritter et al.,
2003; Abbasi et al., 2003). However, it is expected that the initial
trial-and-error steps followed by the first and second steps of inverse
modeling located the search algorithm in the parameter space quite
close to the global minimum of the objective function, making thethird step of the inverse solution less prone to ill-posedness, evenwith nine parameters being optimized simultaneously. Restrictions
were imposed on n1 (1.1 ≤ n1 ≤ 3.5) and θs1 and θr1 in cases of nonconvergence.

where qj*(x, ti) represents the specific measurements of θ and h (mea­surement sets j) at time ti at location x, qj (x, ti, b) represents the cor­
responding model predictions of θ and h for the vector of optimized
parameters b (θs1, θr1, α1, n1, and Ks1), vj and wi,j are weights associatedwith a particular measurement set or point, respectively, mq is thenumber of different sets of measurements (2 in this case, for watercontents and pressure heads), and nqj is the number of measurements
within a particular measurement set. The weighting coefficients wi,jwere all fixed at 1, giving the same weights to all data points in aparticular measurement set. Furthermore, the weighting coefficients vj
were fixed at (Šimůnek and Hopmans, 2002a,b):vj =1
nqj σ 2qj

(5)

which implies that the objective function Φ (equation (4)) is the average
weighted squared deviation between the measured and fitted values
normalized by the measurement variances (σ2).
The Hydrus-1D model additionally reports various other statistics,
such as Mean Weighted Error (ME), Mean Weighted Absolute Error(MAE), Root Mean Square Weighted Error (RMSWE), Akaike Informa­tion Criterion (AIC), and Bayesian Information Criterion (BIC). The Root
Mean Square Weighted Error (RMSWE), Akaike Information Criterion(AIC), and Bayesian Information Criterion (BIC) are defined in Hydrus1D as follows.
√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
Φ(p, b)RMSWE =
(6)
n− 1
(
)
Φ(p, b)
+ 2m
n− 1

(7)

BIC = (AIC − 2m) + mlnn

(8)

AIC = nln

where the sum of squared weighted deviations Φ is calculated using (4),n is the number of data points, and m is the number of optimized pa­
rameters. The Akaike and Bayesian information criteria (Akaike, 1974,1977) penalize for adding fitting parameters, i.e., the model with thesmallest number of parameters (which minimizes AIC and BIC) is apreferred fitting model (Šimůnek and Hopmans, 2002b).
Data for qj*(x,ti) consisted of the time series of θ measured during thedrainage experiments with the 10 HS sensors and h measured with the
Teros 21 sensor at the end of the drainage experiments (locations ofmeasurements are described in Section 2.2). Only the final measured h
was used because Teros 21 readings are insensitive at wetter conditions.
5P. Ivonir Gubiani et al.

Geoderma 429 (2023) 116240

2.4.3. Expanding the unimodal to bimodal water retention curve
Although the simulated θ and h matched their respective measure­ments quite well (RMSWE < 0.05 in most cases), the optimization pro­
cess made the unimodal water retention curve overestimate the soil
residual water content. Optimized values of θr1 ranged between 0.26 and0.44 cm3 cm− 3. These values are inconsistent for the coarse-textured
soils under investigation (Fig. 1). Overestimation in optimized valuesof θr1 could be a consequence of using only the θ and h values measured
in the wet range in the objective function.
We attempted to eliminate this problem by also considering the dry
range θ and h data measured in the laboratory (h < − 5000 cm) with adew point device (Gubiani et al., 2021b) in the objective function. This
made optimized values of θr1 decrease, but the agreement between
simulated and measured θ and h for the drainage experiments became
significantly worse. Therefore, keeping the unimodal θ-h-k relations
with the parameters optimized using θ and h measured in the wet range
would provide accurate drainage predictions, but its use to simulate
fluxes in a drier soil (e.g., root water uptake, infiltration) would mostprobably lead to inaccurate results.
We, therefore, concluded that the Hydrus-1D algorithms for waterflow and parameter optimization yielded unimodal water retention andhydraulic conductivity functions, which were unable to describe soil
hydraulic properties for the entire saturation range. A likely reason isthat the porous system of the soils evaluated in this study is more
complex than the unimodal model would allow, which is possibly aconsequence of the presence of coarse material (Fig. 1). Thus, theunimodal θ-h-k relations were expanded to bimodal functions (Durner
1994):
}
{
1− 1
1− 1
θ = θr2 + (θs1 − θr2 ) w1 [1 + (α1 h)n1 ]n1 + w2 [1 + (α2 h)n2 ]n2
(9)

{

[

(w1 Se1 + w2 Se2 )l w1 α1 1 −K = Ks2

(

n1
n − 1

1 − Se11)1− n1 ]1[+ w2 α 2 1 −

(

n2
n − 1

optimization was also avoided due to the impossibility of evaluating theuncertainty associated with all optimized parameters. This problem isaddressed in the discussion section.
2.5. Evaluation of results
As the optimization performed by Hydrus-1D is the focus, we
numbered the sequence of optimization steps from zero to 5, with
optimized parameters shown in parenthesis: step zero (θs1, θr1, α1, n1,
and Ks1), step 1 (α1 and Ks1), step 2 (θs1 and θr1), step 3 (α1, n1, and Ks1),
step 4 (θs1 and θr1 again), and step 5 (Ks2). A trial-and-error procedure
was used only in step zero. Note that in steps zero to 4, parameters of theunimodal curves (equations (2) and (3)) were optimized, while for thebimodal curves (equations (6) and (7)), only Ks2 was optimized. The
standard errors of parameters w1, w2, θr2, α2, and n2 of the bimodal
curves fitted with the RETC are presented. The performance of inverse
modeling was evaluated using the Root Mean Square Weighted Error(RMSWE) between measured and estimated variables (θ and h), andAkaike (AIC) and Bayesian (BIC) Information Criteria, which asses theadequacy of the model given the increasing number of optimized pa­
rameters. Water retention curves for each location and layer weresummarized in the “averaged” water retention curve, but all estimatedcurves are graphically presented.
3. Results
The sequential optimization gradually reduced the objective func­tion (and RMSWE of measured versus estimated θ and h) to values lower
than 0.05 of the initial value at the final step (Fig. 3). Supplementary
data (Table S.1) provides all values of RMSWE, AIC, and BIC, for all
sequential optimizations carried out in this manuscript. Table S.1 shows)1− n1 ] }22

1 − Se22(10)(w1 α1 + w2 α2 )where θs1, α1, n1, l, and Se1 are parameters described with equations (2)and (3), w1 and w2 (dimensionless) are weighing factors for each sub
curve (conditioned to w1 + w2 = 1), and θr2, α2, n2, Se2, and Ks2 have thesame meaning as their corresponding parameters in equation (2) and(3), but now for smaller textural pores.
Optimizing all parameters of equations (9) and (10) with inverse
modeling would likely be unfeasible and highly dependent on theirinitial guesses. We assumed that the best initial guess could be provided
by keeping θs1, α1, n1, and m1 at their previously optimized values(Section 2.4.2) and estimating w1, w2, θr2, α2, and n2 by the usual pro­cedure of fitting water retention functions (Fig. 2). Thus, parameters w1,
w2, θr2, α2, and n2 were estimated using the RETC software (van Gen­uchten et al., 1991). In addition to θ and h data measured in the dry
range in the laboratory (h < − 5000 cm) with a dew point device(Gubiani et al., 2021b), the wet range θ and h data (h ≥ − 200 cm)
estimated using equation (2) were used to avoid estimates of w1 and w2that would change the final wet branch provided by equation (2). The
limit h = − 200 cm is slightly smaller than the most negative matricpotentials measured in the field at the end of the drainage experiments(between − 140 and − 90 cm).
Finally, the inverse solution routine of Hydrus-1D was run once more
to find the saturated hydraulic conductivity value (Ks2) of equation (10)without any additional alteration in all other parameters. Further

Fig. 3. Root Mean Square Weighted Errors (RMWSE, dimensionless) of simu­
lated versus measured θ and h after each optimization step. Soil hydraulic pa­
rameters θs1, θr1, α1, n1, and Ks1 were defined by trial and error in step zero,
while the Hydrus-1D inverse solution routine was used in all other steps tosequentially optimize α1 and Ks1 (step 1), θs1 and θr1 (step 2), α1, n1, and Ks1(step 3), θs1 and θr1 again (step 4), and finally Ks2 (step 5). For details, see
Sections 2.4.2 and 2.4.3.
6

P. Ivonir Gubiani et al.

Geoderma 429 (2023) 116240

Fig. 4. Measured (solid lines) and simulated (dotted lines) water contents with calibrated parameters at different depths of four locations (rows) during four drainage
experiments (columns).

how sequential optimizations for each location and each drainage
experiment gradually reduced the RMSWE value, as well as both infor­
mation criteria, indicating that each sequential optimization provided abetter description of the experimental data as well as a better fitting
model. As a result, the estimated θ after the fifth optimization step
matched quite well the measured θ during the internal drainage exper­
iments on all dates, locations, and depths (Fig. 4). The greatest

deviations between estimated and measured θ were associated withwater content variations during drainage. The greatest deviation can beseen for Location 2 (Fig. 4), where the lowest and highest amplitude of θ,
considering all layers and dates, were 0.09 and 0.14 cm3 cm− 3. A
decrease in deviations between estimated and measured θ can be seen at
Locations 1, 3, and 4 (Fig. 4), where a wider amplitude of θ occurred.
The lowest and highest amplitude of θ was 0.04 and 0.14, 0.17 and 0.22,
7

P. Ivonir Gubiani et al.

Geoderma 429 (2023) 116240

which was needed to properly connect the wet and dry segments, isnoticeable in all curves (Fig. 7).
4. Discussion
One of the main problems addressed in this manuscript was to findsuitable hydraulic functions for soils with RF. By using θ and h measured
during field drainage experiments, were found parameters of unimodal
functions (equations (2) and (3)) that yielded simulated water contents
quite close to their respective measurements during the entire drainage
experiment. However, as discussed in the methodology section (Section
2.4.3), values of θr1 (between 0.26 and 0.44 cm3 cm− 3) were inconsis­
tent for the coarse-textured soils under investigation (Fig. 1) and would
lead to inaccurate simulations of fluxes in drier soils (e.g., root water
uptake, infiltration, etc.). We then used the laboratory-measured
retention data for h < -5000 cm to correct the dry end of the retention
curve and replaced the unimodal functions (equations (2) and (3)) withthe bimodal functions (equations (6) and (7)) to keep estimated watercontents close to those measured in the wet range.
We can therefore conclude that bimodal functions are more suitable
for soils with RF than unimodal functions taking into account thefollowing evidence: (i) bimodal functions yielded simulated water con­tents quite close to measured ones during the drainage experiments(Fig. 4), (ii) deviations between simulated and measured θ were lower at
the locations with higher contents of RF (1, 3, and 4) than at the location
with a smaller content of RF (location 1) (Fig. 4), (iii) optimized pa­
rameters were physically consistent except for a few high values of Ks2(Fig. 6), (iv) dry ends of retention curves (Fig. 7) matched well datameasured for h < − 5000 cm (from 0.02 to 0.28 cm3 cm− 3), and (v)multimodal functions are commonly used for soils with RF due to thepresence of larges pores and cracks surrounding RF, which contrast with
textural pores in the aggregated fine material among them (Durner,
1994; Ma and Shao, 2008; Ma et al., 2010; Zhang et al., 2016). As θ and hwere not measured for the transition interval between the two porous
systems (around h − 200 cm, Fig. 7), it was impossible to evaluate how
well the estimated parameters (Fig. 6) describe the real porous system.
Future studies can address this issue by adding to the objective function
retention data (θ and h) measured beyond the drainage experiment,
allowing for water losses due to evaporation and root water uptake
when the area is cultivated.
Another challenge addressed in this study was to perform parameter
optimization while avoiding ill-posedness and instability problems of
the inverse solution, which needs more attention when local gradientbased optimization algorithms, such as the Levenberg-Marquardt
scheme, are used (e.g., Šimůnek et al., 2012). Usually, less than fiveparameters are optimized simultaneously to avoid ill-posedness andinstability problems of the inverse solution (Ritter et al., 2003; Abbasi
et al., 2003). Our strategy of sequentially optimizing the vector of 15
parameters (θs1, θr1, α1, n1, and Ks1) for three soil layers successfully
avoided failure of the Hydrus-1D inverse solution due to non-Fig. 5. Measured and estimated pressure head (h) for all locations, dates
and depths.

and 0.13 and 0.16 cm3 cm− 3, respectively, at Locations 1, 3, and 4.
Estimated pressure head also matched quite well their corresponding
measurements at the end of drainage (Fig. 5). These data were presented
in a scatter plot because only the pressure head measured at the and ofdrainage were used (see the last paragraph of Section 2.4.1). The highest
absolute difference between measured and estimated pressure head was
3 cm.
The estimated parameters (Fig. 6) are located in a physicallyconsistent range, except for some extremely high values of Ks2 (6.0E+4,
1.8E+5, 5.4E+5, and 1.4E+6 cm h− 1), while all other Ks2 values werelower than 3.6E+3 cm h− 1. The high and apparently inconsistent valuesof the optimized residual water contents for the unimodal curve (θr1 hadranged from 0.26 to 0.44 cm3 cm− 3) were reduced to expected values,
not higher than 0.12 cm3 cm− 3 (θr2 in Fig. 6), when the bimodal function(equation (6)) was used. These values of θr2, together with w2, α2, and n2(fitted using RETC, Fig. 6), yielded a final descending segment of thewater retention curve (Fig. 7), which is more realistic than the hori­zontal segment of the unimodal function (Fig. 2) for h lower than − 200
cm (θ between 0.26 and 0.44 cm3 cm− 3). Supplementary data
(Table S.2) provide all estimated parameters for each location and eachdrainage experiment and statistical measures such as Standard Error, tvalue, and 95 % Confidence Limits. The four water retention curves ofeach layer differed more toward the wet end (Fig. 7), because parame­ters θs1, α1, and n1 of each curve were optimized with independent data,that is, θ and h measured during drainage experiments on November 30,
2020, December 3, 2020, April 25, 2021, and May 22, 2021. On theother hand, the data used in the dry range (θ and h data measured for h< -5000 cm with a dew point device) were the same for all curves. The
segment indicating a transition between two different pore systems,

Fig. 6. The range of all optimized parameters at all locations, depths, and times. Only w2 is presented because w1 is complementary to w2 (w1 + w2 = 1,
equation (7)).
8P. Ivonir Gubiani et al.

Geoderma 429 (2023) 116240

Fig. 7. Bimodal water retention curves. Gray lines in each graph represent individual retention curves for particular drainage experiments (four dates). Read lines
represent the “averaged” curve. Open circles are the measured data with a dew point device.

convergence and gradually reached the minimum of the objectivefunction (RMSWE of measured versus estimated θ and h).
After the initial trial and error approach (step zero), the sequential
optimization of α1 and Ks1 (step 1), θs1 and θr1 (step 2), α1, n1, and Ks1(step 3), and again θs1 and θr1 (step 4) reduced the RMSWE to < 0.05 of
the initial value (Fig. 3) after step 4, which was the last optimization stepfor parameters of unimodal functions (equations (2) and (3)). In the last
step (step 5), in which only Ks2 for bimodal functions was optimized
using the Hydrus-1D inverse solution, the RMSWE was reduced inseveral cases even further (Fig. 3). Thus, by combining sequential

optimization with bimodal hydraulic functions, in which only a few
parameters are optimized while the others are kept constant (Ritter
et al., 2003; Abbasi et al., 2003), we successfully performed inverse
parameter optimization for layered soils with RF.
The lack of data for the transition interval between the two porous
systems (from around h = − 200 cm to h = − 5000 cm, Fig. 7), as dis­cussed earlier, prevented the optimization of w1, w2, θr2, α2, and n2
directly using the Hydrus-1D inverse solution and thus these parameters
were optimized using RETC. However, even if measured data for thisinterval were available, fitting all parameters of bimodal functions
9

P. Ivonir Gubiani et al.

Geoderma 429 (2023) 116240

described by equations (9) and (10) for multilayered soil profiles would
remain challenging. A three-layered soil profile would result in a vector
with 24 soil hydraulic parameters. It would be unavoidable to fix someparameters to measured values or estimate them independently todeduce the ill-posedness and instability of the inverse algorithm.
Of all parameters, θs could be measured with relatively low uncer­
tainty. However, θs cannot be directly measured in soils with RF thatprevent the sampling of undisturbed cores. θs could be estimated as afunction of the dry bulk density and the mineral density of the solid
fraction when these properties are available. Alternatively, the easiest
way of estimating θs is to assume that its value equals the highest valueof θ measured in the soil profile after saturation, which implies ignoring
any entrapped air during saturation. An effective Ks could be measured
with field methods such as double ring and disc infiltrometers (Baetens
et al., 2009; Verbist et al., 2010). However, the effective Ks measured by
the infiltration experiment disregards the internal variability of Ks inheterogeneous soil profiles. Furthermore, the Ks as a model input often
requires model calibration due to its high sensitivity even whenmeasured values are available (Mertens et al., 2005; Rocha et al., 2006),
which implies little usefulness of measured values. Thus, soils with RF
and high heterogeneity will probably demand indirect determination of
all soil hydraulic parameters.
Differently from earlier studies that evaluated hydraulic proprieties
in disturbed samples created with a mix of RF and finer particles
collected from structured soils (Ma and Shao, 2008; Beibei et al., 2009;
Ma et al., 2010; Parajuli et al., 2017; Naseri et al., 2019), we successfully
optimized parameters of hydraulic functions by using inverse modelingof field drainage experiments, in which the real porous system is pre­served. Using θ and h data collected in the field (during internal
drainage experiments) and laboratory (using a dew point potentiom­eter), we determined that the combination of the Hydrus-1D inverse
solution in sequential mode with parameter estimation using RETC is asuitable strategy to determine hydraulic parameters for complex porous
media like soils with RF. These results highlight the suitability of inverse
modeling to indirectly determine soil hydraulic properties using tran­sient flow experiments (Šimůnek et al., 1998; Lazarovitch et al., 2007;
Jacques et al., 2012; Rashid et al., 2015; Pinheiro et al., 2019; Brunetti
et al., 2019). Furthermore, hydraulic parameter estimation by inverse
modeling in situ can also contribute to a better understanding of theeffects of RF on soil hydrological properties and processes (Zhang et al.,2016).

last step, the RMSWE (measured versus estimated water contents andpressure heads) was reduced to less than < 0.02 in most cases. Using
the optimized parameters of the unimodal van Genuchten-Mualem
function and lab data for pressure heads lower than − 200 cm, allother parameters of the bimodal van Genuchten function (w1, w2, α2,and n2) were estimated by a curve fitting procedure, in which thehigh incoherent values for θr1 were also corrected. The Hydrus-1D
inverse solution was then run again to optimize only Ks2 for thevan Genuchten-Mualem bimodal function with all previously esti­
mated parameters constant. In this way, 27 parameters were suc­cessfully estimated for layered soils with RF.
Declaration of Competing Interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.
Data availability
Data will be made available on request.
Acknowledgements
This study was financed in part by the Coordenação de Aperfeiçoa­
mento de Pessoal de Nível Superior - Brazil (CAPES) - Finance Code 001(via CAPES-PrInt/UFSM).
Appendix A. Supplementary data
Supplementary data to this article can be found online at https://doi.
org/10.1016/j.geoderma.2022.116240.
References
Abbasi, F., Šimůnek, J., Feyen, J., van Genuchten, M.T., Shouse, P.J., 2003. Simultaneous
inverse estimation of the soil hydraulic and solute transport parameters fromtransient field experiments: Homogeneous soil. Trans. ASAE 46 (4), 1085–1095.
Akaike, H., 1974. A new look at statistical model identification. IEEE Trans. Automat.
Contr. AC-19, 716–722.
Akaike, H., 1977. On entropy maximization principle. North-Holland, Amsterdam. In:
Krishnaiah, P.R. (Ed.), Applications of Statistics, pp. 27–41.
Baetens, J.M., Verbist, K., Cornells, W.M., Gabriels, D., Soto, G., 2009. On the influenceof coarse fragments on soil water retention. Water Resour. Res. 45, 1–14. https://doi.
org/10.1029/2008WR007402.
Beibei, Z., Ming’an, S., Hongbo, S., 2009. Effects of rock fragments on water movement
and solute transport in a Loess Plateau soil. Comptes Rendus - Geosci. 341, 462–472.
https://doi.org/10.1016/j.crte.2009.03.009.
Brunetti, G., Šimůnek, J., Bogena, H., Baatz, R., Huisman, J.A., Dahlke, H., Vereecken, H.,
2019. On the Information Content of Cosmic-Ray Neutron Data in the Inverse
Estimation of Soil Hydraulic Properties. Vadose Zo. J. 18, 1–24. https://doi.org/
10.2136/vzj2018.06.0123.
Dikinya, O., 2005. Comparison of the Instantaneous Profile Method and inverse
modelling for the prediction of effective soil hydraulic properties. Aust. J. Soil Res.
43, 599–606. https://doi.org/10.1071/SR04151.
Durner, W., 1994. Hydraulic conductivity estimation for soils with heterogeneous pore
structure. Water Resour. Res. https://doi.org/10.1029/93WR02676.
Gee, G.W., Bauder, J.W., 1986. Particle-size analysis. In: Kluter, A. (Ed.), Methods of Soil
Analysis: Part 1 - Physical and Mineralogical Methods, 2nd ed. Soil Science Society ofAmerica, Madison, pp. 383–411.
Gubiani, P.I., de Almeida, T.A., Mulazzani, R.P., Pedron, F. de A., Suzuki, L.E.A.S.,
Pereira, C.A., 2021a. Shaking settings to reduce the breakdown of entisol fragile
particles in texture analysis. Rev. Bras. Cienc. do Solo 45, 1–14. https://doi.org/
10.36783/18069657rbcs20210066.
Gubiani, P.I., Pereira, C.A., Cauduro, J.S., Campbell, C., Rivera, L., Pigatto, C.S.,
França, J.S., 2021b. Rock size fragments reduction allow including their effect onwater retention properties determined with a dew point potentiometer. Rev. Bras.
Cienc. do Solo 45, 1–11. https://doi.org/10.36783/18069657RBCS20200182.
Hillel, D., Krentos, V.D., Stylianou, Y., 1972. Procedure and test of an internal drainage
method for measuring soil hydraulic characteristics in situ. Soil Sci. 114, 395–400.
https://doi.org/10.1097/00010694-197211000-00011.
Hopmans, J.W., Šimůnek, J., Romano, N., Durner, W., 2018. Inverse Methods. Methods
Soil Anal. Part 4 Phys. Methods 963–1008. doi:10.2136/sssabookser5.4.c40.
Jacques, D., Smith, C., Šimůnek, J., Smiles, D., 2012. Inverse optimization of hydraulic,solute transport, and cation exchange parameters using HP1 and UCODE to simulate

5. Conclusions
This study aimed to determine suitable hydraulic functions for soils
with rock fragments (RF > 2 mm) and find a sequential optimization
strategy to estimate their parameters by inverse modeling using fielddrainage experiments. Four layered soil profiles with different amountsof RF in different layers (from zero to 99 %) were evaluated in the fieldusing internal drainage experiments and dewpoint potentiometers in thelab. The two main conclusions are:1- Bimodal functions should be considered to represent hydraulicproperties (water retention and conductivity) of soils with RF, i.e., insoils with large pores and cracks surrounding RF, which contrast
with textural pores in the aggregated fine material and fine pores
inside porous RF. A disadvantage of these functions is an increase in
the number of calibrated parameters, leading to ill-posed problems
and instability of the optimization algorithm.
2- It is possible to sequentially optimize the van Genuchten parametersof layered soil profiles and successfully avoid failure of the Hydrus1D inverse solution due to non-convergence. In this study, the pa­
rameters of the unimodal van Genuchten-Mualem function(subscript 1) were optimized using the Hydrus-1D inverse solution
repeatedly (steps). α1 and Ks1 were optimized in step 1, θs1 and θr1 instep 2, α1, n1, and Ks1 in step 3, and θs1 and θr1 again in step 4. At the10

P. Ivonir Gubiani et al.

Geoderma 429 (2023) 116240
Romano, N., 1993. Use of an inverse method and geostatistics to estimate soil hydraulicconductivity for spatial variability analysis. Geoderma 60, 169–186. https://doi.org/
10.1016/0016-7061(93)90025-G.
Schoeneberger, P.J., Wysocki, D.A., Benham, E.C, Soil Survey Staff, 2012. Field book fordescribing and sampling soils, Version 3.0. Natural Resources Conservation Service,
National Soil Survey Center, Lincoln, NE, 3. https://www.nrcs.usda.gov/sites/default/files/2022-09/field-book.pdf.
Šimůnek, J., and J. W. Hopmans, 2002. Parameter Optimization and Nonlinear Fitting,in: J. H. Dane and G. C. Topp (eds.) Methods of Soil Analysis, Part 4, Physical
Methods, Chapter 1.7, Third edition, SSSA, Madison, WI, 139-157.
Šimůnek, J., Angulo-Jaramillo, R., Schaap, M.G., Vandervaere, J.P., Van Genuchten, M.
T., 1998. Using an inverse method to estimate the hydraulic properties of crusted
soils from tension-disc infiltrometer data. Geoderma 86, 61–81. https://doi.org/
10.1016/S0016-7061(98)00035-4.
Šimůnek, J., Šejna, M., Saito, H., Sakai, M., van Genuchten, MTh., 2013. The HYDRUS1D software package for simulating the movement of water, heat, and multiple
solutes in variably saturated media. Version 4.17. Riverside: University of California
Riverside; 2013. (HYDRUS Software Series 3).
Šimůnek, J., Hopmans, J.W., 2002a. Chapter 1.7: Parameter optimization and nonlinear
fitting. In: Dane, J.H., Topp, G.C. (Eds.), Methods of Soil Analysis: Part 1. Physical
Methods, 3rd ed. SSSA, Madison, Wisc, pp. 139–157.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2012. Hydrus: model use, calibration, andvalidation. Trans. Asabe 55, 1261–1274.
Soil Survey Staff, 2014. Keys to soil taxonomy. 12th ed. Washington, DC: United States
Department of Agriculture, Natural Resources Conservation Service; 2014.
van Genuchten, M.T., 1980. A Closed-form Equation for Predicting the Hydraulic
Conductivity of Unsaturated Soils. Soil Sci. Soc. Am. J. 44, 892–898. https://doi.org/
10.2136/sssaj1980.03615995004400050002x.
van Genuchten, M. Th., F. J. Leij, and S. R. Yates. 1991. The RETC code for quantifying
the hydraulic functions of unsaturated soils. Report No. EPA/600/2-91/065. Ada,
Okla.: U.S. Environmental Protection Agency, R.S. Kerr Environmental ResearchLaboratory.
Verbist, K., Torfs, S., Cornelis, W.M., Oyarzún, R., Soto, G., Gabriels, D., 2010.
Comparison of Single- and Double-Ring Infiltrometer Methods on Stony Soils.
Vadose Zo. J. 9, 462–475. https://doi.org/10.2136/vzj2009.0058.
Weiler, M., Flühler, H., 2004. Inferring flow types from dye patterns in macroporous
soils. Geoderma 120, 137–153. https://doi.org/10.1016/j.geoderma.2003.08.014.
Zhang, Y., Zhang, M., Niu, J., Li, H., Xiao, R., Zheng, H., Bech, J., 2016. Rock fragments
and soil hydrological processes: Significance and progress. Catena 147, 153–166.
https://doi.org/10.1016/j.catena.2016.07.012.
Zijlstra, J., Dane, J.H., 1996. Identification of hydraulic parameters in layered soils basedon a quasi-Newton method. J. Hydrol. 181 (1–4), 233–250. https://doi.org/
10.1016/0022-1694(95)02909-5.

cation exchange. J. Contam. Hydrol. 142–143, 109–125. https://doi.org/10.1016/j.
jconhyd.2012.03.008.
Lazarovitch, N., Ben-Gal, A., Šimůnek, J., Shani, U., 2007. Uniqueness of Soil HydraulicParameters Determined by a Combined Wooding Inverse Approach. Soil Sci. Soc.
Am. J. 71, 860–865. https://doi.org/10.2136/sssaj2005.0420.
Ma, D., Shao, M., 2008. Simulating infiltration into stony soils with a dual-porosity
model. Eur. J. Soil Sci. 59, 950–959. https://doi.org/10.1111/j.13652389.2008.01055.x.
Ma, D.H., Shao, M.A., Zhang, J.B., Wang, Q.J., 2010. Validation of an analytical methodfor determining soil hydraulic properties of stony soils using experimental data.
Geoderma 159, 262–269. https://doi.org/10.1016/j.geoderma.2010.08.001.
Mavimbela, S.S.W., Van Rensburg, L.D., 2013. Estimating hydraulic conductivity ofinternal drainage for layered soils in situ. Hydrol. Earth Syst. Sci. 17, 4349–4366.
https://doi.org/10.5194/hess-17-4349-2013.
Mertens, J., Madsen, H., Kristensen, M., Jacques, D., Feyen, J., 2005. Sensitivity of soil
parameters in unsaturated zone modelling and the relation between effective,
laboratory and in situ estimates. Hydrol. Process. 19, 1611–1633. https://doi.org/
10.1002/hyp.5591.
Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated
porous media. Water Resour. Res. 12, 513–522.
Naseri, M., Iden, S.C., Richter, N., Durner, W., 2019. Influence of Stone Content on Soil
Hydraulic Properties: Experimental Investigation and Test of Existing Model
Concepts. Vadose Zo. J. 18, 1–10. https://doi.org/10.2136/vzj2018.08.0163.
Parajuli, K., Sadeghi, M., Jones, S.B., 2017. A binary mixing model for characterizing
stony-soil water retention. Agric. For. Meteorol. 244–245, 1–8. https://doi.org/
10.1016/j.agrformet.2017.05.013.
Pinheiro, E.A.R., de Jong van Lier, Q., Inforsato, L., Šimůnek, J., 2019. Measuring fullrange soil hydraulic properties for the prediction of crop water availability using
gamma-ray attenuation and inverse modeling. Agric. Water Manag. 216, 294–305.
https://doi.org/10.1016/j.agwat.2019.01.029.
Rashid, N.S.A., Askari, M., Tanaka, T., Šimůnek, J., van Genuchten, M.T., 2015. Inverse
estimation of soil hydraulic properties under oil palm trees. Geoderma 241–242,
306–312. https://doi.org/10.1016/j.geoderma.2014.12.003.
Rawls, W.J., Brakensiek, C.L., Saxton, K.E., 1982. Estimation of soil water properties.
Trans. - Am. Soc. Agric. Eng. 25 https://doi.org/10.13031/2013.33720.
Ritter, A., Hupet, F., Muñoz-Carpena, R., Lambot, S., Vanclooster, M., 2003. Using
inverse methods for estimating soil hydraulic properties from field data as analternative to direct methods. Agric. Water Manag. 59, 77–96. https://doi.org/
10.1016/S0378-3774(02)00160-9.
Rocha, D., Abbasi, F., Feyen, J., 2006. Sensitivity Analysis of Soil Hydraulic Properties on
Subsurface Water Flow in Furrows. J. Irrig. Drain. Eng. 132, 418–424. https://doi.
org/10.1061/(asce)0733-9437(2006)132:4(418).

11