Performance Evaluation of Models That
Describe the Soil Water Retention Curve
between Saturation and Oven Dryness
ORIGINAL RESEARCH

Muhammed Khlosi,* Wim M. Cornelis, Ahmed Douaik, Martinus Th. van Genuchten,
and Donald Gabriels
The objective of this work was to evaluate eight closed-form unimodal analytical expressions that describe the soil-water retention curve across the complete range of soil water contents. To meet this objective, the eight models were compared in terms of
their accuracy (RMSE), linearity (coeﬃcient of determination, R2, and adjusted coeﬃcient of determination, R2adj), and prediction potential. The latter was evaluated by correlating the model parameters to basic soil properties. Retention data for 137
undisturbed soils from the Unsaturated Soil Hydraulic Database (UNSODA) were used for the model comparison. The samples
showed considerable diﬀerences in texture, bulk density, and organic matter content. All functions were found to provide relatively realistic ﬁts and anchored the curve at zero soil water content for the coarse-textured soils. The performance criteria weresimilar when averaged across all data sets. The criteria were found to be statistically diﬀerent between the eight models only forthe sandy clay loam soil textural class. An analysis of the individual data sets separately showed that the performance criteria were
statistically diﬀerent between the models for 17 data sets belonging to six diﬀerent textural classes. We found that the Khlosi
model with four parameters was the most consistent among diﬀerent soils. Its prediction potential was also relatively good due
to signiﬁcant correlation between its parameters and basic soil properties.

ABBREVIATIONS: OM, organic matter; SSE, sum of squared errors; SWRC, soil water retention curve; UNSODA, Unsaturated Soil Hydraulic
Database.

The unsaturated soil hydraulic properties are key factors governing the partitioning of rainfall and irrigation into soil water
storage, evapotranspiration, and deep drainage. The hydraulicproperties involve the soil water retention curve (SWRC), which
relates the matric potential (ψ) with the soil water content (θ)
and the hydraulic conductivity function. Discrete (θ, ψ) data sets
can be obtained from either laboratory or ﬁeld measurements, orpredicted from other soil properties using pedotransfer functions
or other approaches. Both methods yield discontinuous sets of
θ–ψ data pairs within the range of matric potentials used for themeasurements. For modeling purposes, a continuous and smooth
representation of the SWRC is preferred, which can be obtainedby ﬁtting a closed-form analytical expression to a discrete data set.
To date, various expressions appear in the literature to represent

the SWRC (e.g., Brooks and Corey, 1964; van Genuchten, 1980;
Kosugi, 1999). Most of the retention models are successful in thewet part of the SWRC. The dry part of the SWRC, however, isequally important in a number of water-related processes aﬀectedby water contents well below the residual value, such as deﬂationof soil particles by wind (Cornelis et al., 2004), microbial activity and N mineralization in soils (De Neve and Hofman, 2002),
CH4 oxidation in soils (De Visscher and Van Cleemput, 2003),
and applications in colloid science (Blunt, 2001).
There is, hence, a pressing need to accurately represent the
SWRC for all matric potentials. In the last few years and decades,
several attempts have been made to represent the complete retention curve (Ross et al., 1991; Campbell and Shiozawa, 1992; Rossi
and Nimmo, 1994; Fayer and Simmons, 1995; Morel-Seytoux
and Nimmo, 1999; Webb, 2000; Groenevelt and Grant, 2004;
Khlosi et al., 2006). All of these models, except the equation of
Groenevelt and Grant (2004), were tested on data reported byCampbell and Shiozawa (1992) and Schoﬁeld (1935), who measured water contents far below −1500 kPa. All models performed
relatively well; however, when testing the Rossi and Nimmo(1994) sum model against data sets in which ψ ranged between−1 and −1500 kPa, Cornelis et al. (2005) found that this model
behaved rather poorly compared with the van Genuchten (1980)
and Kosugi (1999) models, despite its physically realistic shape.
In this study, we therefore compared eight closed-form unimodal
analytical expressions to describe the SWRC across the complete
range of soil water contents. The comparison includes expressions
by Campbell and Shiozawa (1992), Rossi and Nimmo (1994),
Fayer and Simmons (1995), Webb (2000), Groenevelt and Grant

M. Khlosi, W.M. Cornelis, and D. Gabriels, Ghent Univ., Dep. of Soil
Management and Soil Care, Coupure Links 653, B-9000 Ghent,
Belgium; A. Douaik, Research Unit on Environment and Conservationof Natural Resources, National Institute of Agricultural Research, BP
415, Ave. de la Victoire, 10000 Rabat, Morocco; M.Th. van Genuchten,
USDA, U.S. Salinity Lab., 450 W. Big Springs, Riverside, CA 92507.
Received 24 May 2007. *Corresponding author (Muhammed.Khlosi@
UGent.be).
Vadose Zone J. 7:87–96
doi:10.2136/vzj2007.0099© 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.
www.vadosezonejournal.org · Vol. 7, No. 1, February 200887⎡ ln (ψ o )− ln ( ψ ) ⎤⎥θ = θ a ⎢⎢⎥⎢⎣ ln (ψ o )− ln (ψ a )⎥⎦(2004), and Khlosi et al. (2006), which were tested using a limited number of data pairs (e.g., UNSODA), as is most often thecase in practice. Three statistical criteria were considered to deﬁne
the best models: accuracy (RMSE), linearity (R2 and R2adj), andprediction potential (the correlation between model parameters
and basic soil properties).

or simply as⎡
ln ( ψ ) ⎤⎥
θ = θ a ⎢⎢1−⎥⎢⎣ ln (ψ o )⎥⎦Available Soil Water Retention Models
A large number of functions have been proposed over theyears to describe the SWRC across the complete range of soil
water contents. Some of these functions are new, while others areextensions of existing models. The extended functions are mostly
modiﬁcations of the popularly used models by Brooks and Corey(1964), van Genuchten (1980), and Kosugi (1999). These functions, which will be referred to in this study as the BC, VG, andKLN models, respectively, are given by⎪⎧⎪ θ s⎪θ = ⎪⎨⎛ ⎞λ⎪⎪ θ r + (θ s − θ r )⎜⎜ ψ b ⎟⎟⎪⎪⎩⎝⎜ ψ ⎠⎟⎟⎡⎤1−(1 n )1⎢⎥
θ = θ r + (θ s − θ r ) ⎢n⎥⎢⎣ 1 + (α ψ ) ⎥⎦

[1]

[2]

⎡
⎤1 m
⎡
⎤
ln
ψ
(
)
1
⎢
⎥
⎥+ A⎢
θ = θ a ⎢⎢1−4⎥⎥
⎢ 1 + (αψ ) ⎥⎢⎣ ln (ψ o )⎥⎦⎣⎦⎡ ⎛ ψ ⎞⎤⎟⎟ ⎥⎢ ln ⎜⎜⎢ ⎜⎝ ψ ⎠⎟⎟ ⎥
1
md
⎥
[3]
θ = θ r + (θ s − θ r )erfc ⎢⎢2
2 σ ⎥⎥⎢⎢⎥⎣⎢⎦⎥
respectively, where θs is the saturated soil water content; θr is theresidual soil water content; ψb is the air-entry matric potential;ψmd is the matric potential corresponding to the median poreradius; λ, α, and n are curve-ﬁtting parameters related to thepore-size distribution; σ is a dimensionless parameter to characterize the width of the pore-size distribution; and erfc denotes
the complementary error function.
Unfortunately, Eq. [1], [2], and [3] have considerable
difficulty in representing the retention of water as the degreeof saturation approaches zero, often giving an unrealistic
path of the retention curve. To overcome this problem, various improvements have appeared in the literature. A firstattempt to cover the complete retention curve was madeby Ross et al. (1991). They modified Campbell’s (1974)
equation, which is identical to the power function of the
BC model with the residual water content taken as zero, toextend the SWRC to oven dryness. Campbell and Shiozawa(1992) and Schofield (1935) measured water contents of soils
ranging from sand to silty clay at matric potentials far below
ψ = −1500 kPa. Inspection of their data suggests a log-linear
relationship between the matric potential and the water content for matric potentials less than approximately −30 and−1000 kPa for sand and silt loam, respectively (the limiting
values in their data set). Their silty clay soil showed an intermediate value. Based on these observations, Campbell andShiozawa (1992) expressed the (θ, ψ) relationship in the low
potential range as

www.vadosezonejournal.org · Vol. 7, No. 1, February 2008

[5]

where θa is a curve-ﬁtting parameter representing the soil watercontent at ψ = −1 m and ψo is the matric potential at oven dryness. Note that Campbell and Shiozawa (1992) expressed theirmatric potentials in units of meters and consequently θa in theirequations corresponds to the soil water content at −10 kPa or−1 m with a ln|ψa| value equal to zero. The value of ψo depends
on the temperature, pressure, and humidity at which the soil isdried. Assuming a logarithmic behavior in the very dry range of
the SWRC is consistent with the adsorption theory of Bradley(1936), which considers adsorbed molecules to build up in a layered ﬁlm in which the net force of electrical attraction diminisheswith increasing distance from the soil particle (Rossi and Nimmo,
1994). Incorporating Eq. [5] in a VG-type model, Campbell andShiozawa (1992) described the SWRC from saturation to ovendryness asfor ψ ≥ ψ bfor ψ < ψ b

[4]

[6]

where A and m are curve-ﬁtting parameters.
Rossi and Nimmo (1994) created a four-parameter sum
model (RN1 model) and a three-parameter junction model (RN2model) to represent the SWRC across the entire range from saturation to oven dryness. Both models are based on the Campbell(1974) model, with the residual water content taken as zero. Theirfour-parameter sum model (RN1), which consists of two functions joined at one point, was written as⎡⎛ ψ ⎞2 ⎤θ = θ s θ I = θ s ⎢⎢1−β ⎜⎜⎜ ⎟⎟⎟ ⎥⎥ for 0 ≥ ψ ≥ ψ i⎝ ψ b ⎠⎟ ⎥⎢⎣⎦⎡⎛ ⎞λ ⎛ ⎞λ⎛ ψ ⎞⎤ψψθ = θ s θ II = θ s ⎢⎢⎜⎜ b ⎟⎟⎟ − ⎜⎜ b ⎟⎟⎟ + γ ln ⎜⎜ o ⎟⎟⎟⎥⎥⎝⎜ ψ ⎠⎟⎥⎢⎣⎜⎝ ψ ⎠⎟ ⎜⎝ ψ o ⎠⎟⎦for ψ i ≥ ψ ≥ ψ o[7]where ψi is the soil matric potential at the junction point where
the two curves join, and β and γ are shape parameters. The termθI represents the Hutson and Cass (1987) parabolic curve thatjoints the Campbell function (1974) at the junction point ψi. The
Ross et al. (1991) correction is included in the expression for θΙΙ.
Further, using data sets from Schoﬁeld (1935) and Campbell andShiozawa (1992), Rossi and Nimmo (1994) showed that at very
low soil water content, the latter becomes proportional to thelogarithm of the soil matric potential, as can be recognized as wellin θΙΙ. Equation [7] contains seven parameters; however, two ofthem can be determined from conditions that ensure continuityof both Eq. [7] and its ﬁrst derivative with respect to ψi. Here wehave chosen to explicitly determine β and γ as analytical functions of ψb, ψi, ψo, and λ (Cornelis et al., 2005):88⎛ ψ ⎞λ ⎛ ψ ⎞λ⎛ ψ ⎞⎛ ψ ⎞λ1−⎜⎜ b ⎟⎟⎟ + ⎜⎜ b ⎟⎟⎟ + ln ⎜⎜ o ⎟⎟⎟⎜⎜ b ⎟⎟⎟ λ⎟⎜ ψ i ⎠⎟⎝⎜ ψ i ⎠⎟ ⎜⎝ ψ o ⎠⎟⎝⎜ ψ i ⎠⎝β = ψ 2b⎡⎛ ⎞⎤⎢1 + 2 ln ⎜⎜ ψ o ⎟⎟⎥ ψ 2i⎟⎢⎝⎜ ψ i ⎠⎟⎥⎦⎣
⎛ ψ b ⎞⎟λ ⎛ ψ b ⎞⎟λ⎛ ψ b ⎞⎟λ⎜2 − 2 ⎜ ⎟⎟ + 2 ⎜⎜ ⎟⎟ −⎜⎜ ⎟⎟ λ⎜⎝ ψ i ⎠⎟⎜⎝ ψ o ⎠⎟ ⎝⎜ ψ i ⎠⎟γ=⎛ψ ⎞1 + 2 ln ⎜⎜ o ⎟⎟⎟⎜⎝ ψ ⎠⎟
i

Fayer and Simmons (1995) further modiﬁed the Brooks–
Corey and van Genuchten functions by replacing the residual
water content with the adsorption equation of Campbell andShiozawa (1992) to obtain
Brooks–Corey:

[8]

θ=
⎧ θs
for ψ ≥ ψ b⎪⎪⎪⎪⎪⎪⎧⎪ θ s −⎪⎫⎪λ⎪⎪⎛⎜ ψ b ⎞⎟⎨⎪ ⎡⎢ln ( ψ ) ⎤⎥ ⎪⎪ ⎡⎤⎪⎪ θ a 1−⎟for ψ < ψ b
+
⎨
⎬
ln
ψ
⎜
(
)
⎥ ⎪⎜⎝ ψ ⎠⎟⎟⎪⎪ ⎢⎣⎢ ln (ψ o )⎥⎦⎥ ⎪⎪ θ a ⎢⎢1−
⎥ ⎪⎪
⎪
ψ
ln
⎪
(
)
o ⎦⎥ ⎭
⎪
⎪
⎩⎪ ⎣⎢
⎩

[9]

By setting ψo arbitrarily at −106 kPa (Ross et al., 1991; Rossi and
Nimmo, 1994) and using Eq. [8] and [9], the number of model
parameters can be reduced to four.
The three-parameter junction model of Rossi and Nimmo(1994) (RN2 model) consists of three functions, which are continuous at the two points where the functions are joined:⎡⎛ ψ ⎞2 ⎤θ = θ s θ I = θ s ⎢⎢1−β′ ⎜⎜ ⎟⎟⎟ ⎥⎥⎝⎜ ψ b ⎠⎟ ⎥⎢⎣⎦⎛ ψ b ⎞⎟λ
θ = θ s θ II = θ s ⎜⎜ ⎟⎟⎝⎜ ψ ⎠⎟⎛ψ ⎞θ = θ s θ III = θ s γ ′ ln ⎜⎜ o ⎟⎟⎟⎜⎝ ψ ⎠⎟

van Genuchten:θ=⎫⎪⎪⎧ θ s −
⎪⎪ ⎡
⎡
⎤ ⎪⎪
⎤1−(1 n )
ln
ψ
(
)
1
⎪
⎪
⎥ +⎨ ⎡⎢⎥⎤
θ a ⎢⎢1−ln ( ψ ) ⎥ ⎬ ⎢
n⎥
⎥ ⎪ ⎢
⎪
⎣⎢ ln (ψ o )⎦⎥ ⎪⎪ θ a ⎢1− ln (ψ )⎥ ⎪⎪ ⎣⎢ 1 + (αψ ) ⎦⎥o ⎦⎥ ⎭⎪⎩⎪ ⎣⎢for 0 ≥ ψ ≥ ψ ifor ψ i ≥ ψ ≥ ψ j[10]where ψi and ψj are the soil matric potentials at the two junction
points, and β′ and γ′ are shape parameters. To describe the shapeof the SWRC near saturation, Rossi and Nimmo combined theparabolic equation proposed by Hutson and Cass (1987) with
the BC model (as described by θII). The equation for θII is a
power law for ψ smaller than the air-entry value ψb. The simple
power law overestimates the water content at very low matricpotentials. For this reason, a third part, θIII, as proposed by Ross
et al. (1991), was added to obtain a water content of zero at ψo.
In this case, there are six parameters other than θs, as well as four
conditions by imposing continuity of the global function and its
ﬁrst derivative at the two junction points. Four parameters (β′,ψi, ψj, and γ′) can be calculated from analytical functions of theremaining two ﬁtted parameters, ψb and λ:

ψ=
2
⎧
⎪
⎪
M (S e − S em ) [( M + 1) S em − M ](S e − S em ) ⎫
⎪
⎪
⎪1−
⎪
−
2
⎪
⎪
S em
S em (1− S em )
ψm ⎨
⎬
⎪
⎪
⎪
⎪
2
⎪
⎪
a
S
S
S
1
+
−
−
(
)
(
)
⎪
⎪
em
em
⎪ MS e⎪⎩⎭aMS =⎡ M ( M + 1) ( M + 1) S − M ⎤
1
em
⎢
⎥
+
2
⎥
S em (1− S em ) ⎢⎢ 2S em
1
−
S
(
em )
⎣
⎦⎥

−1 λ

[12]

⎛ −1 ⎞
ψ j = ψ o exp ⎜⎜ ⎟⎟⎟
⎝ λ ⎠

[13]

⎛ ψ ⎞λγ ′ = λ exp(1) ⎜⎜⎜ b ⎟⎟⎟⎝ ψ o ⎠⎟

[14]

[18]

In this case, there are certain conditions of continuity and smoothness to be satisﬁed to reattach Eq. [17] with the traditional BC
model for the middle potential range.
Webb (2000) proposed a new approach (W model) to combine the VG model with a dry-region expression. This model does
not necessitate reﬁtting of experimental data and consists of two
regions. Region 1 is an adsorption region described with a linear
function on a semilog plot of log(ψ) vs. soil water content (θ).
Region 2 is a capillary-ﬂow region described with the VG model (orThe three free parameters of the RN2 model are then θs, ψb,and λ. Rossi and Nimmo (1994) obtained better accuracy with
their four-parameter sum model, however, than with their threeparameter junction model.
www.vadosezonejournal.org · Vol. 7, No. 1, February 2008[17]

where Se is eﬀective saturation, Se = (θ − θr)/( θs − θr), Sem is
eﬀective saturation at the matching point, ψm is the corresponding matric potential, M = 1/λ, and aMS is deﬁned as[11]⎛ 2 ⎞⎟
ψ i = ψ b ⎜⎜⎜⎟⎝ 2 + λ ⎠⎟[16]

where θa is a curve-ﬁtting parameter representing the soil watercontent at ψ = −1 kPa, and ψo is the matric potential at ovendryness. Equations [15] and [16] are denoted here as FS1 andFS2, respectively.
Morel-Seytoux and Nimmo (1999) extended the BC modelto oven dryness using the three-parameter junction model (RN2).
They divided the matric potential values into three levels: a lowpotential level (from oven dryness to near ﬁeld capacity), a middle
level (ﬁeld capacity to about air-entry matric potential), and ahigh level (air-entry matric potential to zero suction). For the
low-potential level, a slightly modiﬁed form of the RN2 model(θIII) was selected. For the high-potential level, the following
algebraic relation was adopted:

for ψ j ≥ ψ ≥ ψ o⎛ 2 ⎞⎟( 2+λ ) λβ′ = 0.5λ ⎜⎜⎜⎟⎝ 2 + λ ⎠⎟

[15]

89

any other desired function) in which any previous ﬁtting parameteris retained. The linear relation between θ and ψ on a semilog plot
for the dry region was expressed by Webb (2000) as
log ( ψ )− log (ψ m )
+ θm
μ
where
θ=

expressions. Morel-Seytoux and Nimmo (1999) linked up existing SWRCs with a dry-region expression such that existing ﬁtted
curves can be directly used. To utilize other existing ﬁtted curves,
Khlosi et al. (2006) used the adsorption equation of Campbell
and Shiozawa (1992) to modify the Kosugi (1999) model. The
new expression (KCGS model) hence combines the adsorption
equation of the CS model with the KLN model as[19]

log (ψ o )− log (ψ m )[20]
θm
ψm and θm are, respectively, the matric potential and watercontent at the matching point, and μ is the slope of Eq. [19].
Diﬀerent steps are required for determining the water content atmatching point θm. First, we rewrite Eq. [19] as

θ=

μ =−

ψ = 10μ( θ−θ m ) + log( ψ m )⎡ ⎛ ψ ⎞⎤
⎧θ s −⎫⎢ ln ⎜⎜
⎪⎪
⎡1−
⎤ ⎪⎪
⎟⎟ ⎥
⎢ ⎜⎝ ψ ⎠⎟⎟ ⎥
⎢
⎥ ⎪
⎪
1
md ⎥
⎪
θ a ⎢⎢ ln ( ψ ) ⎥⎥ + ⎪⎨ ⎡⎢ln ( ψ ) ⎤⎥ ⎬ erfc ⎢⎢
⎪
⎪
θ
−
1
σ ⎥⎥
2
2
a
⎢
⎢ ln (ψ o )⎥ ⎪⎢ ln (ψ )⎥ ⎪
⎪
⎢
⎥
⎢
⎥
o
⎣
⎦ ⎪
⎪ ⎣
⎪
⎦⎭
⎩
⎣⎢
⎦⎥

Materials and Methods[21]

Sources of Soil Water Retention Curve Data

Second, the VG model is formulated in terms of θ and its slope
is calculated. Third, the slope of the VG model is combined with
Eq. [21] to give the intercept at ψo as
ψ o = 10log( ψ ) − θ m ⎡⎣ d(log( ψ ) d θ m ⎤⎦Soil water retention curve measurements for a selected set ofundisturbed soils from diﬀerent parts of the world for various soil
types were used in this study. The set consisted of 137 undisturbed
soils selected from the UNSODA database (Nemes et al., 2001).
The 137 soils were selected using the following criteria: (i) SWRC
data were available from at least near saturation to, when possible,
near oven dryness (some of the soils had measurements at matricpotentials far below ψ = −1500 kPa), (ii) nearly all soil texture
classes were represented (Fig. 1), and (iii) their basic soil properties
were known (notably clay, silt, sand, and organic matter content
and bulk density). The organic matter content ranged from 0.7 to
214.0 g kg−1, while bulk densities varied from 0.59 to 1.76 Mg
m−3 (Fig. 2).

[22]

where
d [ log( ψ) ]
=
dθm
⎧⎪
⎪⎫⎪
⎪⎪
⎪⎪
⎪⎪
⎪⎪
⎪⎪
1
⎪⎬
− log ⎨
exp(1)n /( n−1) ⎤
⎡ ⎛
⎪⎪
⎪⎪
⎞
⎥
⎪⎪ ( θ m − θ r )(n −1) ⎢1− ⎜⎜ θ m − θ r ⎟⎟
⎪⎪
⎢
⎥
⎟
⎪⎪
⎪⎪
⎜⎝ θ s − θ r ⎠⎟
⎢
⎥
⎣
⎦
⎩⎪
⎭⎪

[23]

Models for the Data Fit

From available SWRC models, we selected eight models:
those proposed by Campbell and Shiozawa (1992), Rossi and

Finally, the value of θm is determined from Eq. [22] by iteration.
Since the value of ψo in Eq. [22] should equal −106 kPa, one cangive an initial estimate of θm and adjust this estimate until theproper value of ψo is obtained.
Recently, a new three-parameter model for the SWRC wasdeveloped by Groenevelt and Grant (2004) (GG model). Thismodel anchors the curve at zero soil water content using thelogarithmic scale in a model for which the pF, deﬁned as log(−ψ),
with ψ expressed in units of centimeters (Schoﬁeld, 1935), is theindependent variable. They found the equation for the model tobe capable of ﬁtting pF curves with remarkable success across thecomplete range from saturation to oven dryness. The soil watercontent hence would be a function of pF as⎡⎛ −k ⎞⎤⎛ −k ⎞
θ (pF ) = k1 ⎢⎢ exp ⎜⎜ 0η ⎟⎟⎟ − exp ⎜⎜ η0 ⎟⎟⎟⎥⎥ for pF > 0⎜⎝ pF ⎠⎟⎥
⎝ 6.9 ⎠⎢⎣⎦
⎛ −k0 ⎞⎟
θ = k1 exp ⎜⎜for pF = 0⎟⎝ 6.9η ⎠⎟[24]

where k0, k1, and η are the three dimensionless parameters and6.9 is the pF value at oven dryness according to the Schoﬁeld
equation (Schoﬁeld, 1935).
Some of the functions described above utilize existing ﬁtted
curves. The Fayer and Simmons (1995) approach uses existingﬁtted curves to estimate the parameters in their modiﬁed SWRC
www.vadosezonejournal.org · Vol. 7, No. 1, February 2008

[25]

FIG. 1. Textural distribution of 137 soils from the UNSODA database(Nemes et al., 2001) used in model development: sandy (sand andloamy sand), loamy (sandy loam, loam, silt loam, and silt), andclayey (sandy clay loam, silty clay loam, clay loam, sandy clay, silty
clay, and clay) soils.
90Comparison Methods

Nimmo (1994) (RN1 and RN2 models), Fayer and Simmons(1995) (FS1 and FS2), Webb (2000), Groenevelt and Grant(2004), and Khlosi et al. (2006).
The parameters of the chosen models were obtained by nonlinear least-squares ﬁtting of the models to the observed SWRCs.
We used, for this purpose, a quasi-Newton algorithm (Press etal., 1992), which is an iterative method starting with some initial
estimates of the parameters. The approach is based on partitioning of the total sum of squares of the observed values into a partdescribed by the ﬁtted model and a residual part of observed
values around those predicted with the model. The goal of thecurve-ﬁtting process is to ﬁnd an equation that maximizes thesum of squares associated with the model while minimizing theresidual sum of squares or sum of squared errors (SSE). The latter
reﬂects the degree of bias and a contribution of random errors,
conventionally expressed asN⎤− θ fit
SSE = ∑ ⎡⎢ θ obs
j (b) ⎥⎦
⎣ j

2

Goodness-of-Fit Statistics

Various statistical measures can be used to compare the ﬁtting accuracy of the SWRC models. In this study, we used asmeasures the RMSE, the coeﬃcient of determination (R2), andthe adjusted coeﬃcient of determination (R2adj), which were calculated for each soil sample. The mean square error (MSE) and
RMSE were derived from the SSE using[27]

RMSE = MSE[28]

respectively. The RMSE (m3 m−3) is an indication of the overall
error of the evaluated function and should approach zero for thebest model performance. The value of R2 reﬂects the proportion of
the total sum of squares (SST) that is partitioned into the modelsum of squares (SSM) since SST is equal to SSM plus SSE:

[26]

j=1

SSM SST − SSE
SSE
=
= 1−
[29]
SST
SST
SST
An R2 value that approaches unity means that the model explains
most of the variability in the observed data, meaning that SSE isvery small compared with SSM.
An additional measure of ﬁt is the R2adj (Neter et al., 1996),
which is designed to take into account the number of parametersin the model. The R2adj better reﬂects how the degree of correlation between observed and ﬁtted data will change as additionalparameters are added to or deleted from the model. The R2adj is
deﬁned by
R2 =where b is a parameter vector containing the p unknowns thatneed to be estimated; j = 1, 2, ..., N; N is the number of soil-water
retention data for each soil sample (ranging between 5 and 27),
θj is the soil water content corresponding to the jth data pair foreach soil, and “obs” and “ﬁt” denote measured and ﬁtted values,
respectively.
The quasi-Newton routine was performed using the mathematical software program MathCad (PTC Corp., Needham, MA).
Initial estimates for the model parameters in the iterative procedure
were selected using reported literature values for the diﬀerent soils,
if available. When not available, we routinely reran the programwith diﬀerent initial parameter estimates to prevent convergenceto local minima in the objective function. To avoid unrealisticallylarge positive (or even negative) values for θa, θr, α, and ψb in theCS, RN1, FS1, FS2, W, and KCGS models, we constrained theirparameters to α > 0 for the CS model, ψb > ψi for the RN1 model,θr ≥ 0 for the W model, and θa by the range of values found forthe same texture class for the CS, FS2, and KCGS models.

2
Radj
=

( N −1) R 2 − ( p −1)

[30]
N−p
where p is the number of model parameters. The R2adj statistic
can take on any value ≤1, with a value closer to 1 indicating abetter ﬁt.
In general, a model with more parameters can ﬁt the observational data better. Although R2adj is generally one of thebest indicators of the quality of the ﬁt when adding additionalparameters to the SWRC model, overparameterization should be
avoided since it results in a nonidentiﬁable model (i.e., a model
leading to sample conﬁguration probabilities identical to those ofa simpler model with fewer parameters), in large variances of theestimated model parameters for similar soils, or in a high degreeof correlation between the parameters (or low parameter uniqueness) if the number of observations is limited, as is often thecase with laboratory-determined SWRCs (Cornelis et al., 2005).
Moreover, it is advantageous to minimize the number of model
parameters when attempting to predict the SWRC from readily
available data using pedotransfer functions.
To check if the R2, R2adj, and RMSE of the eight models are similar or diﬀerent, a statistical test of signiﬁcance is needed. Regarding
the coeﬃcients of determination or, equivalently, their corresponding
coeﬃcients of correlation, the initial hypothesis is that all eight correlation coeﬃcients are equal (Steel and Torrie, 1980):
H 0 : r1 = r2 = ... = r8

FIG. 2. Variation of bulk density and organic matter content in thedata set.
www.vadosezonejournal.org · Vol. 7, No. 1, February 2008

1
SSE
N

MSE =

91

[31]

with r being either R or Radj, and 1 to 8 referring to the eightmodels. The ﬁrst step is then to transform the correlation coefﬁcients into a new variable z such that:or −1, the higher the prediction potential of the parametersin the model. High correlations between model parameters
and basic soil properties are useful keys in developing reliablepedotransfer functions for the model parameters. The statistical signiﬁcance of the correlation coeﬃcients was tested using
the same procedure as for R2 and R2adj (Eq. [31–34]). Results
are reported only for the mean correlation coeﬃcients for thewhole data sets.

1⎛ 1 + r ⎞⎟z = log e ⎜⎜[32]
⎝ 1− r ⎠⎟⎟2
with loge being the natural logarithm. Since we have the same
number of data for the eight models for a given data set, the variable χ2obs is next calculated asp

χ 2obs = ( N − 3 ) ∑ ( z i − z )2Results and Discussion[33]

Evaluation of the Models

i =1

where z is the arithmetic mean of zi. This variable is subsequently compared with a theoretical value (χ2theor) at the 95%
conﬁdence level having (p − 1) degrees of freedom. Hypothesis[31] is rejected, meaning that the eight correlation coeﬃcientsand thus the coeﬃcients of determination are statistically different, when
χ 2obs ≥ χ 2theor

Parameter values for the diﬀerent models and for diﬀerent
soil textural classes (Table 1) were obtained by curve ﬁtting themodels to the entire data set for each soil textural class. These
values can serve as useful initial estimates when attempting touse one of the evaluated expressions.
Table 2 contains results of the statistical measures computed for each expression to compare their goodness of ﬁt at
three levels: each separate data set (137 values for each statisticalmeasure), each soil textural class (11 values), and all of the data
sets combined (one value). Based on the mean value of the three
performance criteria over all data sets, Table 2 (last column) and
Fig. 3 show that the mean RMSE values varied between 0.0105
and 0.0161. The statistical analysis shows that the eight models
do not diﬀer between each other regarding their RMSEs. The
same table and ﬁgure also show that relatively high values were
found for the mean R2 (>0.97) and mean R2adj (>0.96) for allmodels, thus indicating that all equations can be considered valid.
Indeed, the test of signiﬁcance for both R2 and R2adj further suggests that the eight models do not diﬀer. A more detailed analysis
was done by considering each textural class separately (Table 2).
The test of signiﬁcance shows that the eight models diﬀer in termsof RMSE only for the silty clay loam but not for the remaining
10 textural classes.
Regarding R2 and R2adj, the eight models diﬀer statisticallyfor the class sandy clay loam, while they can be considered togive similar results for the 10 remaining textural classes. Becauseof the diﬀerences in results between all data sets combined andthose from each textural class, we further analyzed our data byconsidering each data set separately and focusing mainly on R2adj.
The results show that for 17 out of the 137 data sets, the eightmodels were statistically diﬀerent, while they gave similar R2adj
values for the remaining 120 data sets classiﬁed as silty clay loam
(one data set), silty clay (two data sets), silt (three data sets), loamy
sand (11 data sets), and clay loam (ﬁve data sets). The 17 data
sets were classiﬁed as sand (three data sets out of 18), sandy loam
(one out of 19), silty loam (two out of 41), loam (ﬁve out of 23),
sandy clay loam (three out of four), and clay (three out of 10).
These results reﬂect the fact that mean values sometimes will notreveal diﬀerences when individual data sets are used. To further
examine the goodness of ﬁt of the various models, we plotted the
distributions of RMSE, R2, and R2adj for all soil samples in Fig. 3.
Each boxplot shows the median (solid line), mean (dotted line),the 25 and 75% percentiles (top and bottom of the box), the 10
and 90% percentiles (whiskers), and the outliers (circles). Noticethat the ﬁtting errors for the FS2, W, and KCGS models were
smaller than those for the other models. This is because of thebetter representation of these models for all soil samples. Results[34]

For RMSE, the test of signiﬁcance is possible for variances. For
this purpose, we ﬁrst square RMSE to calculate MSE, which
represents variances. The hypothesis for homogeneity of variances,also called homoscedasticity (Hartley, 1950), is thenH 0 : σ12 = σ 22 = ... = σ82[35]

where σ2 is the MSEs, with 1 to 8 referring to the eight models.
The Hartley test is now applied by ﬁrst computing the observed
value:
σ 2max[36]
σ 2minin which σ2max and σ2min are, respectively, the maximum andminimum values among the eight MSEs. This observed value iscompared with a theoretical one (Htheor) at the 95% conﬁdence
level with p and (N − 1) degrees of freedom. Hypothesis [35] isrejected, meaning that the eight MSEs and thus the RMSEs are
statistically diﬀerent, when
H obs =

H obs ≥ H theor[37]

To conveniently compare the goodness of ﬁt of the eight models,
RMSE, R2, and R2adj were calculated separately for each modeland for each data set. Next, we calculated the mean values ofthese performance criteria for each of the 11 soil textural classes,
as well as the mean values for the whole data sets for the same
performance criteria.
Comparing Models Parameters with Basic Soil Properties

To provide further insight, the r sp, the Pearson coeﬃ cient of correlation between the model parameters and several
basic soil properties, was computed for each soil sample. Meanvalues were calculated also for each model. Th e soil properties considered here were bulk density, ρb, organic matter
content (OM), and sand, silt, and clay content. The Pearson
coeﬃcient of correlation was used as a measure for the prediction potential of a model in that the closer rsp is to either 1
www.vadosezonejournal.org · Vol. 7, No. 1, February 200892

TABLE 1. Average model parameter values for different soil textural classes.†
Model‡
Data sets, no.
Sample size
CS

Parameter§

RN1

RN2

FS1

FS2

W

GG

KCGS

θa
A
α
m
θs

S
18
9
0.058
0.312
0.253
5.446
0.377

LS
11
10
0.165
0.238
0.528
4.887
0.423

ψb

3.718

1.741

ψi
λ
θs

5.403
1.248
0.378

3.455
1.064
0.426

ψb

3.100

1.672

λ
θa
θs

0.653
0.041
0.373

0.289
0.127
0.419

ψb
λ
θa
θs
α
n
θr
θs
α
n
s
k0
k1
η
θa
θs

3.266
0.945
0.066
0.377
0.227
2.534
0.029
0.379
0.232
2.286
112.482
40.192
0.385
4.202
0.082
0.379

ψmd
σ

10.386
1.086

SL
19
11
0.157
0.247
0.193
10.090
0.414

SiL
41
11
0.199
0.212
0.084
8.265
0.427

Si
3
25
0.216
0.189
0.010
0.964
0.402

L
23
12
0.208
0.268
0.628
17.537
0.501

SCL
4
10
0.268
0.177
0.179
21.072
0.458

SiCL
1
27
0.122
0.327
0.084
11.887
0.447

CL
5
9
0.310
0.274
0.473
9.675
0.639

SiC
2
10
0.400
0.143
0.020
6.813
0.521

C
10
9
0.351
0.170
0.173
6.123
0.548

7.734

11.55227.977

0.433

0.609

0.245

0.886

1.272

0.739

17.6120.468
0.415

16.3180.239
0.425

33.337
2.173
0.403

2.636
10.073
0.498

3.167
8.019
0.454

21.399
1.939
0.447

2.460
1.644
0.639

19.3890.026
0.521

32.616
0.956
0.552

9.004

9.442

15.279

3.636

2.150

12.991

0.537

11.98019.479

0.229
0.118
0.417

0.209
0.160
0.430

0.336
0.040
0.401

0.130
0.290
0.499

0.111
0.283
0.471

0.137
0.189
0.456

0.138
0.123
0.638

0.106
0.470
0.527

0.101
0.304
0.569

2.095
0.620
0.137
0.428
1.564
1.863
0.047
0.430
0.772
1.599
52.626
5.612
9.684
1.872
0.148
0.449

4.439
0.288
0.169
0.421
0.350
1.484
0.041
0.421
0.236
1.375
27.580
9.075
0.834
1.771
0.203
0.429

4.988
0.279
0.145
0.429
0.414
1.372
0.013
0.428
0.183
1.293
26.28512.907
9.269
1.805
0.175
0.447

13.741
0.384
0.078
0.403
0.039
2.835
0.028
0.408
0.041
1.627
52.90926.723
0.434
3.147
0.164
0.406

0.723
0.279
0.300
0.506
11.053
1.394
0.011
0.506
0.977
1.166
17.0465.594
22.130
0.760
0.285
0.542

0.553
0.189
0.213
0.473
1.707
1.170
0.000
0.466
0.420
1.132
15.6047.224
104.901
0.497
0.182
0.550

2.802
0.135
0.347
0.460
1.968
1.406
0.063
0.461
0.170
1.205
17.7928.756
41.471
0.962
0.374
0.463

0.497
0.173
0.376
0.647
2.690
1.815
0.091
0.650
1.386
1.251
14.1974.337
56.121
0.497
0.441
0.648

2.642
0.268
0.328
0.528
0.183
1.208
0.055
0.527
0.098
1.174
11.218
7.326
2.742
0.808
0.430
0.540

0.163
0.137
0.390
0.586
503.170
1.450
0.056
0.564
1.472
1.159
14.5318.201
126.554
0.390
0.350
0.590

7.045
2.427

54.6092.148

108.6262.855

59.1181.291

47.2483.890

54.9086.901

177.609
4.588

2.503
1.432

187.658
3.621

947.4214.973† S, sand; LS, loamy sand; SL, sandy loam; SiL, silt loam; Si, silt; L, loam; SCL, sandy clay loam; SiCL, silty clay loam; CL; clay loam; SiC, silty clay; C, clay.
‡ CS, Campbell and Shiozawa (1992) model; RN1, Rossi and Nimmo sum model (1994); RN2, Rossi and Nimmo junction model (1994); FS1, Fayer and Simmons (1995) modified Brooks-Corey model; Fayer and Simmons (1995) modified van Genuchten model; W, Webb (2000); GG, Groenevelt and Grant (2004); KCGS, Khlosi et al. (2006).
§ Curve-fitting parameters: the units of the parameters are based on soil water content θ expressed in m3 m−3 and ψ in kPa.

TABLE 2. Statistical measures of the models for 11 textural classes† and all data sets.
Model‡
CS

RN1

RN2

FS1

FS2

W

GG

KCGS

Significance¶

SI§
RMSE
R2
R2adj
RMSE
R2
R2adj
RMSE
R2
R2adj
RMSE
R2
R2adj
RMSE
R2
R2adj
RMSE
R2
R2adj
RMSE
R2
R2adj
RMSE
R2
R2adj
RMSE
R2
R2adj

S
0.00920.993
0.989
0.00960.993
0.989
0.01470.987
0.983
0.01080.992
0.987
0.00970.993
0.988
0.01040.992
0.988
0.01460.988
0.984
0.0105
0.992
0.987
No
No
No

LS
0.00950.991
0.986
0.00970.992
0.987
0.01500.986
0.981
0.01080.991
0.986
0.00970.992
0.988
0.01060.991
0.986
0.01490.986
0.982
0.0106
0.992
0.986
No
No
No

SL
0.00920.992
0.987
0.01110.988
0.981
0.01490.983
0.978
0.01270.985
0.976
0.01040.990
0.983
0.01090.989
0.983
0.01430.985
0.980
0.0105
0.990
0.985
No
No
No

SiL
0.00990.990
0.984
0.01390.979
0.968
0.01650.975
0.968
0.01560.976
0.963
0.01260.984
0.975
0.01280.983
0.975
0.01490.981
0.975
0.0113
0.988
0.982
No
No
No

Si
0.01000.989
0.984
0.01410.979
0.968
0.01670.974
0.967
0.01590.975
0.963
0.01280.983
0.975
0.01300.983
0.975
0.01520.980
0.975
0.0115
0.988
0.981
No
No
No

L
0.01000.990
0.985
0.01400.979
0.968
0.01620.975
0.968
0.01540.976
0.964
0.01250.983
0.975
0.01300.982
0.974
0.01460.981
0.976
0.0112
0.988
0.982
No
No
No

SCL
0.01000.989
0.984
0.01400.978
0.966
0.01620.974
0.967
0.01510.976
0.965
0.01230.984
0.976
0.01300.982
0.974
0.01440.981
0.976
0.0111
0.988
0.982
No
Yes
Yes

SiCL
0.00390.999
0.998
0.00430.998
0.998
0.00410.998
0.998
0.00500.998
0.997
0.00390.999
0.998
0.00470.998
0.998
0.00860.993
0.992
0.0047
0.998
0.998
Yes
No
No

CL
0.01000.989
0.983
0.01380.978
0.967
0.01600.975
0.968
0.01490.977
0.966
0.01220.984
0.976
0.01290.983
0.974
0.01440.981
0.976
0.0110
0.988
0.982
No
No
No

SiC
0.0101
0.988
0.982
0.01380.978
0.967
0.01590.975
0.968
0.01480.977
0.966
0.01210.984
0.976
0.01270.983
0.974
0.01420.982
0.976
0.0109
0.988
0.982
No
No
No

C
0.01050.986
0.978
0.01400.975
0.962
0.01600.971
0.963
0.0146
0.978
0.966
0.01200.984
0.976
0.01270.983
0.973
0.01430.980
0.975
0.0109
0.988
0.981
No
No
No

TM
0.01050.985
0.977
0.01400.974
0.961
0.0161
0.971
0.963
0.0146
0.978
0.966
0.01200.984
0.976
0.01280.982
0.973
0.01420.980
0.975
0.0109
0.988
0.981
No
No
No

† S, sand; LS, loamy sand; SL, sandy loam; SiL, silt loam; Si, silt; L, loam; SCL, sandy clay loam; SiCL, silty clay loam; CL; clay loam; SiC, silty clay; C, clay; TM,mean values across the whole data sets.
‡ CS, Campbell and Shiozawa (1992) model; RN1, Rossi and Nimmo sum model (1994); RN2, Rossi and Nimmo junction model (1994); FS1, Fayer and Simmons (1995) modified Brooks-Corey model; Fayer and Simmons (1995) modified van Genuchten model; W, Webb (2000); GG, Groenevelt and Grant (2004); KCGS, Khlosi et al. (2006).
§ SI, statistical indices (mean texture class values).
¶ Test of the equality of the three performance criteria for the eight models: No, no statistical difference between the models; Yes, the performance criteria are statistically different for the eight models.

www.vadosezonejournal.org · Vol. 7, No. 1, February 200893θs. Relatively high values were also obtained for the FS1 andW models. The KCGS model showed high correlations for its
additional parameters as well as for θa and θs. Other models suchas CS, RN1, and FS2 showed good correlation with at least one
additional parameter.
These results indicate that the KCGS model performed best.
This model performed as good as the other models in terms of thegoodness of ﬁt, and showed signiﬁcant correlation between all of
its parameters and the basic soil properties. Compared with the
KCGS model, CS expression exhibited two main disadvantages.
One is that the correlation between model parameters and basic
soil properties was far less than that of the KCGS model forall parameters. A second disadvantage is related to the fact thatthe boxplot of the RMSE, R2, and R2adj values showed manyoutliers for CS compared with KCGS. This suggests that CS isless consistent when applied to diﬀerent soils. The FS2 model,
which showed the same ﬁtting performance as the KCGS andCS models, did not represent the dry range well without having(θ, ψ) data below a matric potential of −100 kPa (Khlosi et al.,2006). Although the W model does not necessitate the reﬁttingof observed data, this model does not require many iterations toﬁnd the matching point.

suggest that the KCGS model is the most suitable for describing
the observed data.
Table 3 lists the Pearson correlation coeﬃcients, rsp, between
all parameters and the basic soil properties. The KCGS and FS2
models produced the highest correlation values for θr (or θa)with ρb and clay content, while most models showed high correlation of θs with ρb and clay content. Notice that θs was highly
correlated to ρb and to a lesser extent to clay content while theopposite was true for θa. The other parameters, which mostly
aﬀect the shape of the SWRC, showed much lower correlations
except for the KCGS, W, FS1, RN2, and GG models. This wasto be expected for the RN2 and GG models, since RN2 does
not contain θr (or θa) while GG does not have θr (or θa) andPerformance of the Models

The performance of the eight closed-form analytical expressions is further demonstrated below by ﬁtting them to soil-water
retention data of a relatively coarse-textured soil (loamy sand,ρb = 1.68 g m−3, OM = 3 g kg−1), a medium-textured soil (siltloam, ρb = 1.51 g m−3, OM = 5 g kg−1), and a ﬁne-textured
soil (clay, ρb = 1.47 g m−3, OM = 6 g kg−1). Figure 4 indicatesthat all models performed very well for the loamy sand in thatthey produced comparatively realistic ﬁts and anchored thecurves at zero soil water content. We can notice that the KCGS
model still showed the best ﬁt for all matric potentials. The
correspondence between observed and ﬁtted SWRCs for thesilt loam exhibited deviations for the RN1, RN2, FS1, W, and
GG models. The RN1 and RN2 models showed a poor matchnear saturation. This is because both models use the Hutson
and Cass (1987) expression at high matric potentials, which isless ﬂexible because of its parabolic shape. The RN1 and RN2models further performed rather poorly in the dry range sincethese functions require many data points to provide a smoothmatch at the junction points. The FS1 model mostly missed
the shape of the data near saturation due to its discontinuous
character, which is an inherent feature of the original BC model.
Figure 4 additionally shows that the GG model did not accurately match several points near saturation, not unlike some ofthe other models. This can be attributed to the discontinuity atψ = 1 cm (pF = 0), which decreases the ﬂexibility of the curvein that region. By contrast, the CS, FS2, W, and KCGS models
showed very good ﬁts to the silt loam data, although CS seems
to have a less realistic shape (linear) in the wet region. Slightly
diﬀerent results were obtained for the clay soil. The RN1 and
RN2 models here again showed poor ﬁts, as explained above.
By keeping the oven-dryness pressure (ψo) as a free parameter,
one can improve their ﬁts but with the possibility of producing
oven-dryness values lower than −106 in ﬁne-textured soils. The
GG model now showed a better match to the data, because it
is more ﬂexible when fewer data are presented near saturation.

FIG. 3. Boxplots of (a) RMSE, (b) R2, and (c) adjusted R2 (R2adj)
values for different soil-water retention curve models: Campbell andShiozawa (CS); Rossi and Nimmo four-parameter sum (RN1); Rossi
and Nimmo three-parameter junction (RN2); Fayer and Simmons
modified Brooks–Corey (FS1); Fayer and Simmons modified van
Genuchten (FD2); Webb (W); Groenevelt and Grant (GG); and Khlosi
et al. (KCGS). The box plots summarize the distribution of RMSE, R2,
and R2adj. The horizontal full line in each box signifies the medianvalue and the mean in a dotted line, whereas the bottom and top of
the box represent the 25th and 75th percentile. The whiskers display
the 10th and 90th percentile, while the points indicate the outliers.
www.vadosezonejournal.org · Vol. 7, No. 1, February 200894

TABLE 3. Pearson correlation coefficients between model parameters† and basic soil properties.

for the soil-water retention curve.
The
expressions were evaluatedModel‡
in terms of their goodness of ﬁ t
CS−0.28*
0.64*−0.37*−0.24*
0.12
−0.12
0.17*
−0.33*
using different statistical indiRN1––−0.79*
0.50*−0.21*
0.33*
−0.04
0.07
−0.26*
0.23*
ces. All eight models deﬁned the
RN2––−0.79*
0.51*
0.27*
0.27*
0.52*
−0.48*––
soil water content vs. soil matricFS1−0.26*
0.50*−0.79*
0.56*−0.28*
0.34*
0.51*−0.42*––
potential relationship below theFS2−0.37*
0.69*−0.78*
0.55*
−0.09
0.30*
0.35*
−0.29*
residual water content. Th e perW
0.05
0.11
−0.79*
0.52*
0.24*
−0.12
0.50*
−0.42*
0.60* −0.52*
formance of the models in terms
GG
–
–
–
–
0.16
−0.15
0.27*
0.43*
0.45* −0.56*
of matching the data varied greatly
depending on the degree of aggreKCGS−0.39*
0.67*−0.73*
0.49*−0.18*
0.39*−0.31*
0.41*––
gationor desegregation of the* Significant correlation at the 0.05 level.
data: individual SWRC data sets,† Curve-fitting parameters.
‡ CS, Campbell and Shiozawa (1992) model; RN1, Rossi and Nimmo sum model (1994); RN2, Rossi and Nimmo juncaverages of individual soil texturaltion model (1994); FS1, Fayer and Simmons (1995) modified Brooks-Corey model; Fayer and Simmons (1995)
classes, or the overall mean. Our
modified van Genuchten model; W, Webb (2000); GG, Groenevelt and Grant (2004); KCGS, Khlosi et al. (2006).
§ a corresponds to α, ψb, ψb, ψb, α, α, k0, and ψb in the CS, RN1, RN2, FS1, FS2, W, GG and KCGS models, respectively.
results show that lumping of the¶ All a values are correlated to sand or silt content, except for the values in italics, which are correlated to clay content.
data without considering their tex# b corresponds to m, λ, λ, λ, n, n, k1, and σ in the CS, RN1, RN2, FS1, FS2, W, GG and KCGS model, respectively.
†† All b values in this column are correlated to sand content, except for the values in italics, which are correlated toturalclass provided similar resultsorganic matter (OM).
for the eight models in terms of‡‡ All b values in this column are correlated to silt content, except for the values in italics, which are correlated to clay
content, and the underlined values, which are correlated to ρb.
their ability to ﬁ t observed data.
§§ c corresponds to ψi, s, and η in the RN1, W and GG model, respectively.
When the textural class was taken
into account, however, the eight
For the same reason, the discontinuous character of the FS1models were found to perform diﬀerently for the sandy clay
model did not seem to be problematic for the clay soil. On theloam class, while being not statistically different for theother hand, the CS model showed a biased curve since θs wasremaining 10 textural classes. An analysis for each data set
not considered here. The FS2, W, and KCGS models showed
separately showed that the eight models behaved diﬀerently
excellent ﬁts also for the clay example.
for 17 of the individual 137 data sets, representing six diﬀerent textural classes. The Khlosi et al. (2006) model with four
Conclusions
parameters was found to be the most consistent for the difUsing 137 soil samples from the UNSODA database, we
ferent soils. Moreover, its prediction potential was relatively
compared eight closed-form unimodal analytical expressions
good because of signiﬁcant correlation between its parameters
θr or θavs. ρbθr or θavs. clayθs or Avs. ρbθs or A vs. a§ vs. sand a§ vs. silt b# vs. sand b# vs. silt, c§§ vs. c§§ vs.
clay
clay
or clay¶or clay¶ or OM†† clay, or ρb‡‡ sand

FIG. 4. Observed and fitted soil-water retention curves (SWRCs) for a loamy sand (LS), a silt loam (SiL), and a clay (C). Models are: Campbell
and Shiozawa (CS); Rossi and Nimmo four-parameter sum (RN1); Rossi and Nimmo three-parameter junction (RN2); Fayer and Simmons modified Brooks–Corey (FS1); Fayer and Simmons modified van Genuchten (FD2); Webb (W); Groenevelt and Grant (GG); Khlosi et al. (KCGS); andvan Genuchten (VG). The curve in the case of the GG model is only defined for potential values lower than −0.1 kPa (or −1 cm).
www.vadosezonejournal.org · Vol. 7, No. 1, February 200895

and basic soil properties. Hence, we recommend this analyticalformula for reliable modeling of the whole range of soil watercontents of unsaturated soils, which can vary substantially inbulk density, soil texture, and organic matter content.

References
Blunt, M.J. 2001. Flow in porous media: Pore-network models and multiphaseﬂow. Curr. Opin. Colloid Interface Sci. 6:197–207.
Bradley, R.S. 1936. Polymolecular adsorbed ﬁlms: I. J. Chem. Soc. 139:1467–
1474.
Brooks, R.H., and A.T. Corey. 1964. Hydraulic properties of porous media. Hydrol. Pap. 3. Civ. Eng. Dep., Colorado State Univ., Fort Collins.
Campbell, G.S. 1974. A simple method for determining unsaturated hydraulicconductivity from moisture retention data. Soil Sci. 117:311–314.
Campbell, G.S., and S. Shiozawa. 1992. Prediction of hydraulic properties ofsoils using particle-size distribution and bulk density data. p. 317–328.
In M.Th. van Genuchten et al. (ed.) Indirect methods for estimating thehydraulic properties of unsaturated soils. Univ. of California, Riverside.
Cornelis, W.M., D. Gabriels, and R. Hartmann. 2004. A conceptual model topredict the deﬂation threshold shear velocity as aﬀected by near-surface
soil water: I. Theory. Soil Sci. Soc. Am. J. 68:1154–1161.
Cornelis, W.M., M. Khlosi, R. Hartmann, M. Van Meirvenne, and B. De Vos.
2005. Comparison of unimodal analytical expressions for the soil-water
retention curve. Soil Sci. Soc. Am. J. 69:1902–1911.
De Neve, S., and G. Hofman. 2002. Quantifying soil water eﬀects on nitrogenmineralization from soil organic matter and from fresh crop residues. Biol.
Fertil. Soils 35:379–386.
De Visscher, A., and O. Van Cleemput. 2003. Simulation model for gas diﬀusion and methane oxidation in landﬁll cover soils. Waste Manage. 23:581–
591.
Fayer, M.J., and C.S. Simmons. 1995. Modiﬁed soil water retention functions
for all matric suctions. Water Resour. Res. 31:1233–1238.
Groenevelt, P.H., and C.D. Grant. 2004. A new model for the soil water retention curve that solves the problem of residual water contents. Eur. J. Soil
Sci. 55:479–485.
Hartley, H. 1950. The maximum F-ratio as a short-cut test for heterogeneity ofvariance. Biometrika 37:308–312.
Hutson, J.L., and A. Cass. 1987. A retentivity function for use in soil water
simulation models. J. Soil Sci. 38:105–113.
Khlosi, M., W.M. Cornelis, D. Gabriels, and G. Sin. 2006. Simple modiﬁcation
to describe the soil water retention curve between saturation and oven dryness. Water Resour. Res. 42:W11501, doi:10.1029/2005WR004699.
Kosugi, K. 1999. General model for unsaturated hydraulic conductivity for soils
with log-normal pore-size distribution. Soil Sci. Soc. Am. J. 63:270–277.
Morel-Seytoux, H.J., and J.R. Nimmo. 1999. Soil water retention and maximum capillary drive from saturation to oven dryness. Water Resour. Res.
35:2031–2041.
Nemes, A., M.G. Schaap, F.J. Leij, and J.H.M. Wösten. 2001. Description of
the unsaturated soil hydraulic database UNSODA version 2.0. J. Hydrol.
251:151–162.
Neter, J., M.H. Kutner, C.J. Nachtsheim, and W. Wasserman. 1996. Appliedlinear statistical methods. 4th ed. Richard D. Irwin, Chicago, IL.
Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery. 1992. Numerical recipes in C. Cambridge Univ. Press, Cambridge, UK.
Ross, P.J., J. Williams, and K.L. Bristow. 1991. Equation for extending waterretention curves to dryness. Soil Sci. Soc. Am. J. 55:923–927.
Rossi, C., and J.R. Nimmo. 1994. Modeling of soil water retention from saturation to oven dryness. Water Resour. Res. 30:701–708.
Schoﬁeld, R.K. 1935. The pF of the water in soil. p. 38–48. In Trans. Int. Congr.
Soil Sci., 3rd, Oxford, UK. July–Aug. 1935. Vol. 2.
Steel, R., and J. Torrie. 1980. Principles and procedures of statistics: A biometrical approach. 2nd ed. McGraw-Hill, New York.
van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.
Webb, S.W. 2000. A simple extension of two-phase characteristic curves to include the dry region. Water Resour. Res. 36:1425–1430.

www.vadosezonejournal.org · Vol. 7, No. 1, February 200896