Hydrol. Earth Syst. Sci., 13, 453–465, 2009www.hydrol-earth-syst-sci.net/13/453/2009/© Author(s) 2009. This work is distributed under
the Creative Commons Attribution 3.0 License.

Hydrology and
Earth System
Sciences

An alternative deterministic method for the spatial interpolation ofwater retention parameters
H. Saito1 , K. Seki2 , and J. Šimůnek3
1 Department of Ecoregion Science, Tokyo University of Agriculture and Technology, Fuchu, Tokyo 183-8509, Japan
2 Faculty of Business Administration, Toyo University, Tokyo 112-8606, Japan
3 Department of Environmental Sciences, University of California, Riverside, CA 92521, USA

Received: 1 August 2008 – Published in Hydrol. Earth Syst. Sci. Discuss.: 3 September 2008
Revised: 31 March 2009 – Accepted: 31 March 2009 – Published: 7 April 2009

Abstract. There are two approaches available for mapping
water retention parameters over the study area using a spatial
interpolation method. (1) Retention models can be first fitted to retention curves available at sampling locations prior
to interpolating model parameters over the study area (the FI
approach). (2) Retention data points can first be interpolated
over the study area before retention model parameters are
fitted (the IF approach). The current study compares the performance of these two approaches in representing the spatial
distribution of water retention curves. Standard geostatistical interpolation methods, i.e., ordinary kriging and indicatorkriging, were used. The data used in this study were obtained
from the Las Cruces trench site database, which contains
water retention data for 448 soil samples. Three standard
water retention models, i.e., Brooks and Corey (BC), vanGenuchten (VG), and Kosugi (KSG), were considered. For
each model, standard validation procedures, i.e., leave-oneout cross-validation and split-sample methods were used toestimate the uncertainty of the parameters at each sampling
location, allowing for the computation of prediction errors(mean absolute error and mean error). The results show that
the IF approach significantly lowered mean absolute errors
for the VG model, while also reducing them moderately for
the KSG and BC models. In addition, the IF approach resulted in less bias than the FI approach, except when the BC
model was used in the split-sample approach. Overall, IF
outperforms FI for all three retention models in describing
the spatial distribution of retention parameters.

Correspondence to: H. Saito(hiros@cc.tuat.ac.jp)

1

Introduction

The predictions of soil moisture distributions in the vadosezone or estimates of contaminant arrival time to groundwater
rely heavily on robust estimates of the spatial distribution ofsoil hydraulic parameters. Spatial interpolation techniques
have been used to estimate the unknown values of these parameters at unsampled locations from available observations.
It has been widely accepted that soil hydraulic parameters
are spatially correlated to a greater or lesser extent. Becauseof this, techniques that take such information into accountmust be used for spatial interpolation. Among many available techniques, including the inverse distance method or the
linear interpolation method, only a least-square interpolation
technique called kriging accounts for spatial correlations between variables. Kriging is now commonly used for mapping
soil physical, chemical, and/or hydraulic properties. Kriging
estimates not only the values of an attribute at unsampled locations, but also their uncertainties in terms of an error variance known as kriging variance, when the underlying geostatistical model is correct (Goovaerts, 1997).
Unlike water flow in saturated systems, predictions ofvariably-saturated water flow in soils depend not only onknowledge of saturated hydraulic conductivities, but also on
knowledge of the water retention and unsaturated hydraulicconductivity characteristics that are usually described by various functional forms. Although experimental data for bothretention and hydraulic conductivity functions are required
for efficient parameterization, in many cases, only retention
data is available. Retention curves are relatively easy to collect when in comparison with unsaturated hydraulic conductivities, which are generally difficult and time-consuming to

Published by Copernicus Publications on behalf of the European Geosciences Union.

1
2
3

Figures454

H. Saito et al.: Spatial interpolation of water retention parametersFit-first Interpolate-later (FI)Interpolate-first Fit-later (IF)Water Retention Curve data

4
5
6

Fit models

Interpolate water contents

Interpolate parameters

Fit models

78Flow and transport simulation
9
10
11

Fig. 1. A difference between the fit-first interpolate-later (FI) andinterpolate-first fit-later (IF) approaches used in this study.
Figure 1:

12measure. Typically, various pore size distribution models(e.g., Mualem, 1976; Burdine, 1953) are then used to predict unsaturated hydraulic conductivity functions from theretention data. In this study, only water retention curve data
will be analyzed, since we have assumed that only this information is available, while information about unsaturated hydraulic conductivities is lacking. The water retention parameters for a particular soil are usually obtained by first measuring a series of pressure head and water content data pairs
from core samples in a laboratory and then by fitting the constructed discrete curve using a simple analytical model, suchas the well-known van Genuchten model (van Genuchten,1980). Water flow in soils is then usually obtained using
a numerical model that simulates
variably-saturated water24
flow, and uses analytical models of soil hydraulic properties. Predictions of water flow are thus greatly improved
when estimated parameters for a given soil hydraulic model
adequately represent the highly nonlinear relationships between water contents, pressure heads and hydraulic conductivities. Although accurate estimation of water retention
curves (WRC) or their model parameters is rarely a goal, it
is one of the most crucial steps in modeling water flow andsolute transport in the vadose zone.
It is well known that soil water retention parameters are
not only soil-type dependent, but also spatial-location dependent. There are two possible approaches to map water retention parameters over the study area using a spatial interpolation method. (1) Retention models are first fitted to experimental retention curves available at sampling locations, andthen the obtained model parameters are interpolated over the
study area. (2) Retention data can first be interpolated overthe study area, and then the retention model is fitted to thisinterpolated data. The former approach is called the FI (f itfirst and interpolate-later) approach, and the latter is referred
to as the IF (interpolate-first and f it-later) approach (Fig. 1).
When non-linear processes are involved in either step, results
differ depending upon the approach chosen (Heuvelink andPebesma, 1999; Leterme et al., 2007).
The FI approach treats the estimated water retention parameters as meaningful physical parameters. However, sincesome water retention parameters are derived in a purely emHydrol. Earth Syst. Sci., 13, 453–465, 2009

pirical way and do not have a physical meaning, the evaluation of these parameters has to be done with great care, or the
results may be misleading and meaningless (Romano, 2004).
However, there are still a number of studies that rely solely
on water retention parameters, because one can dramaticallyreduce the number of parameters that need to be analyzed.
For example, Zhu and Mohanty (2002) averaged the widelyused van Genuchten “parameters” to simulate large-scale infiltration and evaporation processes. In Oliveira et al. (2006),
stochastic fields of “parameters” were generated using conditioning parameters to simulate variably-saturated water flow
in soils. In some studies, the number of parameters is even reduced to one or two by using the scaling method (e.g., Shouse
et al., 1995; Oliveira et al., 2006). The question then arises
whether or not we can consider these derived (i.e., averaged
or randomly generated) parameters to be equivalent to experimentally measured soil physical properties, such as volumetric water contents. In other words, to map water retention
parameters, which approach, FI or IF, should we take?
Although there are a great number of studies that have analyzed and/or modeled the spatial distribution of soil hydraulicproperties, including water retention model parameters (e.g.,
Wendroth et al., 2006; Unlu et al., 1990), we do not know ofany studies that would rigorously investigate this question.
Tomasella et al. (2003) compared two approaches to map retention parameters, in a manner similar to our study. In one
approach, called the point-based approach, water retention
data (the relationship between water contents and water potentials) were predicted using pedotransfer functions (PTF)
prior to fitting retention models, while in the other approach,called the parametric approach, retention model parameters
were directly predicted using PTFs. By using a comprehensive soil water retention database of Brazilian soils, theyfound that the point-based approach outperformed the parametric approach. Sinowski et al. (1997) also compared two
approaches to map water retention parameters using PTFs. In
one approach, soil properties that were used as input parameters of PTFs were interpolated before applying PTFs. In theother approach, water retention parameters obtained throughPTFs were interpolated directly. Their comparison showedthat the former approach outperformed the latter approach.
However, the use of PTFs in predicting retention data introduces additional uncertainty (Schaap et al., 2001) that may
be larger than uncertainties due to the use of fitting and interpolation techniques.
This study therefore aims to quantify how well the estimated soil hydraulic model parameters reproduce the actual
water retention curve at a given location using one of two approaches discussed above (i.e., IF and FI approaches). The
two approaches are evaluated in terms of prediction errors
using standard validation techniques.

www.hydrol-earth-syst-sci.net/13/453/2009/

H. Saito et al.: Spatial interpolation of water retention parameters

455

1
2

0

3
4
5

Table 1. Water retention, Se , models used in this study.

1
2
z (m)

0 cm
0.42
0.4
0.38
0.36
0.34
0.32
0.3
0.28
0.26
0.24
0.22

3

Models4
5

Brooks and Corey (BC)

6

6

0

5

10

x (m)

15

20

25

Expression(
1h≤hb
Se (h) =  h −λh>hbhhb [cm]

b

λ[–]

7
8

Parameters

Fig. 2. A location map of measured saturated water content data
Figure 2:(h=0 cm) at the Las Cruces Trench site. The vertical axis represents
the depth from the surface.

van Genuchten (VG)

Kosugi (KSG)

2
2.1

1
Se (h)= 1+|αh|
n m
]
[

n[–]
α [cm−1 ]m (=1−1/n) [–]no√ hl )Se (h) = 12 erfc ln(h/2σhl [cm]σ [–]Materials and methods
Water retention data and models

The soil hydraulic data set (i.e., water retention curves) usedin this study is available through the Las Cruces Trench Site
database (Wierenga et al., 1989). The database was generated as part of a comprehensive field study conducted insouthern New Mexico, near Las Cruces, for validating andtesting numerical models of water flow and solute transport
in the unsaturated zone (Wierenga et al., 1991). A 24.6-mlong by 6.0-m deep trench wall was excavated, and a total of 450 undisturbed samples, 7.6 cm internal diameter and7.6 cm long, were taken from nine layers, 50 samples weretaken from each layer, for every 0.5 m. Soil samples (undisturbed and disturbed) were then analyzed in the laboratory
for soil properties, such as bulk density, saturated hydraulicconductivity, and soil water retention curve. Soil water retention curves were determined at 448 sampling locations, ui(i=1, 2. . .448) (Fig. 2). At each location, the water contents,
θ ui ; hj , were measured at eleven pressure heads, hj (j=1,
2. . . 11) of 0, −10, −20, −40, −80, −120, −200, −300,
25−1000, −5000, and −15 000 cm
H2 O. While undisturbed
soil cores were used for the wet range (−300 cm to 0 cm),
disturbed soil samples were used with a standard pressure
plate apparatus for the dry range (−15 000 cm to −1000 cm).
More details regarding the experimental procedures can befound in Wierenga et al. (1989). The soil characterizationand infiltration experiments at this site have been analyzed ina number of studies, including recent articles by Rockhold etal. (1996), Oliveira et al. (2006), and Twarakavi et al. (2008).
Figure 2 shows the location map of measured saturated
water contents when the pressure head is equal to 0 cm or h1
at the site. The saturated water content data are, as expected,
spatially heterogeneous with mean, maximum, and minimum
values of 0.322, 0.529, and 0.218, respectively. Although
not shown here, water contents measured for other pressure
heads are also spatially heterogeneous. Figure 3 shows sample histograms of water contents for eleven pressure heads.
The mean of each distribution decreases as the pressure headdecreases, as expected. Not all histograms show highly
skewed or asymmetric distributions, indicating that no datatransformation is required. Figure 4 depicts the experimental
www.hydrol-earth-syst-sci.net/13/453/2009/

water retention curves obtained for seven different depths (z)at x=10.75 m. There are some discontinuities between thewater contents observed for pressure heads of −300 cm, andthose observed for −1000 cm, where water contents increaseas pressure heads decrease. This is likely due to the different
measurement procedures used; undisturbed soil cores wereused for water content measurements from −10 to −300 cm
pressure heads, while disturbed soil cores were used for pressure heads ranging from −1000 to −15 000 cm (Hills et al.,1993). Such inconsistency should have an impact on the results of model fitting, as all data are equally weighted. Investigating the impacts of this inconsistency on final parameter
estimates by FI and IF approaches requires rigorous analysis
of how errors are propagated between different approaches.
Not surprisingly, all seven curves are quite different. Although this variation in the vertical direction is expected to
be much more pronounced than in the horizontal direction,
variations in the horizontal direction also exist and cannot
be ignored, as can been seen from Fig. 2. This confirms theimportance of analyzing the spatial distribution of water retention curves at this site.
Among many available retention curve models, the Brooks
and Corey (BC) model (Brooks and Corey, 1964), the van
Genuchten (VG) model (van Genuchten, 1980), and the Kosugi (KSG) model based on the lognormal pore size distribution (Kosugi, 1996) are the most widely used analytical expressions for representing the dependence of the water content on the capillary pressure head for unimodal pore systems. Detailed discussions of the commonly used soil hydraulic models can be found in Leij et al. (1997). Analytical expressions of the three closed-form models used in thisstudy are listed in Table 1. The parameter Se is an effective
saturation defined as follows:
Se (h) =θ − θrθs − θr(1)

where θ s and θ r are the saturated and residual water contents(m3 m−3 ), respectively. As can be seen from Table 1, in addition to θ s and θ r , all three models require two other shape (orHydrol. Earth Syst. Sci., 13, 453–465, 2009

1

H. Saito et al.: Spatial interpolation of water retention parameters2maximum .53
median .32
minimum .22

.100

.050

.300

.400
0 cm

h=-40cmmaximum .39
median .24
minimum .05

Frequency

9

.100

.100

.200
40 cm

.300

h=-200cm

maximum .32
median .15
minimum .04

.080

.100

h=-5000cm
.200

.200
200 cm16

maximum .17
median .10
minimum .04

Frequency

.100

.080

.200

.300
20 cm

.400

Number of Data 448
mean .17
std. dev. .05
maximum .34
median .17
minimum .05

.120

.080

.000
.000

.300

.100

h=-1000cm

Number of Data 448
mean .14
std. dev. .04.160

maximum .31
median .14
minimum .04

.120

.200
120 cm.300

Number of Data 447
number trimmed 1
mean .12
std. dev. .03
maximum .21
median .12
minimum .05

.080

.040

.050

.100

.150
300 cm

.200

.250

.120

Number of Data 447
number trimmed 1
mean .09
std. dev. .02.080

maximum .15
median .09
minimum .03

.000
.000

.050

.100
.150
1000 cm

.200

.040

.050

.000
.000

17

.200
80 cm

h=-15000cm

Number of Data 447
number trimmed 1
mean .10
std. dev. .02

.150

15

.100

.080

.000
.000

.300

maximum .41
median .27
minimum .07

.040

.040

.000
.000

14

.080

h=-300cm

Frequency

13

.120

.160 h=-120cm

Number of Data 448
mean .19
std. dev. .05

.120

.040

Number of Data 448
mean .27
std. dev. .04

.000
.100

.500

maximum .36
median .19
minimum .05

.000
.000

.400

Number of Data 448
mean .15
std. dev. .04.120Frequency

12

.400

.040

.000
.000

11

.300
10 cm

.120

.050

10

.200

h=-80cm

Number of Data 448
mean .23
std. dev. .05

.200

.160

.040

.000
.100

.500

.150

8

.100
.050

.000
.200

7

maximum .42
median .29
minimum .10.150

Frequency

6h=-20cm

Number of Data 350
number trimmed 98
mean .29
std. dev. .04.200Frequency

5

Frequency.150

4h=-10cm

Number of Data 448
mean .32
std. dev. .03

Frequency

h=0cm

Frequency.200

Frequency

3

Frequency

456

.050

.100
.150
5000 cm

.200

.000
.000

.040

.080
15000 cm

.120

.160

Fig. 3. Sample histograms of water contents for eleven pressure heads (sequentially, h=0, −10, −20, −40, −80, −120, −200, −300, −1000,−5000, and18 −15 000 cm).
1
19
Figure 3:
empirical) parameters, leading to a total of four parameters2
representing
20 the retention curve models.

0.4

Water content [cm3 cm-3]

Application of these water retention functions requires
3 reliable estimates of the parameters for a soil of interest. While
4 pa0.3
θ s can be obtained independently, the other retention
rameters usually have to be estimated indirectly, by fitting
5
the analytical function to the experimental water retention
26
0.2
data using an optimization approach, such as a non-linear
6
least-squares minimization approach (e.g., implemented in
z=0.06m
the RETC code). In most optimization procedures, perforz=0.71m
7
z=2.71m
0.1
mance is improved if the initial estimates are close to thez=3.21mz=3.98m“true” solution. In other words, if initial estimates are
8 too
z=4.89mz=5.79m
far from the “true” solutions, the solution may not converge
0
-1
0
1
2
3
4
9 the
to the global minimum, but to the local minimum during
10
10
10
10
10
10
Suction [cm]optimization process. Seki (2007) recently developed a pro10
gram code that uses a full-automatic procedure to estimate
Fig. 4. Water retention curves obtained for seven different depths at
soil water retention model parameters. The program11autox=10.75 m.
matically selects initial estimates based on observations, sothat the user does not have to make his/her selections.
12 DeFigure 4:
Hydrol. Earth Syst. Sci., 13, 453–465, 2009www.hydrol-earth-syst-sci.net/13/453/2009/

2
3

H. Saito et al.: Spatial interpolation of water retention parameters

.100

.050

.400
qsb

BC model: q_r

8

.200

maximum .14
median .06
minimum .00.150

Frequency

.080

.100
qrb

13

.150

.250

Number of Data 448
mean 15.04
std. dev. 8.59

.200

maximum 64.46
median 14.68
minimum .04

.150
.100

.000
10.0

Frequency

16

30.0
hb

BC model: lambda
.160

15

20.0

40.0

.200

.050

.100
qrv

.150

VG model: alpha

.150

.080

1.00
1.50
lambdab2.00

maximum .16
median .09
minimum .01

.150
.100

.000
.000

.200

.050

.100
qrl

.150

LN model: h.200

Number of Data 448
mean .12
std. dev. 1.34

.160

Number of Data 448
mean 87.57
std. dev. 76.22

maximum 28.16
median .04
minimum .00.120

maximum 816.49
median 64.97
minimum .61

.080

.000

.050

.100
.150
alphav

.200

100.

200.

LN model: sigma

Number of Data 448
mean 1.71
std. dev. .36
maximum 3.60
median 1.64
minimum 1.14

.100

.000
1.00

0.

300.

hl

.120

Number of Data 448
mean 1.58
std. dev. .54
maximum 5.78
median 1.51
minimum .65

.080

.040

.050

.50

Number of Data 448
mean .09
std. dev. .02

.200

.040

VG model: nmaximum 1.77
median .41
minimum .11

.500

.050

.100

.200

.120

.000
.00

.400

LN model: q_r

Number of Data 448
mean .08
std. dev. .03.150

Number of Data 448
mean .46
std. dev. .25

.040

17

.300
qsl

.100

.000
.000

50.0

Frequency

.0

.080

.000
.200

.500

.050

.050

14

.400

maximum .15
median .08
minimum .00

.000
.000

.200

Frequency

Frequency

12

.050

BC model: h

11

.300

.050

.000
.000

maximum .51
median .32
minimum .22.040

VG model: q_r

Number of Data 448
mean .06
std. dev. .03

Number of Data 448
mean .32
std. dev. .03

.120

qsv

.040

10

.080

.000
.200

.500

.120

9

.120

maximum .47
median .32
minimum .22

LN model: q_s.160

Frequency.300

Frequency.160

Number of Data 448
mean .32
std. dev. .03

.040

.000
.200

7

.160

Frequency

maximum .51
median .31
minimum .22.150
Frequency

6

VG model: q_s

Number of Data 448
mean .32
std. dev. .03

Frequency

5

BC model: q_s

Frequency.200

Frequency

4

457

2.00

3.00
nv

4.00

.000
.00

1.00

2.00
sigmal3.00

Fig. 5. Sample
histograms of water retention curve model parameters for BC (left), VG (middle), and KSG (right) models. The saturated
18
and residual water contents are in the top two rows, two shape parameters are in the bottom two rows.
19

Figure 5:

tails of thisapproach will not be discussed in this paper, but
2.2 Geostatistical interpolation theory
20
it should be noted that the full-automatic approach worked
To estimate an unknown value of a given soil property at unvery well in most cases when applied to a number of water
sampled locations, a geostatistical least-square interpolation
retention data from the UNSODA database (Leij et al., 1996).
technique, known as kriging, was used in this study. ConSuch an approach is especially attractive when a great num28
sider the problem of estimating the value of a soil attribute
ber of retention curves need to be fitted simultaneously. Allz (e.g., water content or soil hydraulic parameters) at an unthree models were fitted to 448 retention data points using
sampled location u, where u is a vector of spatial coordinates.
the full-automatic procedure to obtain water retention paramThe available information consists of values of the variable
eters at all sampling locations. Figure 5 shows histograms ofz at N locations ui , i=1, 2. . . N. All univariate kriging estiretention parameters, including θ s and θ r . Similar to the wamatesare variants of the general regression estimate z∗ (u)ter content histograms shown in Fig. 3, θ s and θ r histograms
defined as:show no significant skewness and are more-or-less symmetric, except for θ r of the BC and VG models, where aboutn(u)X
12 and 2.5% of data have a value of 0, respectively. As forz∗ (u) − m (u) =λα (u) [z (uα ) − m (uα )](2)
the shape parameters (e.g., α v and nv ), most are positively
α=1skewed. The parameter hb even shows a bimodal distribuwhere λα (u) is the weight assigned to datum z(uα ) andtion.
m(u) is the trend component of the spatially varying attribute (Goovaerts, 1997). In practice, only those observations closest to u are retained, that is the n(u) data within a
www.hydrol-earth-syst-sci.net/13/453/2009/

Hydrol. Earth Syst. Sci., 13, 453–465, 2009

458

H. Saito et al.: Spatial interpolation of water retention parameters

given neighborhood or window W (u) centered on u, while
the influence of those farther away are discarded (Saito andGoovaerts, 2000). One of the most common kriging estimators is ordinary kriging (or OK), which estimates an unknown
value as a linear combination of neighboring observations:z∗ (u) =

n(u)
X

λOK
α (u) z (uα )

(3)

α=1

In OK, unlike the constant value used in simple kriging,the mean (or trend) at each estimation location (i.e., local
mean, m(u)) is implicitly re-estimated. OK weightsλOK
α aredetermined so as to minimize the error or estimation variance σ 2 (u)= Var{Z ∗ (u) − Z(u)} under the constraint of anunbiased estimate (Goovaerts, 1997). These weights are obtained by solving the system of linear equations known as theordinary kriging system: n(u)P OKλβ (u) γ (uα − uβ ) + µOK (u) = γ (uα − u)







β=1

(4)

λOK
β (u) = 1X
1 n(u)[z (uα ) − z (uα + h)]22N (h) α=1n(u)Xλα (u) γ (uα − u) − µ (u)

(5)

(6)

α=1

When histograms show positively skewed distributions
and/or are not unimodal, OK might not be the best approach to estimate values at unsampled locations (Saito andHydrol. Earth Syst. Sci., 13, 453–465, 2009

λα (u) i (uα ; zk )

(7)

where


where N(h) is the number of data pairs for a given separation vector h. The choice of the model is limited to functions
that ensure a positive definite covariance function matrix ofthe left-hand-side of the kriging system (Eq. 4). Spatial correlations often vary with direction, and such a case requires
one to compute semivariograms in different orientations to
fit anisotropic (direction-dependent) semivariogram models.
In this study, differences in horizontal and vertical variations
observed in retention data were taken into account by considering anisotropic semivariograms. Details of model fitting can be found in Deutsch and Journel (1998), Goovaerts(1997), and Kitanidis (1997). In addition to an estimate of
the unknown z value, OK provides an error variance that iscomputed as
2
σOK
(u) =

n(u)
X
α=1

where µOK (u) is the Lagrange parameter, γ (uα -uβ ) is thesemivariogram between observations at uα and uβ , andγ (uα -u) is the semivariogram between the datum locationuα and the location being estimated, u. The semivariogram
γ (h) models the variability between observations separated
by a vector h. The only information required by the system(Eq. 4) are the semivariogram values, γ , for different separation distances. These values are readily derived from thesemivariogram model fitted to experimental values (i.e., linear model of regionalization):γ̂ (u) =[F (u; zk | (n))]∗ = [Pr {Z (u) ≤ zk | (n)}]∗ =α = 1, ..., n(u)

β=1

n(u)
P

Goovaerts, 2000). In such case, indicator kriging (IK), in
which all data are transformed into either 0 or 1 depending
upon the exceedence of a given threshold (Journel, 1983),
can be used. Indicator kriging allows one to estimate theprobabilities of exceeding a series of thresholds at unsampled locations, leading to the construction of conditional cumulative distribution functions (ccdf) of the variable of interest. Once the ccdf is obtained, an optimal estimate and theuncertainty associated with the estimated value can be quantified by taking the mean (E-type estimate) and variance ofthe ccdf, respectively.
The conditional probabilities at u for a series of thresholds zk discretizing a range of variation of z are estimated asa linear combination of n(u) surrounding indicator transformation of data i(uα ; zk ):i (uα ; zk ) =1 if z (uα ) ≤ zk
0 otherwise

k = 1, ..., K(8)

Kriging weights λα (u; zk ) are then computed from the indicator kriging system, similar to the OK system (Eq. 4),
where all experimental semivariograms γ (u) are replaced byindicator semivariograms γ I (u). A major drawback of IK isthat in order to construct ccdfs at any location u, K indicator
semivariograms, usually 9 to 16 (Deutsch and Journel, 1998),
have to be computed, and K indicator kriging systems haveto be solved. Consequently, IK is much more computationally demanding than OK. It is common to use nine deciles ofthe sample distribution as thresholds. More details on indicator kriging can be found in Goovaerts (1999).
2.3

Prediction performances: FI vs. IF

To compare the FI and IF approaches in terms of reproduction of WRC at i locations ui , Se (h; ui ), two standard validation techniques were used: a leave-one-out cross-validation,in which one observation at a time is temporarily removed
from the dataset and re-estimated from remaining data, anda split-sample method, where the dataset is divided into two
non-overlapping observation and validation sets (Fig. 6). The
number of observation locations in each set were 199 (Np )and 249 (Nv ), respectively. In the split-sample method, theobservation set was used to estimate values at all locations in
the validation set so that estimation errors could be calculated
explicitly. Figure 6 depicts the location map of both observation and validation sets used in the split-sample method.
For complete full validation, different split scenarios had to
be tried. In this study, only the scenario shown in Fig. 6 wasused. This scenario was considered because it is the most
likely and realistic scenario where samples were taken layerby-layer, but not randomly. In addition, the main purpose ofusing the split-sample method was to investigate the impact
www.hydrol-earth-syst-sci.net/13/453/2009/

1

H. Saito et al.: Spatial interpolation of water retention parameters

2

0

459
estimated parameters for the pressure head hj , andθ ui ; hj is the observed water content at the location
ui for the pressure head hj .

Sample set
Validation set

1

4

2
z (m)

3

3
4
5

2.3.2

6

5
6
7

0

5

10

x (m)

15

20

25

1. Unknown water content values corresponding to elevenpressure heads were estimated using OK at 448 locations for the leave-one-out cross-validation and 249 locations for the split-sample method.

Fig. 6. Locations used for sample (closed squares) and validation(open squares) sets in the split-sample
procedure.
Figure 6:8of the sample size. Note that, in this study, the semivariogram models used in kriging were those fitted to semivariograms for all 448 observations. Details about the FI and IF
approaches used in the study are given below.
2.3.1

2. Using the automatic fitting procedure (Seki, 2007), parameters for three WRC models (i.e., BC, VG, andKSG) were obtained at each estimated location. To account for estimation errors when fitting retention models, each h-θ pair was weighted based on the reciprocalof the kriging variance after all kriging variances werestandardized by the total sills of the semivariograms.
More weights are therefore given to those data points
for which the error variances were small.

FI approach

1. Parameters for three water retention functions (BC, VG,
and KSG) were first obtained for all 448 water retention
curves Se (h; ui ) using the automatic fitting procedure(Seki, 2007).

3. For each model, the discrepancy between observed water retention curves and those calculated from parameters obtained in the previous step was obtained by computing MAE and ME, similar to what was done for the
FI approach.

2. For each soil hydraulic parameter, unknown values
were estimated at all 448 locations for the leaveone-out cross-validation, and at 249 locations for the
split-sample method. Following two interpolation approaches were compared:(a)Ordinary kriging (OK) was used for allparameters.

(b)

While OK was used for θ s and θ r ,indicator kriging 29(IK) was used for the
shape parameters.

4. The impact of the number of data pairs used to construct WRC when estimating parameters was investigated by repeating 1 through 3, using only the following six pressure heads: 0, −20, −80, −200, −1000, and−15 000 cm. This approach is referred to as IF6 in theremainder of the paper, while the approach considering
all eleven pressure heads is referred to as IF11. The
number of data pairs used in the IF approach is arbitrary. The optimum number depends on many factors,
including soil type. In this study, instead of finding the
optimum number of data pairs for this particular dataset, we have investigated how prediction errors change
when the number of data pairs is reduced to about halfof the original number.

The former is referred to as the FIOK approach, while
the latter will called the FIIK approach for the remainderof the paper.
3. For each model and approach, a discrepancy between
observed water retention curves Se (h; ui ) and those calculated from estimated parameters Se (h; ui ) was obtained at each location by computing differences in water contents at the eleven corresponding pressure heads.
Mean absolute error (MAE) and mean error (bias ME)
are used in this study in the following way:

MAE =

N
1 X

N i=i

"

11
1 X

11 j =1#θ̂ ui ; hj − θ ui ; hj

(9)

"
#
N
11 


1 X
1 X
θ̂ ui ; hj − θ ui ; hj
ME =
(10)
N i=i 11 j =1
where N is 448 for the leave-one-out cross-validationand 249 for the split-sample method, θ̂ ui ; hj is
the water content at the location ui calculated fromwww.hydrol-earth-syst-sci.net/13/453/2009/

IF approach3

Results and discussions

The results are divided into two sections. The first section
summarizes experimental semivariograms computed fromwater retention data and model parameters. In the second
section, different approaches (IF and FI) are compared interms of describing the spatial distribution of soil hydraulicparameters.
3.1

Semivariograms

Figure 7 shows the experimental semivariograms of WRCmodel parameters with fitted geometric anisotropy models.
Semivariograms of water contents corresponding to eleven
Hydrol. Earth Syst. Sci., 13, 453–465, 2009

460

H. Saito et al.: Spatial interpolation of water retention parameters10.00150.00210.00150.001

0.00140.001

0.0005

5

0
0

θsb
5

Lag [m]

γ

4

γ

3

γ

2

0.00070.0005

θsv
0
0

10

θsl

5

0
0

10

5

Lag [m]10

Lag [m]60.00150.00090.00060.001

0.00060.0004

γ

9

γ

8

γ

7

0.00050.00030.0002

θrb

10

0
0

5

θrv
0
0

10

θrl

5

Lag [m]

0
0

10

5

Lag [m]10

Lag [m]

11
90

4.5

7200

60

3

4800

1.5

2400

14

30

γ

γ

13

γ

12

αv

hb

15

0
0

5

0
0

10

5

10

5

10

0.09

0.18

0.39

0.06

0.12

0.26

γ

Lag [m]γ

Lag [m]γ

Lag [m]

16

hl
0
0

0.03

0.06

0.13

17
18
19
20

λb
0
0

σl

nv

5

10

0
0

Lag [m]5

10

Lag [m]

0
0

5

10

Lag [m]21

Fig. 7. Experimental semivariograms of water retention curve model parameters with fitted geometric anisotropy models for horizontal(circles and the
directions.
Semivariograms were calculated from model parameters at
22 solid line) and vertical (triangles and the dashed line)
Figure7:448 sampling locations. Subscripts b, v, and l correspond to BC (the left column), VG (the middle column), and KSG (the right column)
models, respectively.

pressure heads, h1 − h11 , with fitted geometric anisotropy
models are depicted in Fig. 8. All these semivariograms
were calculated using 448 observations. Except for α v , thereis clear anisotropy in all semivariograms with major spatial continuity observed in the horizontal direction. Whilehorizontal semivariograms, are generally well structured, inmost cases vertical semivariograms fluctuate a lot and are notsmooth. There are not a sufficient number of pairs in the vertical direction to obtain well structured semivariograms. Existence of soil horizons would also not result in spatial correlation of soil variables in the vertical direction at the scale

Hydrol. Earth Syst. Sci., 13, 453–465, 200930 of observations. All experimental horizontal semivariograms

were fitted using either exponential or spherical models, andthey all display a clear nugget effect. Ranges and sill values
vary depending on the variable. For example, as expected,
semivariograms for the saturated (θ s ) and residual (θ r ) volumetric water contents are similar for all three models (Fig. 7,
top two rows). However, other parameters have semivariograms of different shapes, which was expected as the spatial
distribution of these parameters varies. Geometric anisotropy
models, in which the same nugget and sill values are used for
different directions, were fitted for vertical semivariograms.

www.hydrol-earth-syst-sci.net/13/453/2009/

H. Saito et al.: Spatial interpolation of water retention parameters461

10.00150.0021

0.003

0.001

0.0014

0.002

γ

γ

3

γ

2

0.00050.0007

0.001

4

h=0cm
0
0

5

5

h=-10cm
0
0

10

5

10

0.00360.00360.003

70.00240.0024

0.002

γ

8

0.00120.00120.001

6h=-40cm9005h=-80cm0010Lag [m]

10

5

5

10

Lag [m]γ

Lag [m]γ

Lag [m]h=-20cm00h=-120cm0010Lag [m]5

10

Lag [m]0.00240.00180.00120.00160.00120.0008

γ

γ

12

γ

11

0.00080.00060.000413h=-200cm
0
0

14

5

h=-300cm
00100.00060.0006160.00040.0004γ

170.00020.000215h=-5000cm
18

0
0

5

5

10

Lag [m]γ

Lag [m]10h=-1000cm
0
0

5

10

Lag [m]h=-15000cm
00Lag [m]5

10

Lag [m]19

Fig. 8. Experimental semivariograms of water retention data (water contents for a particular pressure head) with fitted geometric anisotropy
models for horizontal (circles and the solid line) and vertical (triangles and the dashed line) directions. Semivariograms were calculated from20
Figure 8:retention data obtained at 448 sampling locations.

Nugget effects may be underestimated using this approach,since the number of data pairs in the vertical direction is not
sufficient. Using a geometric anisotropy model instead of azonal anisotropy model, which can be used to account fordifference in spatial correlation between vertical and horizontal directions, is justified and needed since kriging systems requires the same nugget effect in different directions(Gringarten and Deutsch, 2001). In practice, unless there areconclusive physical explanations for phenomena, it is common to use geometric anisotropy models.

www.hydrol-earth-syst-sci.net/13/453/2009/31

As for the semivariograms of water contents at particular pressure heads (Fig. 8), the ranges of water contents nearsaturation in the horizontal direction (i.e., a major range) are
generally larger than those for drier conditions. This suggests that the spatial continuity of larger pores is more profound than that of smaller pores. When pressure heads are
smaller than −40 cm, the shape of the semivariograms becomes almost identical. From the Laplace equation for capillary rise, the radius of a pore corresponding to a pressure
head of −40 cm can be calculated to be about 0.035 mm (assuming that the contact angle is equal to 0◦ ). Therefore, it can

Hydrol. Earth Syst. Sci., 13, 453–465, 2009

1

462

H. Saito et al.: Spatial interpolation of water retention parameters
1

2

0.4

Water content (cm3 cm-3)

3

5
6
7
8
9

0.3

4

0.2
5

0.1

6

0 -1
10

100

101

11

1

0.05

FI_OK
FI_IK
IF11
IF6

MAE

0.04
0.03
0.02

8
9

0.04
0.03
0.02

0.00
BC

6

0.05

0.01

0.01
0.00

7

8

(b)

(a)

3

5

103

7
104

Fig. 9. Water retention data observed(square), fitted using the KSG
Figure 9:9model (solid black line), and kriged (triangle) at (x, y)=(7.25, 2.16).
In the IF approach, the KSG model was fitted to the constructedwa10
ter retention curve (red). In the FI approach, KSG model parameters
were kriged using ordinary kriging to obtain the FI-KSG curve.
11

2

4

102

Suction [cm]

10

MAE

4

2
Fit
FIOK KSG
IF11 KSG
Observations
Kriged WRC 3

VG

KSG

BC

Model

VG

KSG

Model

Fig. 10. The impact of the choice of the approach on mean absolute errors (MAE) of retention curve
predictions for (a) the leaveFigure 10:
one-out cross-validation and (b) the split-sample procedures. Three
WRC models (BC, VG, and KSG)32were compared.

10

be expected that spatial distributions of pores smaller than
0.035 mm are similar in this soil.
3.2

Prediction performance

Using both the FI and IF approaches, water retention curvemodel parameters can be obtained at any location. Figure 9shows differences between the FI and IF approaches in estimating a water retention curve using the KSG model at a selected location, (x, z)=(7.25, 2.16), where observed retention
data (squares) are available. Due to the exactitude propertyof kriging, data at this location were not included in kriging (i.e., the leave-one-out cross-validation). While the blueline is obtained with parameters predicted by the FIOK approach, the red line is calculated with parameters obtainedby the IF11 approach. In the FIOK approach, model parameters were estimated directly from33 surrounding conditioning
parameters using OK. In the IF11 approach, the water retention curve was first constructed by interpolating water content values corresponding to eleven pressure heads, h1 − h11 .
Constructed water retention data are shown in Fig. 9 using
red symbols. The KSG model (the solid red line in Fig. 9)
was then fitted to the constructed retention data using theHydrol. Earth Syst. Sci., 13, 453–465, 2009

Water Content [m3/m-3]

0.4

Fit (VG model)
Observations

0.3

0.2

0.1

(x, y) = (12.75, 3.21)
0 -3
10

10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

Suction [cm]Fig. 11. Observed water retention
data (square) with the VG model
Figure 11:
fitted at (x, y)=(12.75, 3.21).

aforementioned non-linear least squared method to obtainmodel parameters. While the model obtained using the FIOK
approach could not reproduce the S-shape trend observed inretention data well, the IF11 approach could capture the trendof WRC reasonably well. This confirms that, depending on
the chosen approach, resulting WRC models can be significantly different. In the following section, both approaches
were compared in a more comprehensive manner using theleave-one-out cross-validation and the split-sample methods.
Figure 10 depicts the mean absolute prediction errors(MAE) calculated for all four approaches, i.e., FIOK , FIIK ,IF11, and IF6, using the leave-one-out cross-validation andthe split-sample procedures for the three WRC models. In
general, both IF approaches resulted in smaller MAE thanthe FI approaches, regardless of which WRC model wasused. Decreases in prediction errors were much larger when
the VG model was used, mainly because the prediction performance of the VG model for34 the FI approach was muchworse than that for the other two models. At one location,(x, z)=(12.75, 3.21), the fitted parameter α v (=28.16 cm−1 )was three orders of magnitude greater than in the rest of thedomain, where the mean of fitted α v is 0.12 cm−1 and the
median is 0.04 cm−1 (see the histogram of α v in Fig. 5).
Although the fit was acceptable at this particular location(Fig. 11), the large α v value (an outlier) affected the estimation of α v values at surrounding locations, leading to a
relatively poor performance of the FIOK approach for the VG
model. Indicator kriging was expected to reduce the impactof outliers in the FI approach, as all data are transformed firstto binary data, either 0 or 1, depending upon exceedence ofthe threshold value. However, when E-type estimates werecomputed, the observed extreme value was taken as a possible maximum value. Because of that, the indicator approach
could not improve the prediction. The problem associatedwith the outlier can be avoided if the IF approach is used. For
the other two WRC models, BC and KSG, average prediction errors decreased by about 10–15% when the IF approachwww.hydrol-earth-syst-sci.net/13/453/2009/

H. Saito et al.: Spatial interpolation of water retention parameters

FI_OK
FI_IK
IF11
IF6

0.03
0.02
ME
0.01

5
6
7
8

2

0.01

3

-0.01
BC

VG

KSG

0

2

0.02

0.00

0.00
-0.01

0.03

BC

VG

-0.02

Model

Model

4

4

KSG

5
Fig. 12. The impact of the choice of the approach on meanerrors(ME) of retention curve predictionsfor
Figure12:(a) the leave-one-out cross6 models
validation and (b) the split-sample procedure. Three WRC(BC, VG, and KSG) were compared.

6

(b)

0

5

10

0

5

10

0

5

10

x (m)

15

20

25

15

20

25

15

20

25

0

7
2

was used in comparison with the prediction errors for
the FI8
approach, where there were no obvious outliers in initiallyevaluated parameters.
9
For the FIOK approach, mean absolute errors vary depending on the WRC model used. On the other hand,
10 prediction errors do not depend upon the model used for the IF
approach. This is not a surprising result since three
11 models
were fitted to the same reconstructed WRCs in the IF ap12 for theproach. Therefore, the prediction errors calculated
IF approach are mainly the summary of model fitting performance, but do not directly reflect the accuracy of 13
interpolation as those for the FI approach do.
14 whether
With regard to MAE, there is no distinct trend in
the IF11 approach is better than the IF6 approach is bet15
ter. While IF6 outperforms IF11 when the VG model
wasused for both the leave-one-out cross validation and the split16
sample method, IF11 had slightly smaller mean absoluteer35
rors or almost the same errors when compared to IF6 for the17
BC and KSG models. This indicates that reducing the number of θ-h pairs by about half does not necessarily worsen
18
prediction performances, at least for this particular data set,
as long as the IF approach is adopted. This result is19quite important, since the main reason for preferring the FI approachover the IF approach in most studies is that the numberof20
variables one needs to analyze can be significantly reduced.
The number of parameters required in the WRC models usedin this study is four, which is significantly smaller than theoriginal eleven θ-h pairs used to construct water retention
curves. To use all pairs in the IF approach, one needs to perform kriging eleven times, which requires significantly more
effort than performing it only four times in the FI approach.
However, it was just shown here that when the number ofdata pairs is reduced to six, the IF approach is still better
than the FI approach. With a little extra work, we can improve the description of the spatial distribution of water retention curves using the IF approach. This may not be true
when soils with multimodal pore systems (e.g., aggregatedsoils) are considered. Such soils may require models thatcan distinguish between different pore systems, such as theDurner model (Durner, 1994). To obtain parameters for suchmodels, we need to have more θ-h pairs than required by traditional unimodal models. An appropriate number of pairs
www.hydrol-earth-syst-sci.net/13/453/2009/

z (m)

4

(a)

0.04

ME

3

1

(b)

(a)
0.04

z (m)

2

463

4

6
x (m)

(c)
6

z (m)

1

4

2

0
x (m)

Fig. 13. Locations where IF11 outperformed FI in terms of mean
absolute errors (i.e., smaller MAE) are depicted by open squares,
13: where FI outperformed IF
while closed squares indicateFigure
locationsfor (a) BC, (b) VG, and (c) KSG models.

thus varies depending upon the soil type and can be assessedusing the IF approach. This is another advantage of the IF
approach over the FI approach.
As for the mean error, there is no obvious trend of overor under-estimation, except for the VG model where mean
36
errors are always positive (meaningthat predicted water contents are, on average, greater than observed values). When
comparing different approaches, ME is much smaller (i.e.,
has less bias) for the IF approaches than for the FI approaches, except when the BC model was considered in the
split-sample method. Like that for MAE, this result confirms the superiority of the IF approach over the FI approach.
While for the BC and KSG models, IF11 resulted in less bias
than IF6, for the VG model, IF6 had much smaller biases
than IF11 for both the leave-one-out cross-validation and thesplit-sample method. This indicates again that reducing thenumber of h-θ pairs from eleven to six in the IF approach
does not necessarily worsen the overall performance.
Figure 13 shows spatial distributions of locations where
IF11 outperformed FIOK (open squares), and vice versa
Hydrol. Earth Syst. Sci., 13, 453–465, 2009

464

H. Saito et al.: Spatial interpolation of water retention parameters

Table 2. Percentages (%) of locations where mean absolute errors(MAE) for the IF approach were smaller than those for the FI approach in the leave-one-out cross-validation.
Model
MAE

BC
IF11
71.7

VG
IF6
70.1

IF11
81.0

KSG
IF6
82.1

IF11
62.1

IF6
59.4

(closed squares), in terms of MAE in the leave-one-out crossvalidation. There is no specific trend observed, such as acluster of open or closed squares. This means that the prediction of water retention curves was spatially unbiased. Proportions of open squares for both IF11 and IF6 (spatial distributions are not shown) are summarized in Table 2 for allcases. For all models, the IF approaches outperformed the
FI approaches at more than 60% of the sampling locations.
For the VG model especially, the IF approaches resulted insmaller prediction errors at more than 80% of the sampling
locations. Although the difference between IF11 and IF6 issmall, IF11 was slightly better than IF6 for BC and KSG,
while IF6 was better than IF11 for VG. Overall, these percentages confirm the superiority of the IF approach over the
FI approach in predicting WRCs at any given location.

4

VG model because an extreme value in initially estimatedparameters led to the poor estimate of that parameter at surrounding locations in the FI approach. There is always a riskof including parameter outliers in the estimation process aslong as the FI approach is used. In the IF approach, sucha problem rarely occurs. In addition, the IF approach performed better than the FI approach even when the numberof data pairs used to construct retention curves was reduced
from eleven to six. This makes the IF approach comparable
to the FI approach in terms of workload, because the numberof variables used in the FI approach is four, while that in the
IF approach would be only six. The FI approach has been
preferred mainly because one can reduce the computationaltime and workload. However, as this study shows, with alittle bit of extra work, one can achieve a much better estimation of the spatial distribution of water retention curves byusing the IF approach. In the future, the IF approach should
be applied to estimate unsaturated hydraulic conductivity as
well.
Overall, this study shows the superiority of the IF
approach over the FI approach in evaluating the spatial
distribution of water retention curves, which is the initial but
critical step in modeling variably-saturated water flow andsolute transport in field-scale heterogeneous soils. This study
shows that interpolating retention model parameters, such nand α in the VG model, should be done with extreme care.
In future studies, the actual impact of the FI approach onwater flow and solute transport still needs to be investigated.

Conclusions
Edited by: W. DurnerThis study compared the performance of two interpolation
approaches, a commonly used fit-first interpolate-later (FI)approach and a proposed interpolate-first fit-later (IF) approach, to estimate the spatial distribution of water retention curve parameters. In the former approach, the retention curve model parameters are first evaluated at sampled
locations from observed retention data, and then these parameters are interpolated directly for locations without measurements. In the latter approach, retention data are first interpolated for locations without measurements, and then retention curve model parameters are fitted to these estimatedretention curves. Three retention models, Brooks and Corey,
van Genuchten, and Kosugi models were considered. As for
spatial interpolation, a least-squared interpolation technique,known as kriging, was used. For variables that do not show
asymmetric distributions, ordinary kriging was used, while
for those with asymmetric distributions, indicator kriging
was used as well. Standard validation techniques (the leaveone-out cross-validation and the split-sample methods) wereused to compare the two approaches in terms of reproducing
retention curves.
Both validation results showed that, regardless of the retention model selected, the IF approach had lower prediction errors (mean absolute errors) when compared to the FI
approach. Improvements were especially substantial for theHydrol. Earth Syst. Sci., 13, 453–465, 2009

References
Brooks, R. J. and Corey, A. T.: Hydraulic properties of porous media, Hydrology Paper no. 3, Civil Engineering Department, Colorado State Univ., Fort Collins, Colorado, USA, 1964.
Burdine, N. T.: Relative permeability calculations from pore-size
distribution data. Transactions of the American Institute of Mining, Metall. Petrol. Eng., 198, 71–77, 1957.
Deutsch, C. V. and Journel, A.: GSLIB: Geostatistical Software
Library and User’s Guide, Oxford Univ. Press, New York, USA,
1998.
Durner, W.: Hydraulic conductivity estimation for soils with heterogeneous pore structure, Water Resour. Res., 32(9), 211–223,
1994.
Hills, R. G., Wierenga, P., Luis, S., McLaughlin, D., Rockhold, M.,
Xiang, J., Scanlon, B. and Wittmeyer, G.: INTRAVAL – Phase IImodel testing at the Las Cruces Trench Site, Tech. rep., US Nuclear Regulatory Commission – NUREG/CR-6063, Washington,
DC, USA, 131 pp., 1993.
Goovaerts, P.: Geostatistics for Natural Resources Evaluation. Oxford Univ. Press, New York, USA, 483 pp., 1997.
Goovaerts, P.: Geostatistics in soil science: state-of-the-art and perspectives, Geoderma, 89, 1–45, 1999.
Gringarten, E. and Deutsch, C. V.: Teacher’s Aide, Variogram Interpretation and Modeling, Math. Geol., 33(4): 507–534, 2001.

www.hydrol-earth-syst-sci.net/13/453/2009/

H. Saito et al.: Spatial interpolation of water retention parameters
Heuvelink, G. B. M. and Pebesma, E. J.: Spatial aggregation andsoil process modelling, Geoderma, 89, 47–65, 1999.
Journal, A. G.: Nonparametric estimation of spatial distributions.
Math. Geol. 15(3): 445–468, 1983
Kitanidis, P. K.: Introduction to Geostatistics: Applications in Hydrogeology, Cambridge University Press, New York, USA, 1997.
Kosugi, K.: Lognormal distribution model for unsaturated soil hydraulic properties, Water Resour. Res. 32(9), 2697–2703, 1996.
Leij, F. J., Alves, W. J., van Genuchten, M. T., and Williams, J.
R.: The UNSODA-Unsaturated Soil Hydraulic Database, User’s
manual Version 1.0. Report EPA/600/R-96/095. National Risk
Management Research Laboratory, Office of Research and Development, US Environmental Protection Agency, Cincinnati,
OH, USA, 103 pp., 1996.
Leij, F. J., Russell, W. B., and Lesch, S. M.: Closed-Form Expressions for Water Retention and Conductivity Data, Ground Water,
35(5), 848–858, 1997.
Leterme, B., Vanclooster, M., van der Linden, A. M. A., Tiktak, A.,
and Rounsevell, M. D. A.: The consequences of interpolatingor calculating first on the simulation of pesticide leaching at theregional scale, Geoderma, 137, 414–425, 2007
Mualem, Y.: A new model for predicting the hydraulic conductivityof unsaturated porous media, Water Resour. Res., 12, 513–522,
1976.
Oliveira, L. I., Demond, A. H., Abriola, L. M., and Goovaerts,
P.: Simulation of solute transport in a heterogeneous vadose zone describing the hydraulic properties using a multistep stochastic approach, Water Resour. Res., 42, W05420,
doi:10.1029/2005WR004580, 2006.
Rockhold, M. L., Rossi, R. E., and Hills, R. G.: Application ofsimilar media scaling and conditional simulation for modelingwater flow and tritium transport at the Las Cruces Trench Site,
Water Resour. Res., 32(3), 595–609, 1996.
Romano, N.: Spatial structure of PTF estimates, in: Development ofPedotransfer Functions in Soil Hydrology, edited by: Pachepsky,
Y. A. and Rawls, W. J., 295–319, Elsevier Science, 2004.
Saito, H. and Goovaerts, P.: Geostatistical interpolation of positively skewed and censored data in a dioxin contaminated site,
Environ. Sci. Technol., 34(19), 4228–4235, 2000.
Schaap, M. G., Leij, F. J., and van Genuchten, M. T.: Rosetta: acomputer program for estimating soil hydraulic parameters withhierarchical pedotransfer functions, J. Hydrol., 251, 163–176,
2001.

www.hydrol-earth-syst-sci.net/13/453/2009/

465

Seki, K.: SWRC fit – a nonlinear fitting program with a water retention curve for soils having unimodal and bimodal pore structure,
Hydrol. Earth Syst. Sci. Discuss., 4(1): 407–437, 2007.
Shouse, P. J., Russell, W. B., Burden, D. S., Selim, H. M., Sisson,
J. B., and van Genuchten, M. T.: Spatial variability of soil water
retention functions in a silt loam soil, Soil Sci., 159(1), 1–12,
1995.
Sinowski, W., Scheinost, A. C., and Auerswald, K.: Regionalisationof soil water retention curves in a highly variable soilscape, II.
Comparison of regionalization procedures using a pedotransfer
function, Geoderma 78, 145–159, 1997.
Tomasella, J., Pachepsky, Y., Crestana, S., and Rawls, W. J.: Comparison of two techniques to develop pedotransfer functions forwater retention, Soil Sci. Soc. Am. J., 67, 1085–1092, 2003.
Twarakavi, N. K. C., Saito, H., Simunek, J., and van Genuchten, M.
T.: A New Approach to Estimate Soil Hydraulic Parameters using Only Soil Water Retention Data, Soil Sci. Soc. Am. J., 72(2):471–479, 2008.
Ünlü, K., Nielsen, D. R., Biggar, J. W., and Morkoc. F.: Statisticalparameters characterizing the spatial variability of selected soil
hydraulic properties, Soil Sci. Soc. Am. J., 54, 1537–1547, 1990.
van Genuchten, M. T.: A closed-form equation for predicting thehydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am.
J., 44, 892–898, 1980.
van Genuchten, M. T., Lesch, S. M., and Yates, S. R.: The RETC
code for quantifying the hydraulic functions of unsaturated soils,
Version 1.0, US Salinity Laboratory, USDA, ARS, Riverside,
CA, USA, 85 pp., 1991.
Wendroth, O., Koszinski, S., and Pena-Yewtukhiv, E.: Spatial association between soil hydraulic properties, soil texture and geoelectric resistivity. Vadose Zone Journal, 5, 341–355, 2006.
Wierenga, P. J., Toorman, A. F., Hudson, D. B., Vinson, J., Nash,
M., and Hills, R. G.: Soil physical properties at the Las Cruces
trench site, NUREG/CR-5441, Nucl. Regul. Comm., Washington
DC, USA, 1989.
Wierenga, P. J., Hills, R. G., and Hudson, D. B.: The Las Cruces
trench site: characterization, experimental results, and onedimensional flow predictions, Water Resour. Res., 27, 2695–
2705, 1991.
Zhu, J. and Mohanty, B. P.: Spatial averaging of van Genuchten hydraulic parameters for steady-state flow in heterogeneous soils: anumerical study, Vadose Zone Journal, 1, 261–272, 2002.

Hydrol. Earth Syst. Sci., 13, 453–465, 2009