Navin K. C. Twarakavi*
Dep. of Environmental Sciences
Univ. of California
Riverside, CA 92521

Hirotaka Saito
Tokyo University of Agriculture and Technology
Dep. of Ecoregion Science
Fuchu, TokyoJapan 183-8509

Jirka Šimunek
Dep. of Environmental Sciences
Univ. of California
Riverside, CA 92521

M. Th. van Genuchten
U.S. Salinity Laboratory
450 W. Big Springs Rd.
Riverside CA 92507

SOIL PHYSICSA New Approach to Estimate Soil HydraulicParameters Using Only Soil Water Retention Data
Traditional approaches for estimating soil hydraulic parameters (such as the RETC code)
perform well when experimental data for both the retention curve and hydraulic conductivity function are available; however, unsaturated hydraulic conductivity data are often
unavailable. The objective of this work was to develop an approach to estimate robust soil
hydraulic parameters from water retention curve data alone. The proposed approach, called
the Multiobjective Retention Curve Estimator (MORE), is based on the Multiobjective
Shuffled Complex Evolution Metropolis (MOSCEM-UA) algorithm and estimates an optimal parameter set by simultaneously minimizing two objective functions each representing
water content and relative unsaturated hydraulic conductivity residuals. To address the lackof observed unsaturated hydraulic conductivity data, MORE transforms both predicted andobserved water contents into the hydraulic conductivity space using a pore-size distribution
model (e.g., the Mualem model) and also optimizes the soil hydraulic parameters in thisfictitious space. We applied MORE to two cases. In the first case, MORE was used to estimate soil parameters using only retention curve data for 12 random soils selected from the
UNSODA database. While the soil hydraulic parameters estimated using RETC and MORE
fit the retention curve similarly, the MORE approach consistently decreased the error in
fitting unsaturated hydraulic conductivities by as much as 5% compared with RETC. The
second case involved using the parameters fitted using the MORE and RETC approaches tomodel a field-scale experiment. Compared with RETC, the error in predicted water contents
decreased by 25% using parameters predicted by MORE. The MORE approach was shownto fit robust soil hydraulic parameters; however, the approach is relatively slower and moretime consuming than RETC.
Abbreviations: MORE, Multiobjective Retention Curve Estimator; MOSCEM-UA, Multiobjective
Shuffled Complex Evolution Metropolis.

T

he variably saturated flow process in soils is a highly nonlinear and dynamic phenomenon. Commonly used numerical models for variably saturated water flow use the classical
Richards equation (Richards, 1931). The Richards equation ishighly nonlinear due to the dependence of unsaturated hydraulic conductivity and water content on the capillary pressure
head. Consequently, modeling variably saturated water flow
involves solving the Richards equation along with a nonlinear soil hydraulic model (Gardner, 1958; Brooks and Corey,
1964; van Genuchten, 1980; Leij et al., 1997) that characterizes the relationships between the water content, capillarypressure head, and unsaturated hydraulic conductivity. Of allavailable soil hydraulic models, the van Genuchten–Mualem
model (van Genuchten, 1980) is perhaps the most widely usedmodel in simulation of vadose zone flow processes and, thereSoil Sci. Soc. Am. J. 72:471–479
doi:10.2136/sssaj2006.0342
Received 3 Oct. 2006.
*Corresponding author (navink@ucr.edu).
© Soil Science Society of America
677 S. Segoe Rd. Madison WI 53711 USA
All rights reserved. No part of this periodical may be reproduced ortransmitted in any form or by any means, electronic or mechanical,
including photocopying, recording, or any information storage andretrieval system, without permission in writing from the publisher.
Permission for printing and for reprinting the material contained
herein has been obtained by the publisher.
SSSAJ: Volume 72: Number 2 • March–April 2008

fore, used in this study for representing the complex relationships between the hydraulic conductivity, the water content,
and the pressure head.
Due to the multidimensionality of the parameter space, acommon dilemma during parameter estimation in hydrologic
processes is that there are a number of parameter sets that can fit
the data equally well. The ultimate goal during the soil hydraulic parameter estimation is to obtain parameter estimates thatare the most representative of flow processes in soils. An application of the van Genuchten–Mualem model in vadose zone
flow modeling requires estimates of the following soil hydraulicparameters: the so-called shape parameters (α and n), saturated
hydraulic conductivity (Ks), saturated water content (θs), andresidual water content (θr). Of the required parameter set, Ks
and θs can be estimated experimentally. The remaining parameters (α, n, and θr) have to be estimated indirectly by fitting thevan Genuchten–Mualem model to the experimental retention
and unsaturated hydraulic conductivity data using an optimization approach, such as the RETC code (van Genuchten et al.,1991). It has been traditionally assumed that a good estimate of
the necessary soil hydraulic parameters is achieved when the estimated parameter set can accurately represent both the retention
curve and hydraulic conductivity–pressure head relationships aswell as the transient flow processes.
The RETC code is a widely used computer program for estimating the soil hydraulic parameters that represent the retention
curve and the unsaturated hydraulic conductivity–pressure head471

relationships (van Genuchten et al., 1991). The RETC code uses
a nonlinear least-squares optimization approach to estimate theunknown model parameters from observed retention, conductivity, or diffusivity data. The aim of this curve-fitting process
is to find the set of parameters that best minimizes the sum ofsquared residuals, SSQ, between model predictions and data (vanGenuchten et al., 1991). The SSQ reflects the degree of model
bias (lack of fit) and the contribution of random errors. Note that
the SSQ objective function is statistically synonymous with the
RMSE. The RETC code minimizes the SSQ iteratively by meansof a weighted least-squares approach based on Marquardt’s maximum likelihood method (Marquardt, 1963).
Previous studies have shown that RETC can be efficiently
used to fit the retention curve and the unsaturated hydraulicconductivity function when experimental data are available forboth properties (e.g., van Genuchten et al., 1991; Yates et al.,1992). While experimental data for the retention curve can becollected relatively easily, accurate measurements of unsaturated hydraulic conductivities are generally difficult and time
consuming (Schaap and van Genuchten, 2005). Soil physicists,
however, often find themselves in the predicament where theymust estimate the soil hydraulic parameters that adequately
represent the entire gamut of pressure head–water content–
hydraulic conductivity relationships using only the retention
curve data. When only retention curve data are available,
RETC sometimes fits the retention curve at the cost of poorer
characterization of the unsaturated hydraulic conductivity. For
example, the van Genuchten–Mualem model, when fitted only
to the experimental retention curve data, often results in a good
fit of the retention curve but does not necessarily ensure a good
description of the hydraulic conductivity–pressure head relationship (see below for illustrations). This is mainly due to thevarying sensitivities of the retention curve and the unsaturated
hydraulic conductivity function to the parameter set.

Fig. 1. Schematic of the location of the Pareto optimal solution of amultiobjective problem in the objective space. Circles represent
various nonoptimal parameter sets that do not provide the bestsolution. Points A and B represent parameter sets that provide
the best solution for individual objectives F1 and F2, respectively.
The line connecting Points A and B is the nondominated “Pareto
optimal” set of solutions. Point C corresponds to the Pareto solution that may be regarded as the “best” possible solution.
472

In the past, many improvements have been suggested to
the RETC code (e.g., Hollenbeck et al., 2000; Abbaspour et al.,2001); however, these improvements do not address conditions
where only the soil water retention data are available. Similar
to the RETC approach, these improved methods often resultin parameter sets that inadequately represent the unsaturated
hydraulic conductivity–pressure head relationship.
Traditional methods, such as RETC, for estimating thesoil hydraulic parameters are formulated as a single-objective
optimization problem. The classical single-objective optimization approach is based on the central assumption that thechosen single-objective function can adequately and correctly
represent all errors that may be incurred due to the parameterset used in the model (Vrugt et al., 2003). Gupta et al. (1999)
have convincingly pointed out that this may not be true, andthat the single-objective function can sometimes contributestructural errors during model calibration that may even exceed
measurement errors.
In this study, we developed a multiobjective optimization approach, called MORE, for estimating the soil hydraulicparameters that is designed to overcome the aforementioned
weaknesses of the single-objective optimization in RETC.
Multiobjective optimization is not an alien concept for soil
physicists and has been explored in many previous studies (e.g.,
Neuman, 1973; Vrugt and Dane, 2005; Schoups et al., 2005;
Mertens et al., 2006). The basic idea behind multiobjective
optimization methods is finding parameter sets that simultaneously minimize multiple objective functions. The solution tothis problem, however, will generally no longer be a single “best”
parameter set but will consist of a group of possible solutions
due to tradeoffs among different objective functions (Vrugt etal., 2003). A key characteristic of this so-called “Pareto set ofsolutions” is that each parameter set within the Pareto set of solutions is better than the others for at least one objective function,
but none is better than the others for all objective functions.
Figure 1 gives an illustration of this concept. The Pareto
set of solutions would also include the solutions for the individual objectives. For example, while Point A represents parameters that provide the best solution for the objective F1, Point B
represents the optimal parameter set for the objective F2. Manyalgorithms are available for estimating the “optimal” parameter
sets through multiobjective optimization, such as MOCOMUA (Yapo et al., 1998) and MOSCEM-UA (Vrugt et al.,2003). One may select a “best” possible solution (represented
by Point C) from the Pareto set of solutions by minimizing thenormalized Euclidean measures of distance of the Pareto set
to the origin (Vrugt et al., 2003). It is possible to estimate thesame “best” possible solution using a weighted single-objectivefunction; however, the relative values of weights would dominate the final “best” possible solution. It has been convincingly
argued that the “best” solution has the smallest normalized distance. Throughout this study, the point on the Pareto set withthe minimal normalized Euclidean distance to the origin has
been selected as the “best” possible solution.

MATERIALS AND METHODS
Proposed Approach
For a one-dimensional scenario, the Richards equation is describedmathematically as follows:
SSSAJ: Volume 72: Number 2 • March–April 2008⎤∂θ ∂ ⎡∂h= ⎢ K (h ) − K (h )⎥ −S∂t ∂z ⎣∂z⎦[1]

where θ is the volumetric water content, h is the soil water pressure
head (L), t is time (T), z is the distance from a reference datum (L),K(h) is the unsaturated hydraulic conductivity (L T−1) as a function ofh or θ, and S represents flow from sinks and sources.
The van Genuchten–Mualem model, which is used to represent
the functions K(h) and θ(h) in the flow models, is described asθs − θ r
⎧
⎪θr +
n m
1+ α hθ( h ) = ⎨
⎪
⎩ θs

(

Se (h ) =

)

h <0

[2a]

h ≥0

θ( h ) − θ r
θs − θ r[2b]m
K ( h ) = K s S e l ⎡ 1− (1− S e 1/ m ) ⎤⎣⎦m = 1−1 n

2

[2c]
[2d]

n >1

where θ(h) is the volumetric water content at pressure head h, θs and θr
are the saturated and residual volumetric water contents, respectively;Se(h) is the scaled water content at pressure head h (L), Ks is the saturated hydraulic conductivity (L T−1), α (L−1) and n are van Genuchten’s
shape parameters, and l is a tortuosity or pore-connectivity parameter
estimated by Mualem (1976) to be 0.5. Equations [1] and [2] havebeen used with or without additional modifications to model variably
saturated water flow and have been successfully extended into two andthree dimensions (e.g., Šimunek et al., 1999; Preuss et al., 2004).
The MORE approach utilizes two objective functions to obtain asimultaneous fit of the retention curve and the unsaturated hydraulicconductivity function. The first objective function reflects the lack of fitof observed water contents, while the second objective function assesses
the lack of fit of observed unsaturated hydraulic conductivities. The
goal of the MORE methodology is to find the parameter set such thatboth objective functions, O(b), described in Eq. [3] are minimized:⎧ RMSE ⎡ θ i , θˆ i ( b) ⎤⎣⎦⎪O( b) = ⎨⎪ RMSE log ⎡⎣ ki , (θ i , b )⎤⎦ , log ⎡ kˆi θˆ i , b ⎤
⎣
⎦
⎩

{

( )}

[3]

where b [ = α, n, θs, and θr] is the set of optimized soil hydraulicparameters, θi are the observed water contents, θˆ i ( b) are the fitted
water contents as a function of the parameter set b, kˆi ( θˆ i , b) are thefitted relative hydraulic conductivities as a function of the parameterset b and fitted water contents θˆ i ( b) , and ki(θi,b) are the apparent
observed relative hydraulic conductivities as a function of the parameter set b and observed water contents θi. In the absence of observed
relative conductivities, ki(θi), apparent observed relative hydraulic conductivities, ki(θi,b), are estimated using the parameter set b
and observed water contents, θi, from Eq. [2a–2d]. It is assumed
that ki(θi,b) closely reflects the actual relative hydraulic conductivity,ki(θi), that could have been estimated experimentally. Similarly, fitted relative hydraulic conductivities, kˆi ( θˆ i , b) , are estimated using fitted water contents, θˆ i ( b) , and the parameter set b. The optimization
based on Eq. [3] thus simultaneously minimizes water content residuSSSAJ: Volume 72: Number 2 • March–April 2008

als in both the water content and the relative unsaturated hydraulicconductivity objective function spaces. A pore-size distribution model(e.g., Mualem, 1976) is used to project both fitted and observed watercontents into the hydraulic conductivity space. Similar to the RETC
approach, the RMSE is used in MORE as the objective function. Notethat the RMSE objective function for the relative unsaturated hydraulicconductivities considers log-transformed values as unsaturated hydraulic conductivities that vary across many orders of magnitude.
Several assumptions are made in the MORE approach. Note thatthe water content data are projected into both the retention curve andhydraulic conductivity objective spaces in MORE. The second major
assumption is that pore-size distribution models, such as the Mualem
model, accurately transform water contents into relative hydraulicconductivities. This assumption is especially crucial, as the accuracyof the second objective function is entirely dependent on the validityof the pore-size distribution model for the soil under consideration.
The MOSCEM-UA algorithm is used to solve the multiobjective
optimization problem posed in MORE. The algorithm in the MORE
approach is essentially the same as the MOSCEM-UA algorithm (fora detailed description of the approach, see Vrugt et al., 2003). The
only additions in the MORE approach are in the multiobjective vector estimation stage described above. Figure 2 shows the flow chart forcalculating the multiobjective vector used in the MORE approach fora given parameter set and experimental retention data when the van
Genuchten–Mualem model is considered.
Only a brief explanation of the procedure used in the MOSCEMUA algorithm is presented here. The MOSCEM-UA algorithm starts
with an initial population of points from the feasible parameter space
defined by the user. For each individual of the population, values ofthe multiple objective functions are computed and the population isranked and sorted based on their fitness. The fitness of each individual

Fig. 2. A flowchart for calculating the multiobjective vector for agiven parameter set in the Multiobjective Retention Curve Estimator (MORE) approach.

473

0.76
0.75

RMSE(log(k))

0.74
0.73
0.72
0.71
0.7
0.69
0.68
0.67

0.02

0.021 0.022 0.023 0.024 0.025 0.026 0.027

RMSE( θ)Fig. 3. Tradeoffs in the objective space of the Multiobjective Retention
Curve Estimator (MORE) algorithm for the silt soil (UNSODA code
no. 4670). Gray points represent different parameter sets considered by the MORE approach. Circles represent the final nondominated Pareto set after 10,000 function evaluations. The star symbol represents the parameter set estimated using RETC and thesolid square symbol corresponds to the “best” possible parameterset estimated using the MORE approach; k is the relative hydraulicconductivity and θ is the volumetric water content.

is based on how much the multiple objectives are minimized. The
population (100 individuals were used in this study) is then partitioned into complexes (five in this study) and in each complex, parallel sequences of individuals are launched and new candidates arederived from the structure of the sequences and associated complexes.
A Metropolis–Hastings algorithm is used to decide if the candidate
points are acceptable and in case of acceptance, the worst memberof the current complex is replaced. The procedure is continued until
a prescribed number of iterations is reached (here 10,000), after
which they are combined and new complexes are formed through theprocess of shuffling. This process converges to the Pareto set. Fromthe estimated Pareto set, the nondominated Pareto set is selected. A
Pareto set of parameters is said to be nondominated if no parameterset exists that has the lowest values of all objective functions. The nondominated Pareto set is illustrated by a thick line in Fig. 1. From thenondominated Pareto set, the parameter set superior to other possible
solutions is chosen as the “best” possible parameter set (represented
by Point C in Fig. 1). The MORE approach was implemented in the
MATLAB environment (The MathWorks, 1999) and was primarily
based on the MOSCEM-UA algorithm (Vrugt et al., 2003). Figures 3
and 4 are practical examples of Fig. 1 and further illustrate the parameter sets obtained by the MOSCEM-UA algorithm.
The performance of MORE was evaluated using two case studies of varying complexity. The first case study compared the performance of RETC and MORE in estimating parameters for a variety ofsoils of different textures and structures randomly selected from the
UNSODA database. The second case study involved estimating thesoil hydraulic parameter set from experimental retention curve data
for the Las Cruces one-dimensional variably saturated flow experiment (Wierenga et al., 1991). While the first case study involved
laboratory measurements of core samples with smaller heterogeneityand fewer experimental data, the second case study involved a fieldscale experiment with larger heterogeneity and more experimental
data. HYDRUS-1D (Šimunek et al., 2005) was used to model the
one-dimensional variably saturated water flow with parameters estimated from RETC and MORE in the second case study. To compare
the observed data with those predicted using the MORE and RETC
approaches, we used the RMSE statistic. The RMSE is a commonly usederror statistic that is the square root of the mean of the squared residuals.

Case 1: Selected Soils from the UNSODA Database

Relative hydraulic conductivity (−)

Water content (−)

In comparing the performance of the MORE and RETC
approaches for estimating parameters from soil water retention experimental data, we randomly chose 12 soils of different textures and structures from the UNSODA database (Leij et al., 1996; Nemes et al., 2001).
These selected soils were chosen such that retention curve and unsaturated hydraulic conductivity data are available. Among these soils in the
UNSODA database, many have similar retention curve and unsaturated
hydraulic conductivity data. Most of these soils with similar retention
curves and unsaturated hydraulic conductivity data shared the same sampling location. Care was taken to avoid such repetitions. Also, soils withobvious sampling and experimental errors were ignored. The selectedsoils, covering a wide range of textures and struc1
0.5tures, are listed in Table 1. Although the MORE
0.9
0.45
(b)
and RETC approaches were compared for their(a)0.8
performance for all 12 selected soils, only one of
0.4
0.7
these randomly selected soils is used here to illus0.35
0.6
trate in more detail the MORE approach.
0.3
Even though the selected soils had experi0.5
0.25
mentaldata for both soil water retention and
0.4
unsaturated
hydraulic conductivity, only the
0.2
0.3
retention curve data were used to estimate
0.15
0.2
the soil hydraulic parameters. The estimated0.1
0.1
parameter set was then used to see how well0
theMORE and RETC approaches were able
0.05
0
5
0
5
10
10
10
10
to characterize the unsaturated hydraulicPressure head (cm)Pressure head (cm)conductivity by comparing the experimental
Fig. 4. Fitted (a) water content–pressure head and (b) relative unsaturated hydraulic conductivunsaturated hydraulic conductivity data with
ity–pressure head relationships obtained using parameter estimates from the Multiobjective
Retention Curve Estimator (MORE) (solid line) for the silt soil (UNSODA code no. 4670). predicted values. In other words, the experiThe points represent the experimental data. The dotted line is the solution based on RETC. mental unsaturated hydraulic conductivities
Gray lines represent the several parameter sets corresponding to the nondominated Pareto
set. Only retention data were used in parameter optimization by both RETC and MORE.

474

SSSAJ: Volume 72: Number 2 • March–April 2008

were used only to evaluate the robustness of the “best” parameter sets
estimated using MORE and RETC.

Case 2: Las Cruces Trench Experimentsoil hydraulic parameters from only the retention data. The
optimized retention parameters were as follows: θs = 0.441, θr= 0.000, α = 0.005 cm−1, and n = 1.41.
The MORE approach was then used to obtain the soil
hydraulic parameters using multiobjective optimization. Figure
3 shows the evolution of the “best” parameter set in the objective space using the MORE approach. Figure 3 is similar to
Fig. 1, which illustrates the concept of the Pareto optimality.
Figure 3 also shows the location of the RETC solution andthe “best” possible solution obtained by MORE in the objective space. From the nondominated Pareto set of solutions, the“best” possible parameter set (illustrated by Point C in Fig. 1)
was chosen as follows: θs = 0.429, θr = 0.030, α = 0.004 cm−1,and n = 1.56.
The estimated parameter sets obtained by RETC and
MORE were then used to predict the retention curves andhydraulic conductivity–pressure head functions using Eq. [2].
Figure 4 compares observed water contents and relative hydraulic conductivities with the corresponding soil hydraulic functions obtained from MORE and RETC. It is clear that RETC
fits water contents much better than it predicts relative hydraulic conductivities. Relative unsaturated hydraulic conductivities
are clearly underestimated when parameters estimated only
from the retention curve are used. An underprediction of relative hydraulic conductivities would translate into a predictionof flow rates slower than in reality during model simulations.
On the other hand, the MORE approach shows a significantly
improved correspondence between measured and predicted
relative unsaturated hydraulic conductivities (Fig. 4b), despite
the fact that measured hydraulic conductivities were not usedin the optimization process. Although the MORE approach
predicts the hydraulic conductivity function significantly better than the RETC approach, the retention curves fitted by
the two approaches differ only slightly (Fig. 4a). Figures 5aand 5b show contour plots of the RMSE between the observed
and estimated water contents and relative hydraulic conductivities, respectively, for different combinations of the α andn parameters. Figure 5 also shows the location of the RETC
and MORE solutions. One can immediately observe the dissimilarity of the two error surfaces. For example, the hydraulic

The second case study involved modeling of the one-dimensional
infiltration for the Las Cruces trench site experiment (Wierenga et al.,
1989, 1991). The Las Cruces trench site experiment was a comprehensive field study conducted in southern New Mexico. The primary
purpose of the study was to develop a data set for validating and testing numerical models. For this purpose, the study site was heavily
instrumented with tensiometers, neutron probes, and solute samplers
for collecting water contents, pressure heads, and solute concentrations. More than 500 soil samples (undisturbed and disturbed) weretaken at the experimental site and analyzed in the laboratory for bulk
density and to find the saturated hydraulic conductivity and soil water
retention curve. The experimental infiltration study involved application of water to a 4-m-wide area using closely spaced drips with anaverage surface flux of 1.82 cm d−1 for the 86 d of the experiment.
During the experiment, no surface ponding was observed. To reducethe disruption of the experimental conditions by rain and evaporation, the irrigated area and the surroundings were covered by a pond
liner. To evaluate the data set, Wierenga et al. (1991) developed asimple one-dimensional numerical model assuming uniform soil
conditions. Similar to HYDRUS-1D, the model was based on the
Richards equation (Eq. [1]) for one-dimensional variably saturated
water flow. For modeling purposes, Wierenga et al. (1991) assumed
the equivalent saturated hydraulic conductivity (Ks) as a geometricmean of laboratory-measured values from all soil samples. The retention curve parameters were estimated from the laboratory data using
an approach similar to RETC (Wierenga et al., 1991). Water contents
were observed along the depth on Days 19 and 35 after the experimented was started.
Similar to the work done by Wierenga et al. (1991), the onedimensional simulation of the infiltration in the Las Cruces experiment was performed with the MORE-estimated soil hydraulic parameters. Initial pressure heads (hi = −100 cm) in the soil profile were
the same as those used in Wierenga et al. (1991). While a constantwater flux was used as the upper boundary condition (q0 = 1.82 cmd−1), free drainage was assumed at the lower boundary. Also, similarboundary conditions and equivalent saturated hydraulic conductivity(Ks = 270.1 cm d−1) were used. The retention
Table 1. Summary of 12 soils selected from the UNSODA database for the first case study.
curve parameters were estimated using again
Note that the silt soil (no. 4670, in italics) was used for a thorough analysis of thethe RETC and MORE approaches. Except forMultiobjective Retention Curve Estimator (MORE) approach.
the retention curve parameters, all remaining
Experimental data points
parameters and initial and boundary condi- UNSODA
TextureStructure
Soil water Unsaturated hydraulicno.
tions were kept the same. This allowed us toretention
conductivity
compare the effectiveness of the RETC and——————— no. ———————MORE approaches in estimating parameter
4672
silt loamfine to moderate polyhedric25
25
sets for variably saturated flow models.

RESULTS AND DISCUSSION
Case 1: Selected Soils from the
UNSODA Database
One selected soil (undisturbed coresample of silt collected from Hannover,Germany; UNSODA code no. 4670) was
initially used to analyze and demonstrate
the MORE approach. For the selected soil,
the RETC code was used to estimate the

2242
1400
4011
1382
4661
4670
1490
1361
1371
1270
2740

SSSAJ: Volume 72: Number 2 • March–April 2008

sandy loam
granularclay
angularloamy sandgranularsandy clay loam angularsand
single grain
silt
blockysilt loam
blockysilty clay
crumb angularsilty clay loam angular–slightly weldedclay loam
angular

11
28
24
46
25
25
39
45
48
5
6

12
10
11
10
25
25
10
10
10
8
13

475

n (−)

n (−)

3
0.1

0.
10

09

3
04
0.

0.

02

1.5

(b)

7

1.4
0.03

1.2

1.5
0.13

1.3

0.07

0.23

1.2

1.1

0.34

1.1
2

4

6

α (1/cm)

0

1.6
1.4

1.3

.07

1.7

05

1.6

1.8

07
0.

0.

1.7

0.27

0.

1.8

0

1.9

0.20

(a)

1.9

8

0
x 10

−3

2

4

6

α (1/cm)

8
−3

x 10

Fig. 5. Contour plots of the RMSE between the observed and fitted (a)water contents and (b) log of relative hydraulic conductivities as afunction of shape parameters α (cm−1) and n for the silt soil (UNSODA code no. 4670). The star symbol represents the parameter set
estimated using RETC and the solid square symbol corresponds tothe “best” possible parameter set estimated using the Multiobjective Retention Curve Estimator (MORE) approach. Gray points are
the final nondominated Pareto set estimated by the MORE approach.
Only retention data were used in parameter optimization by both
RETC and MORE.

Fig. 6. Tradeoffs in the objective space of the Multiobjective Retention
Curve Estimator (MORE) algorithm for selected soils. Gray points
represent different parameter sets considered by the MORE approach. Circles represent the final nondominated Pareto set after 10,000 function evaluations. The star symbol represents theparameter set estimated using RETC and the solid square symbol
corresponds to the “best” possible parameter set estimated using
the MORE approach; k is the relative hydraulic conductivity andθ is the volumetric water content.
476

conductivity is relatively more sensitive to the n parameter than
to the α parameter. Figure 5 further illustrates that retention
and hydraulic conductivity functions have different minima.
While the estimated RETC parameter set is located well inside
the global minimum for water contents, it is not necessarily inthe region of the global minimum for relative hydraulic conductivities. Similar observations have been reported in previous
studies (e.g., Yates et al., 1992). Yates et al. (1992) analyzed the
RETC code for a range of soils from the UNSODA databaseand noted that estimates of the unsaturated hydraulic conductivity were, on average, less accurate. Because of the often poor
characterization of the unsaturated hydraulic conductivity, one
may expect that RETC may lead to parameter estimates thatare less representative of the flow characteristics of the medium.
Compared with the RETC solution, the MORE solution isbetter placed in the region of global minima for both the watercontent and relative hydraulic conductivity. We emphasize here
again that no measured hydraulic conductivity data were usedin optimizing the soil hydraulic parameters and that these measured data were used only to evaluate the solutions obtained by
the two methods.
It should be noted that the MOSCEM-UA algorithm on
which the MORE approach is based has a set of user-input
values. These include the number of complexes, the numberof function evaluations, and the size of the population. In our
analysis of the various soils selected from the UNSODA database, approximately 10,000 function evaluations were sufficient
to attain a robust Pareto set of solutions. The sensitivity of the
MORE code to the number of complexes and population size
was analyzed and is reported in Table 2. Similar to the observation made by Vrugt et al. (2003), the number of Pareto points
generated by the MORE approach is relatively insensitive tothe population size and the number of complexes. The “best”
possible solution was also relatively insensitive to the population
size and number of complexes. Similar observations were madeduring the sensitivity analysis of the remaining selected soils fromthe UNSODA database.
The MORE approach seemed to consistently provide better
results than the RETC code for a variety of soils. Figure 6 shows
the tradeoff in the objective space of the MORE approach for several selected soils. One may note that the Pareto set is relatively uniformly distributed. Any gaps that may be observed in the Pareto
set distribution (Soil no. 1490, for example) are due to poor representation of the data by the model. Table 3 lists the “best” parameter sets estimated by the MORE and RETC approaches, while
Table 4 lists the RMSE estimates between fitted and observed
water contents and predicted and observed unsaturated hydraulic
Table 2. The total number of nondominated PARETO points after
10,000 function evaluations with the Multiobjective Retention Curve Estimator (MORE) approach for the silt soil (UNSODA no. 4670) as a function of population size 10, 20, 50,
or 100 and the number of complexes.
Complexes

1
2
5
10

Nondominated Pareto points
20
50
100
––––––––––––––no.––––––––––––––
84
148
114
83
153
110
88
108
98
100

10

SSSAJ: Volume 72: Number 2 • March–April 2008

conductivities. Even though unsaturated hydraulic Table 3. The final parameter sets (residual and saturated volumetric water contents,
θr and θs, respectively, and shape parameters α and n) estimated by the RETC
conductivity data were not used in either the RETC
and Multiobjective Retention Curve Estimator (MORE) approaches for the ranor the MORE approaches, the parameter sets gendomly selected soils from the UNSODA database.
erated by MORE showed an improved characterRETC
MORE
UNSODAization of the hydraulic conductivity–pressure head
soil
no.
α
α
θ
θ
θ
θ
n
n
relationship. The RMSE for fitted retention curves
r
s
r
s
−1
−1
was almost identical for both approaches, whilecm
cm
the RMSE for the predicted hydraulic conductiv- 4672
0.000
0.391
0.007
1.15
0.000
0.392
0.007
1.16
ity function was often much lower for the MORE 2242
0.161
0.548
0.024
4.75
0.140
0.551
0.024
4.09
approach than for the RETC approach. Only in 1400
0.000
0.449
0.004
1.09
0.000
0.448
0.003
1.09
one case did the RETC approach provide a lower 4011
0.057
0.415
0.021
1.70
0.052
0.415
0.022
1.66
RMSE than the MORE approach.
1382
0.000
0.344
0.017
1.11
0.000
0.345
0.015
1.12

Case 2: Las Cruces Trench Experiment

4661
4670
1490
1361
1371
1270
2740

0.037
0.000
0.000
0.165
0.208
0.066
0.000

0.408
0.441
0.399
0.433
0.407
0.169
0.710

0.102
0.005
0.008
0.004
0.003
0.002
0.107

1.82
1.41
1.19
1.26
1.38
1.67
1.17

0.036
0.030
0.007
0.192
0.219
0.064
0.000

0.408
0.429
0.400
0.430
0.405
0.169
0.707

0.101
0.004
0.008
0.004
0.002
0.002
0.099

1.81
1.56
1.20
1.32
1.44
1.63
1.17

The RETC and MORE algorithms were used
to estimate the retention curve parameters by fitting the retention curve data for the Las Cruces
Trench experiment. Figure 7 shows the tradeoffs
between different possible parameter sets in theobjective space when fitting the retention curve
and the unsaturated hydraulic conductivity funcRETC. In other words, it may be concluded that the parametertion. The RETC solution and the “best” possible solution fromset estimated using MORE is more representative of the flow
MORE are also compared in the objective space. Application ofprocesses than the one estimated using RETC. One may also
the RETC code to the retention curve data for the undisturbed
observe that the errors in predicted water contents are higher
and disturbed soils (>500 soil samples) resulted in the followfor Day 19 than for Day 35. It is interesting to see an improved
ing retention curve parameter values: θs = 0.321, θr = 0.083,
prediction of the wetting front with the MORE parameters.
α = 0.055 cm−1, and n = 1.51. Note that the estimated values
A simple error analysis of the predicted water contents wasfrom the RETC code are the same as those used by Wierenga
performed for Days 19 and 35 (Table 5). The error statistics
et al. (1991). The MORE algorithm led to the following “best”
indicate a considerable improvement in predicted water conparameter set: θs = 0.323, θr = 0.084, α = 0.038 cm−1, andtents using the parameter set from MORE. The errors in pren = 1.65. The MORE solution is considerably different from
dicted water contents were reduced by 25% using the MORE
the RETC solution in terms of the estimated parameter values.
approach. The cumulative distribution functions of the errors
Of particular interest is the value of n, which has a significant
in predicted water contents for Days 19 and 35 are plotted in
influence on the shape of the unsaturated hydraulic conductivFig. 9. Also, the maximum error in predictions was reduced fority. It should be noted that in Fig. 7, while RMSE(θ) comparesboth days from 0.13 to 0.08 by using the soil hydraulic paramthe observed and predicted water contents, RMSE(k) compares
eters from MORE. The cumulative distribution function for
the apparent observed and predicted relative hydraulic conducMORE is initially steeper than the RETC, indicating prevativities (which were calculated from the optimized parameter
lence of smaller deviations between observed and calculated
set and observed and predicted water contents). Vrugt et al.
water contents for the latter.
(2003) suggested that it is desirable
Table 4. Root mean squared errors between observed water contents or unsaturated hydraulicto have a smooth nondominatedconductivities and corresponding fitted values calculated using the parameter sets estimatedPareto set, which is clearly the caseby the RETC and Multiobjective Retention Curve Estimator (MORE) approaches.
in Fig. 7. The estimated retention
RMSE
curve parameters obtained with the
RETC
MORE
RETC and MORE algorithms were
UNSODA no.
Unsaturated hydraulicUnsaturated hydraulicthen used to predict water content
Water content
Water content
conductivity
conductivity
profiles at two different times for
0.004
3.100
0.004
2.970
which observation data were avail- 4672
2242
0.011
0.536
0.011
0.359
able. The modeling was done using
0.003
0.852
0.003
0.816
HYDRUS-1D (Šimunek et al., 1400
4011
0.010
0.902
0.010
0.832
2005). Figure 8 compares predicted
1382
0.007
1.297
0.008
1.218
water content profiles calculated
4661
0.008
1.661
0.008
1.652
with parameters obtained using
4670
0.019
0.809
0.021
0.764
the RETC and MORE approaches
1490
0.006
0.815
0.006
0.825
against observed water contents for
1361
0.005
0.493
0.005
0.409
Days 19 and 35. The parameter set
1371
0.004
0.298
0.004
0.250
estimated using MORE provides a
1270
0.002
0.717
0.002
0.679
better estimate of water contents
2740
0.031
3.789
0.011
3.590
than the parameter set obtained by
SSSAJ: Volume 72: Number 2 • March–April 2008

477

Cumulative probability

1(a) Day 19
0.8
RETC
MORE
RETC−MO

0.6

RETC

0.4
0.2
0

0

0.02

0.04

0.06
0.08
Absolute error

0.1

0.12

Fig. 7. Evolution of the tradeoff curve in the objective space of theMultiobjective Retention Curve Estimator (MORE) algorithm for
the Las Cruces trench experiment. Gray points represent different
parameter sets considered by the MORE approach. Circles represent the final nondominated Pareto set after 10,000 function
evaluations. The solid square symbol corresponds to the “best”
possible parameter set estimated using the MORE approach; k is
the relative hydraulic conductivity and θ is the volumetric watercontent. The solution set estimated by RETC is less optimal thanthe MORE solution and is beyond the limits of the plot.

Cumulative probability

1(b) Day 35
0.8
0.6

MORE
RETC

0.4
0.2
0

0

0.02

0.04

0.06

0.08

0.1

0.12

Absolute error
Fig. 9. Cumulative distribution functions for absolute errors in predicted water contents for (a) Day 19 and (b) Day 35 for the Las
Cruces trench experiment.

CONCLUSIONS
We have presented a new approach, called MORE, for estimating soil hydraulic parameters when only retention data are
available. This approach provides the optimal parameter set bysimultaneously and independently minimizing two objectivefunctions that represent water content and relative hydraulicconductivity residuals. When observed hydraulic conductivities
are not available, the MORE approach transforms both predicted and observed water contents into the hydraulic conductivity space using Mualem’s pore-size distribution model and also
optimizes the soil hydraulic parameters in this fictitious space.
By doing so, it accounts for the effect the fitted retention function exerts on predicted unsaturated hydraulic conductivities.
We have demonstrated the applicability of the MORE algorithm using two case studies, both of which showed an improve-

ment in descriptions of unsaturated hydraulic conductivities
and soil water flow characterization. The MORE approach is apromising approach for improving the estimation of soil hydraulic parameters, despite some minor drawbacks. For example,
MORE tends to be slower than RETC due to the higher computational demand.
The performance of both approaches was also compared forthe van Genuchten–Burdine model and for conditions of different data availability. It was found that when both retention andhydraulic conductivity data are available, the performance of RETC
and MORE is similar. In scenarios where only retention data are
available (which is the most common case in modeling studies), the
MORE approach estimates the soil hydraulic parameters better than
RETC since, in addition to retention data, it also considers theirtransformation into the hydraulic conductivity space.
During the development of MORE, a number of philosophical questions were encountered that may be of interest for
future research. One of the key issues was the uniqueness ofthe soil hydraulic parameters and its dependence on the typeof modeling problem. In modeling scenarios involving multiple components, such as movement of water, solute, and heat,
it would be interesting to analyze whether multiple objectives
representing different components need to be considered during the parameter estimation. For example, in a coupled water
Table 5. Summary statistics (RMSE and mean absolute error, MAE)of errors in water contents predicted by the RETC and Multiobjective Retention Curve Estimator (MORE) approaches atdifferent times for the Las Cruces trench experiment.
RMSE

Parameter
Fig. 8. Water content measured and predicted using HYDRUS-1D withthe soil hydraulic parameters optimized by the RETC and Multiobjective Retention Curve Estimator (MORE) approaches for (a) Day
19 and (b) Day 35 for the Las Cruces trench experiment.
478

RETC

MORE

MAE
RETC

MORE

Water content, Day 19

0.040

0.033

0.030

0.026

Water content, Day 35

0.039

0.030

0.029

0.022

SSSAJ: Volume 72: Number 2 • March–April 2008

and heat transport modeling (such as Saito et al., 2006), thethird objective representing the thermal conductivity may need
to be included in MORE. The MORE approach has been initially written in the MATLAB environment (The MathWorks,
1999) and will be integrated into the RETC code in the near
future. Interested readers may contact the authors for using the
MORE code.
REFERENCES
Abbaspour, K.C., R. Schulin, and M.Th. van Genuchten. 2001. Estimating
unsaturated soil hydraulic parameters using ant colony optimization.
Adv. Water Resour. 24:827–841.
Brooks, R.H., and A.T. Corey. 1964. Hydraulic properties of porous media.
Hydrol. Pap. 3. Colorado State Univ., Fort Collins.
Gardner, W.R. 1958. Some steady state solutions of unsaturated moisture flow
equations with applications to evaporation from a water table. Soil Sci.
Soc. Am. J. 85:228–232.
Gupta, H.V., L.A. Bastidas, S. Sorooshian, W.J. Shuttleworth, and Z.L. Yang.
1999. Parameter estimation of a land surface scheme using multicriteria
methods. J. Geophys. Res. 104:19491–19503.
Hollenbeck, K.J., J. Šimunek, and M.Th. van Genuchten. 2000. RETMCL:
Incorporating maximum-likelihood estimation principles in the RETC
soil hydraulic parameter estimation code. Comput. Geosci. 26:319–327.
Leij, F.J., W.J. Alves, M.Th. van Genuchten, and J.R. Williams. 1996. The
UNSODA—Unsaturated soil hydraulic database. User’s manual, Version
1.0. Rep. EPA/600/R-96/095. Natl. Risk Manage. Res. Lab., USEPA,
Cincinnati, OH.
Leij, F.J., W.B. Russell, and S.M. Lesch. 1997. Closed-form expressions forwater retention and conductivity data. Ground Water 35:848–858.
Marquardt, D.W. 1963. An algorithm for least-squares estimation of nonlinear
parameters. J. Soc. Ind. Appl. Math. 11:431–441.
Mertens, J., R. Stenger, and G.F. Barkle. 2006. Multiobjective inverse modeling forsoil parameter estimation and model verification. Vadose Zone J. 5:917– 933.
Mualem, Y. 1976. New model for predicting hydraulic conductivity ofunsaturated porous media. Water Resour. Res. 12:513–522.
Nemes, A., M.G. Schaap, F.J. Leij, and J.H.M. Wosten. 2001. Description of
the unsaturated soil hydraulic database UNSODA version 2.0. J. Hydrol.
251:151–162.
Neuman, S.P. 1973. Calibration of distributed parameter groundwater flow
models viewed as a multi-objective decision process under uncertainty.
Water Resour. Res. 9:1006–1021.
Pruess, K. 2004. The TOUGH codes—A family of simulation tools for multiphase

SSSAJ: Volume 72: Number 2 • March–April 2008

flow and transport processes in permeable media. Vadose Zone J. 3:738–746.
Richards, L.A. 1931. Capillary conduction of liquids through porous media.
Physics 1:318–333.
Saito, H., J. Šimunek, and B.P. Mohanty. 2006. Numerical analysis of coupled water,
vapor and heat transport in the vadose zone. Vadose Zone J. 5:784–800.
Schaap, M.G., and M.Th. van Genuchten. 2005. A modified Mualem–van
Genuchten formulation for improved description of the hydraulicconductivity near saturation. Vadose Zone J. 5:27–34.
Schoups, G.H., J.W. Hopmans, C.A. Young, J.A. Vrugt, and W.W. Wallender.
2005. Multi-criteria optimization of a regional spatially-distributed
subsurface water flow model. J. Hydrol. 311:20–48.
Šimunek, J., M. Šejna, and M.Th. van Genuchten. 1999. The HYDRUS-2D
software package for simulating two dimensional movement of water, heat,and multiple solutes in variably saturated media. Version 2.0. IGWMC-TPS53. Int, Ground Water Modeling Ctr., Colorado School of Mines, Golden.
Šimunek, J., M.Th. van Genuchten, and M. Šejna. 2005. The HYDRUS1D software package for simulating the movement of water, heat, andmultiple solutes in variably saturated media. Version 3.0. Dep. ofEnvironmental Sciences, Univ. of California, Riverside.
The MathWorks. 1999. MATLAB, Version 5.3. MathWorks, Natick, MA.
van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.
van Genuchten, M.Th., S.M. Lesch, and S.R. Yates. 1991. The RETC codefor quantifying the hydraulic functions of unsaturated soils. Version 1.0.
U.S. Salinity Lab., Riverside, CA.
Vrugt, J.A., and J.H. Dane. 2005. Inverse modeling of soil hydraulicproperties. p. 1003–1120. In M.G. Anderson et al. (ed.) Encyclopediaof hydrological sciences. John Wiley & Sons, Chichester, UK.
Vrugt, J.A., H.V. Gupta, L.A. Bastidas, W. Bouten, and S. Sorooshian. 2003.
Effective and efficient algorithm for multiobjective optimization of hydrologic
models. Water Resour. Res. 39(8):1214, doi:10.1029/2002WR001642.
Wierenga, P.J., R.G. Hills, and D.B. Hudson. 1991. The Las Cruces trenchsite: Characterization, experimental results, and one-dimensional flow
predictions. Water Resour. Res. 27:2695–2705.
Wierenga, P.J., A.F. Toorman, D.B. Hudson, J. Vinson, M. Nash, and R.G. Hills.
1989. Soil physical properties at the Las Cruces trench site. NUREG/CR5441. Nuclear Regulatory Commission, Washington, DC.
Yapo, P.O., H.V. Gupta, and S. Sorooshian. 1998. Multi-objective global
optimization for hydrologic models. J. Hydrol. 204:83–97.
Yates, S.R., M.Th. van Genuchten, A.W. Warrick, and F.J. Leij. 1992. Analysisof measured, predicted, and estimated hydraulic conductivity using the
RETC computer program. Soil Sci. Soc. Am. J. 56:347–354.

479