MulƟmodel SimulaƟon of Water Flowin a Field Soil Using Pedotransfer FuncƟonsO

R

A. K. Guber,* Ya. A. Pachepsky, M. Th. van Genuchten, J. Simunek, D. Jacques,
A. Nemes, T. J. Nicholson, and R. E. Cady

Calibra on of variably saturated flow models with field monitoring data is complicated by the strongly nonlinear dependency of the unsaturated flow parameters on the water content. Combining predic ons using various independent
models, o en called mulƟmodel predicƟon, is becoming a popular modeling technique. The objec ve of this study was
to compare diﬀerent methods of mul model simula on of the field soil water regime using pedotransfer func ons(PTFs). We solved the Richards flow equa on using HYDRUS-1D with parameter sets derived from 19 published PTFs
and compared diﬀerent methods of combining the simula on results from the 19 individual models by (i) using only thebest model, (ii) using equal weights, (iii) regressing measured values to the results of the individual models, (iv) using
singular-value decomposi on (SVD) in the regression, (v) using Bayesian model averaging, and (vi) using weights derived
from Akaike criteria. Data on soil water contents and basic soil proper es at five depths along a 6-m transect in a layered loamy soil were used to calibrate the Richards equa on and to develop the input for the PTFs. The SVD mul model
was the best method, with an accuracy of about 0.01 m3 m−3 at the 35-cm depth and about 0.005 m3 m−3 at greater
depths for 30 d of monitoring and 13 mo of tes ng. This indicates that mul modeling in combina on with monitoringof the soil water regime can be a viable approach to simula ng water flow in the vadose zone.

A: BMA, Bayesian model averaging; PTF, pedotransfer function; SVD, singular-value decomposition.

S

ﬂow simulations have a multitude of applications at various scales, from irrigation scheduling toobtaining global climate change estimates. Modeling of ﬂow invariably saturated soils requires water retention and hydraulicconductivity parameters that are impractical to measure for largescale projects. As opposed to saturated ﬂow, the nonlinearity of
the unsaturated ﬂow constitutive (hydraulic) properties seriously
complicates calibration of the variably saturated ﬂow models
against ﬁeld monitoring data.
Pedotransfer functions are routinely used to relate thehydraulic parameters to readily available data on soil properties

that can be found on soil maps or extracted from soil survey
reports. Since pedotransfer functions are empirical regressiontype relationships, their accuracy outside of their development
region is essentially unknown. A wealth of pedotransfer information has recently accumulated in nearly all parts of the world.
Unfortunately, no good method currently exists to decide whichpedotransfer function model should be used for a speciﬁc site orapplication (Pachepsky and Rawls, 2004).
Climate predictions faced similar uncertainties in model
selection in the 1980s. To deal with these uncertainties, multimodel prediction emerged as a popular technique in climate
prediction (Barnston et al., 2003; Palmer et al., 2000, 2004;
Shukla et al., 2000). The objective of multimodel prediction
is to reduce modeling errors by combining forecasts of various
independent models. Since its introduction by Bates and Granger(1969), multimodel prediction has been subject to much debatethat can be summarized into questions: (i) is a multimodel prediction better than the single best forecast, and (ii) what is thebest approach to weigh predictions obtained with the diﬀerent
models? Regarding the ﬁrst question, several studies (Kharin andZwiers, 2002; Krishnamurti et al., 2000; Palmer et al., 2004)
have reported systematic improvements in the forecasts withmultimodel ensemble prediction, while others have argued thatthe beneﬁt of multimodel prediction in seasonal predictions
is marginal compared with that of using the best single model(Doblas-Reyes et al., 2000; Graham et al., 2000; Peng et al.,2002). Regarding the second question, some studies have demonstrated that proper weight selection leads to relatively better
forecasts compared with simple averaging of the predictions from

A.K. Guber and J. Simunek, Dep. of Environmental Sciences, Univ. of
California, Riverside, CA 92521; A.K. Guber and Y.A. Pachepsky, USDA-ARS,
Environmental Microbial Safety Lab., BARC East, 173 Powder Mill Rd.,
Beltsville, MD 20705; M.Th. van Genuchten, Dep. of Mechanical Engineering, COPPE/LTTC, Federal Univ. of Rio de Janeiro, UFRJ, Rio de Janeiro,
RJ, 21945-970, Brazil; D. Jacques, Belgium Nuclear Research Centre,
SCK-CEN, Mol, Belgium; A. Nemes, USDA-ARS, Crop Systems and Global
Change Lab., Beltsville, MD; and T.J. Nicholson and R.E. Cady, U.S. Nuclear
Regulatory Commission, Oﬃce of Nuclear Regulatory Research, Radioacve Waste, Rockville, MD. Received 17 Aug. 2007. *Corresponding author(Andrey.Guber@ars.usda.gov).
Vadose Zone J. 8:1–10
doi:10.2136/vzj2007.0144© Soil Science Society of America
677 S. Segoe Rd. Madison, WI 53711 USA.
All rights reserved. No part of this periodical may be reproduced or transmi edin any form or by any means, electronic or mechanical, including photocopying,recording, or any informa on storage and retrieval system, without permission
in wri ng from the publisher.

www.vadosezonejournal.org · Vol. 8, No. 1, February 20091

Regression-Based Methods

several models (Krishnamurti et al., 2000; Pavan and DoblasReyes, 2000; Yun et al., 2003), while others have shown thatweighting the forecasts from individual models is not better thansimple averaging (Kharin and Zwiers, 2002; Peng et al., 2002;
Doblas-Reyes et al., 2005). The various studies suggest that theanswer to these two questions may well depend on the numberof data points available for hindcasts (i.e., for evaluation of themodels with data from the past).
The eﬃciency of multimodel predictions was demonstrated
for streamﬂow forecasts by Regonda et al. (2006). Multimodel
prediction methods are now slowly also being used in groundwater modeling. For example, Ye et al. (2004) suggested using theweighted results of several spatial variability models for unsaturated fractured tuﬀ to run ﬂow simulations for situations where
standard information criteria provide an ambiguous ranking of
the models such that it does not justify selecting one of themand discarding all the others. Methods of weighting predictions
obtained from diﬀerent groundwater ﬂow models have been discussed in detail by Poeter and Anderson (2005). Recently, thefeasibility of applying multimodel simulations to water ﬂow inthe vadose zone was demonstrated (Guber et al., 2006).
The objective of this study was to compare diﬀerent methodsof multimodel simulation of the ﬁeld soil water regime using PTFs.
We applied multimodeling using a large database of information on
unsaturated soil water dynamics along a 6-m transect.

If the measured values of the variable to be simulated areY, then the weights wi are determined by treating Eq. [1] as amultiple linear regression equation with S − Y as the dependent
variable and Fi − Fi as the independent variable, and ﬁtting
the linear combination of multimodel predictions ΣiN=1 wi(Fi −Fi ) to the observations Y − Y . Such a method of obtaining
multimodel forecasts has been called superensemble forecasting(Krishnamurti et al., 2000).
The main problem of using regression to search for wi is theexistence of relatively high correlations between Fi values. Having
correlated independent variables in regression, or multicollinearity,
does not preclude using the resultant regression for predictions
within the range of observation (Neter and Wasserman, 1974);
however, it leads to very inaccurate regression coeﬃcients that are
not easily interpreted and a matrix for the system of equations forcalculating weights that is numerically close to singular (Kharin
and Zwiers, 2002). Yun et al. (2003) showed that improvements
in superensemble forecasts can be achieved by applying an SVD
technique to solve the system of equations for the coeﬃcients wi.
Bayesian Model Averaging

Following Raftery et al. (2003), a basic premise is theassumption that for any given forecast, there is a “best” model;
while we do not know what that model is, uncertainty about thebest model can be quantiﬁed using Bayesian model averaging.
Let the bias-corrected value to simulate be Y − Y , and the biascorrected simulations fi = Fi − Fi . Simulation fi is then associatedwith a probability density function (PDF), gi = (Y − Y |fi), whichis interpreted as the PDF of Y − Y conditioned on fi given thatfi is the best simulation. The Baysian model averaging (BMA)
predictive model is

Methods for MulƟmodel SimulaƟons
We consider the case where several PTF models exist to estimate parameters in the constitutive relationships (soil hydraulicproperties) needed for application of the governing Richards
equation to particular variably saturated ﬂow problems. The
multimodel simulation may then mean either (i) combining
the output of several PTFs for the hydraulic parameters into asingle parameter set and then running the ﬂow model with thisparameter set, or (ii) running the ﬂow model with outputs ofindividual PTFs for the hydraulic parameters and then combining the obtained outputs of the ﬂow model. In this study, we
considered multimodel simulations in the sense of the second
option, which is commonly used in meteorological predictions.
Combining simulation results from N models is performed asNp (Y −Y | f 1 , f 2 ,..., f N )= ∑ w i g i (Y −Y | f i )where wi is the posterior probability of the simulation fi being
the best one and is based on the performance of model i in the
training period. The wi are probabilities that add up to one.
Raftery et al. (2003) indicated that it is often reasonable to use
the normal distribution N(fi,σi2) to approximate the conditional
PDF gi = (Y − Y |fi). In that case, the BMA predictive mean is theconditional expectation of Y given the simulations, namelyN

S = Y + ∑ w i (Fi − Fi )

[1]

i =1

N

E ⎡⎣Y −Y | f 1 , f 2 ,..., f N ⎤⎦ = ∑ w i f iwhere S is the multimodel simulation result, Y is the mean of theobserved values during the training period, i is the model (i = 1, 2, …,
N), N is the total number of PTF models, Fi is the simulation result
from model i, Fi is the mean of the simulation results obtained withmodel i during the training period, and wi are the weights on theforecasts of individual models. Equation [1] shows that during thetraining period, the multimodel prediction equation relates deviations from the average of the observed values (S − Y ) with deviations
from the average of the simulated values (F − Fi ).
Existing methods to combine predictions according to Eq.
[1] diﬀer in the way in which the weights wi are obtained. Suchmethods have been reviewed by Clemen (1989), Burnham andAnderson (2002), Armstrong (2001), and Jolliﬀe and Stephenson(2003), among others. The review below brieﬂy outlines the conceptual foundations of some of those methods and includes theapproaches selected for comparison in our study.
www.vadosezonejournal.org · Vol. 8, No. 1, February 2009

[2]

i =1

[3]

i =1

The weights wi and the variances σi2 are found using the maximum likelihood method. The expectation maximization, or EM
algorithm (McLachlan and Krishnan, 1997) can be used for thenormal conditional PDFs (Raftery et al., 2003). The EM algorithm is a method for ﬁnding the maximum likelihood estimatorwhen the problem can be recast in terms of “missing data.” Thisalgorithm is iterative, and alternates between two steps, the E (orexpectation) step, and the M (or maximization) step. Assumingthat the conditional PDF gi = (Y − Y |fi) is computed as
gi =

2

2⎤
⎡
1
⎢ ( f i −Y + Y ) ⎥exp ⎢−⎥⎢⎥
σ i 2π
2σ 2i
⎣⎢
⎦⎥

[4]

the algorithm starts with initial guesses for the weights wi andvariances σi. At the E step, estimates of the missing values zji(m)are computed as

Experimental Data()
∑ g (Y i −Y f ji , σmi −1 )g Y j −Y f ji , σmi −1z (jim ) = N

RMSEs of individual models. Singular-value decomposition was
performed using routines from Press et al. (1992).

The experimental ﬁeld was located at Bekkevoort, Belgium,in a meadow at the bottom of a 4% slope. The soil was classiﬁed
as a Eutric Regosol (FAO, 1975) or a Udiﬂuvent (Soil Survey
Staﬀ, 1999). The top 1 m included three soil horizons: an Ap
horizon between 0 and 25 cm, a C1 horizon between 25 and 55
cm, and a C2 horizon between 55 and 100 cm. A trench, 1.2 m
deep and 8 m long, was excavated at the ﬁeld site. Soil texture
was measured with the pipette method following pretreatment
with sodium hexametaphosphate. Textural classes were loam at
the 15-, 35-, and 55-cm sampling depths and silty loam at 75and 95-cm depths (Table 1). The grass cover was removed fromthe experimental area. A plastic sheet covered the side of thetrench along which instrumentation was installed. Volumetricwater contents were measured with time domain reﬂectometry(TDR). Sixty TDR probes (two rods 25 cm long, 0.5-cm rod
diameter, 2.5-cm rod spacing) were installed along the trench at
12 locations spaced 50 cm apart laterally, at ﬁve depths of 15, 35,
55, 75, and 95 cm. The TDR measurements were performed witha Tektronix (Beaverton, OR) 1502B cable tester. The automatedsystem of Heimovaara and Bouten (1990) was used to control,
retrieve, store, and analyze measurements of the travel time of theelectromagnetic waves along the TDR rods. The time diﬀerence
between two measurements for the same probe was 2 h. Rainfall
was measured and recorded continuously near the trench acrossa catch area of 200 cm2 (Fig. 1).
The trench was ﬁlled after all devices were installed. A thinlayer of gravel (1–2 cm) was evenly distributed over the study area
to (i) decrease the erosive eﬀect of rain impact on the bare soil
surface, (ii) minimize evaporation from the soil surface, and (iii)
decrease the growth of weeds on the experimental plot. Weeds
were regularly removed from the site during the summer. Field
measurements started on 11 Mar. 1998 (Day 0). We used 250 dof data. Additional details about the site and the experiment are
given by Jacques (2000) and Jacques et al. (2001).

[5]

i =1

where j is the number of the observation within the training
period (j = 1, 2, …, K), i is the number of the model, and m is thenumber of the iteration. The M step then consists of estimating
weights and variances asw i(m ) =1 K (m )∑ z jiK j =1σ(im ) =∑ j =1 z (jim ) (Y j −Y − f ji )K
∑ j =1 z (jim )

[6]
2

K

[7]

In general, Bayesian model averaging is designed and has beenused to generate probabilistic forecasts. For that reason, thisapproach appears to have wider applicability in probabilistic forecasting than deterministic regression-based averaging(Rajagopalan et al., 2002).
Using Informa on Theory

Following Burnham and Anderson (2002), Poeter andAnderson (2005) suggested using the Akaike information criterion to derive the weights for individual models aswi =exp (− 0.5Δ i )[8]∑ l =1 exp (− 0.5Δl )Nwhere
Δ i = AICc i − AICc min[9]

AICci is the value of the Akaike criterion for the ith model, andAICcmin is the minimum AICc value of all models in the set. The
AICc value for a model is given by

Pedotransfer Func ons and Flow Simula ons

The literature was searched for pedotransfer functions toestimate the soil water retention and hydraulic conductivity properties from soils data available at the site. We used only PTFs
that had been developed from relatively large (>200 samples)
databases. To estimate soil water retention, we selected 19 PTFs
developed in diﬀerent regions. Equations for the PTFs were takenfrom the appendices of Guber et al. (2006) and Pachepsky etal. (2006). Five of the PTFs (those by Campbell and Shiozawa[1992], Mayr and Jarvis [1999], Rawls and Brakensiek [1985],

2k (k + 1)
[10]
n − k −1
where n is the number of observations, σ2 is the residual variance estimated as the sum of squared residuals divided by n, andk is the number of estimated parameters. This methodology is
eﬃcient if models require estimating their parameters before theprediction averaging.

( )AICc = n log σ 2 + 2k +Mul model Simula ons

The following multimodel prediction methods were compared in our study: (i) using only
the best model, (ii) assigning equal weights to allmodels, (iii) using the superensemble, (iv) using
the superensemble with SVD to ﬁnd weights, (v)using Bayesian model averaging, and (vi) using
information theory. Initial estimates of the BMA
weights were taken proportional to the inverse
www.vadosezonejournal.org · Vol. 8, No. 1, February 2009T

1. Average values of soil proper es at the monitoring depths.
Depth

cm
15
35
55
75
95

No. ofsamples
7
8
5
3
4

3

Soil textural frac on>50 μm 50–20 μm 20–10 μm 10–2 μm <2 μm——————————% (w/w) ——————————
58.6
19.3
6.4
4.5
11.1
56.7
18.9
7.8
3.2
13.3
57.3
17.6
6.6
3.7
14.8
49.6
21.2
7.9
4.4
17.4
43.8
30.3
7.4
4.5
14.0

Bulk
density

Organic C

g cm−3
1.42
1.54
1.53
1.53
1.53

% (w/w)
2.2
0.8
0.4
0.3
0.6

Saxton et al. [1986], and Williams et al. [1992]) estimated the parameters of the equation of Brooks andCorey (1964):⎧⎪ ⎛ h ⎞λθ − θ r ⎪⎪⎪ ⎜⎜⎜ b ⎟⎟⎟ , h > hb=⎨ ⎝ h ⎠φ− θ r ⎪⎪h ≤ hb⎪⎪⎩1,[11]

while four other equations (Varallyay et al., 1982;
Vereecken et al., 1989; Wösten et al., 1999) estimated the parameters of the van Genuchten (1980)
equationθ − θr1=n ⎤mθs − θr ⎡⎢⎣1 + (α h ) ⎥⎦F . 1. Observed (symbols) and simulated (using the calibrated model, lines)
daily average water contents. Observed values are averages across the transect
for each depth and each day.
T

2. Enumera on of pedotransfer func ons (PTFs) inves gated in this work.
No.

1
2
3
4
5
6
7
8
9
10
11
12
13
15
16
17
18
18
19

PTF
Rawls et al. (1982)
Saxton et al. (1986)
Williams et al. (1992)
Campbell and Shiozawa (1992)
Rawls et al. (1983)
Verekeen et al. (1989)
Gupta and Larson (1979)
Baumer (1992)
Varallyay et al. (1982)
Canarache (1993)
Peterson et al. (1968)
Rajkai and Várallyay (1992)
Tomasella and Hodne (1998)
Rawls and Brakensiek (1985)
Hall et al. (1977)
Williams et al. (1992)
Wösten et al. (1999)
Mayr and Jarvis (1999)
Rose a (Schaap et al., 2001)

Clay

Silt

+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

+

+
+

PTF predictors
Bulk
Organic
Sanddensity ma er
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

+

+
+
+
+
+
+

+
+
+
+
+
+

Itop†

+

+
+
+
+
+

+
+

† Variable equal to zero for subsoil and one for topsoil.

www.vadosezonejournal.org · Vol. 8, No. 1, February 2009

4

+

[12]

where θ is the volumetric water content, h is the pressure
head (taken here positive for unsaturated conditions),
φ is the porosity, θr is the residual water content, hb
is bubbling pressure, λ is the pore size distribution
index, θs is the saturated water content, and α, m, andn are empirical shape-deﬁning parameters. Nine pedotransfer functions (Baumer, 1992; Canarache, 1993;
Gupta and Larson, 1979; Hall et al., 1977; Petersen
et al., 1968; Rajkai and Várallyay, 1992; Rawls et al.,
1982, 1983; Tomasella and Hodnett, 1998) estimatedwater contents at several ﬁxed pressure heads. The van
Genuchten parameters θs, θr, α, and n in these cases
were evaluated by ﬁtting Eq. [12] to the water retention
points obtained from those PTFs. The value of parameter m was calculated as m = 1 − 1/n. The residual watercontent θr was set to 0.001 m3 m−3 and the saturated
water content θs to the porosity φ for PTFs that evaluated values of the water content at only two pressure
heads (15,000 and 330 cm). Two PTFs (Vereecken etal., 1989; Varallyay et al., 1982) used van Genuchten
parameters assuming m = 1. To use these PTFs for the
ﬂow simulations, water contents were calculated at
the pressure heads used by those researchers to derive
the PTFs, after which we estimated θs, θr, α, and n byﬁtting Eq. [12], with m = 1 − 1/n, to the calculated
water retention data. We also used the Rosetta software(Schaap, 2004) to generate van Genuchten parametersfrom texture and bulk density. Pedotransfer functions
were enumerated as shown in Table 2. A FORTRAN
code to estimate water retention with the pedotransfer
functions of this study, except Rosetta, is available on
request from the corresponding author.
The saturated hydraulic conductivity, Ksat, has also
been estimated from a PTF. The literature contains only
three Ksat PTFs that were developed or tested usinglarge databases (Rawls et al., 1998; Wösten et al., 1999;
Schaap et al., 2001). We focused on the multimodelingwith diﬀerent water retention PTFs, and used only the
PTF developed by Rawls et al. (1998).
The HYDRUS-1D software (Simunek et al., 1998)
was used to run the simulations. This software gives
options to run simulations either with the Brooks–
Corey water retention Eq. [11] or the van Genuchten

T
3. Calibrated van Genuchten (Eq. [12]) and van Genuchten–Mualem (Eq. [14]) parameters for the Richards flow equa on, and associated RMSEs for diﬀerent soil horizons.
Depth
cm
0–25
25–45
45–65
65–85
85–105θrθs————m3 m−3 ————
0.005 ± 0.220
0.363 ± 0.001
0.086 ± 1.209
0.362 ± 0.008
0.154 ± 1.150
0.367 ± 0.004
0.108 ± 1.443
0.388 ± 0.008
0.010 ± 0.453
0.386 ± 0.004

αcm−1
0.00110 ± 0.0004
0.00184 ± 0.0040
0.00123 ± 0.0044
0.00133 ± 0.0042
0.00155 ± 0.0016

[13]

m⎤
⎡
K (h ) = K sat S el ⎢1− 1− S e1/m ⎥
[14]
⎢⎣
⎥⎦
for the van Genuchten–Mualem model (van Genuchten, 1980),
where

Se =

)

θ − θrθs − θr1.80 ± 0.19
1.27 ± 1.14
1.53 ± 1.35
1.49 ± 1.44
1.20 ± 0.20

l

RMSE

0.01 ± 0.83
0.50 ± 36.4
5.01 ± 54.4
0.5 ± 24.9
18.4 ± 31.5

m3 m−30.00880.00550.00380.00530.0043Observed and simulated soil water contents are shown in
Fig. 1. The topsoil exhibited more weather-related variations in
the water content than the subsoil. Still, the relatively dry periods
between Days 55 and 76 and between Days 140 and 166 are
reﬂected in lower soil water contents at depths down to 75 cm.
Calibrations with the Richards ﬂow model for the ﬁve-layer
proﬁle resulted in relatively high accuracy (Fig. 1). The RMSE
values decreased with depth (Table 3); calibrated parameters alongwith their standard errors and linear tolerance intervals are shownin Table 3. The parameter θr was found to have the lowest accuracy. This was not surprising since most of the observed watercontents and pressure heads were in the relatively wet range (datanot further shown). In contrast, values of θs were very reliable.
The parameters α and n could be deﬁned reliably only for the top
layer. For all other layers, the tolerance intervals were very wide,
although the estimates themselves were acceptable. As indicatedby Hill (1998), this suggests that the information content of data
for the deeper layers was insuﬃcient for conclusive evaluation.
Since the top layer experienced more drying (Fig. 1), the simulations were found to be more sensitive to the retention shapeparameters α and n. By comparison, most or all of the data inthe deeper layers remained close to saturation, leading to unreliable estimates for α and n. The same is true for the hydraulicconductivity parameters Ksat and l. Still, the satisfactory agreement
between the observed and calculated data does indicate applicability of the Richards equation to the ﬂow processes at this site.
The HYDRUS-1D model (Simunek et al., 1998) was used
next to simulate ﬂow at the ﬁeld site using hydraulic parameters obtained with the PTFs. An example of the results of thesesimulations is shown in Fig. 2 for the 15-cm depth. The scatter in the simulated water contents was found to be substantial.
In general, the simulated soil water contents were lower thanthe observations and also decreased somewhat faster during
dry periods. Temporal variations in the soil water content were

for the model of Brooks and Corey (1964), and(Ksat
cm d−10.055 ± 0.006
50 ± 409
0.207 ± 0.300
0.167 ± 0.329
4.16 ± 7.64

Results

Eq. [12]. The particular equation was selected depending on the
invoked PTF. The unsaturated hydraulic conductivity (K) function was calculated using
K (h ) = K sat S e2/λ+l +2

n

[15]

and l is a pore-connectivity parameter. A separate simulation run
was performed with HYDRUS-1D for each PTF. Water retention
was estimated by applying the PTFs to the soil properties of eachindividual layer at the Bekkevoort ﬁeld site (0–25, 25–45, 45–65,
65–85, and 85–150 cm; see Table 1).
We evaluated the weighting methods in terms of their accuracy (i.e., errors in reproducing the training or hindcast data sets),and reliability (i.e., errors in reproducing the test data sets). The
training (i.e., determining the weights) was done with daily watercontents using moving windows that were from 30 to 150 d long.
All data outside the windows were used to test the multimodel
prediction. Average and standard deviations of RMSEs foreach window size were obtained separately for the testing andtraining sets.
The Richards water ﬂow model was calibrated using complete
time series of measured water contents, measured precipitation
data, and evaporation rates estimated from water budget computations (Pachepsky et al., 2006). The calibration was used toevaluate the applicability of the Richards equation to the observed
ﬂow processes, and to compare the accuracy of the calibrated
model with those of the multimodel predictions. The inverse
solution option of HYDDRUS-1D was used for this purpose.
The soil proﬁle was subdivided into ﬁve layers: 0 to 25, 25 to45, 45 to 65, 65 to 85, and 85 to 105 cm. The four parametersθr, θs, α, and n in Eq. [12] and the two additional parameters
Ksat and l in Eq. [14] were estimated for each layer by ﬁtting the
Richards ﬂow model to daily water content values averaged acrossthe trench at each depth. Initial estimates of the parameters wereobtained by using data on soil texture (Table 1) and the Rosetta
PTF included in the HYDRUS-1D software. A total of 30 parameters were estimated during model calibration.

F . 2. Daily average water contents at the 15-cm depth; lines are simula ons with 19 pedotransfer func ons; circles are observed values.
www.vadosezonejournal.org · Vol. 8, No. 1, February 20095

more pronounced with some PTFs than with others;
however, inspection of Fig. 2 shows that the diﬀerent
PTFs produced almost parallel time series that were
shifted in various degrees with respect to each other.
Computed correlation coeﬃcients between the timeseries supported this ﬁnding. A color-coded matrix ofcorrelation coeﬃcients between the calculated watercontent time series using diﬀerent pedotransfer functions
is shown in Fig. 3 for the 15-cm depth. Correlation coefﬁcients between simulated soil water contents exceeded
0.95 in >48% of the cases, while very few correlationcoeﬃcients were <0.5. Similar results were obtained forthe other depths (data not shown).
The accuracy and reliability of diﬀerent methods
to build the multimodel prediction are compared inFig. 4. The superensemble was found to be the worst
weighting method, characterized by high accuracy
and very low reliability (data not shown in Fig. 4
because they are outside the axis scales). Assigningequal weights (simple averaging) was the second worst
method. This method also had testing RMSE values
larger than the training RMSEs, although the diﬀerence was not as dramatic as with the superensemble
method. The accuracy of the superensemble was theworst of all the methods. Bayesian model averaging
fared only slightly better than simple averaging. Using
only the best model was a reasonably good approachin terms of accuracy. The reliability of this method,
however, was almost the same as with simple averaging
and Bayesian averaging. The results of the information
theory method were no diﬀerent from using the bestmodel. Using SVD in the superensemble appeared tobe the best method in that both the accuracy and reliability RMSE values were almost two times less thanthose of all other methods considered in this study. An
example in Fig. 5 of applying the SVD to the whole
250-d observation period as a training period illustrates the performance of this technique.
The results in Fig. 4 further show that training
with 30 d of data leads to lower reliability than when
60 or 90 d of data are used in the training. We founda weak trend of increasing accuracy with the duration of training. Figure 4 also compares multimodel
results with the results obtained with the Richards
model calibrated on data from the entire observation
period. The error bars for the calibrated model wereobtained from the statistical distributions of RMSE
for the same training and testing periods that wereused for the multimodel predictions. The data in Fig.
4 show that the accuracy of the superensemble with
SVD is comparable to the accuracy and reliability of
the calibrated Richards equation, even when weights
are derived from 30 d of observations. The overall
accuracy and reliability of the SVD multimodel was
slightly less, or comparable to, the accuracy of thecalibrated model.
The reliability of the SVD method was foundto be not random in time. The RMSE time series
shown in Fig. 6 indicate that the SVD multimodel
www.vadosezonejournal.org · Vol. 8, No. 1, February 2009

F . 3. Spearman correla ons between me series of water contents at the 15-cm
depth from simula ons obtained with diﬀerent pedotransfer func ons (PTFs).

F . 4. Root mean squared errors of mul model simula ons of the training andtes ng data sets; results are for the best in training model (red), the arithme caverage of all individual models (green), the superensemble with singular-value
decomposi on (pink), Bayesian averaging (blue), the use of the informa on theory (cyan), and calibrated Richards equa on (brown).
6F . 6. Dependencies of root mean square errors in the training(solid line) and tes ng (dashed line) data sets on the beginningand dura on of the training period for the superensemble withsingular-value decomposi on method for training periods of (a) 30,(b) 60, and (c) 90 d.

F . 5. Results of applying the superensemble with the singularvalue decomposi on to mul modeling of the whole observa on
period. Daily average water contents at 15-cm depth; lines are
simula ons and circles are observed values.

was slightly less accurate, and slightly more reliable, in the fall
toward the end of the experiments. This was observed for all
depths (data partly shown in Fig. 6). Spearman correlation coeﬃcients between SVD-simulated and measured soil water contents
were between 0.6 and 0.8 at the 15-cm depth, between 0.5 and
0.8 at the 35-cm depth, and between 0.4 and 0.8 depending on
the starting day of the training period.
Time series of weights of the individual models as a functionof the beginning of the training period are shown in Fig. 7 for theﬁrst nine of 19 individual PTFs at the 15-cm depth. Inspectionof the weights shows that none of the models was consistently
ignored. Depending on the training period, the same PTF could
have larger or smaller input into the multimodel prediction.
While the weights of a few PTFs were sometimes negative, inmost cases the PTF predictions were summed up with positive
weights. Similar results were obtained for the other four depths(data not shown).

F . 7. Weights of the simula ons with the first nine pedotransfer
func ons in the superensemble mul model with singular-value decomposi on obtained for the 15-cm depth with 30-d training periods.

was found to be by far the best method in our study (Fig. 4). Thisconclusion cannot be generalized. We believe that a preliminary
assessment of all methods should be performed for speciﬁc applications, at least until some experience is accumulated in using
multimodels for simulations of ﬂow in variably saturated soils.
Also, the saturated hydraulic conductivity, Ksat, was estimated froma single PTF in this study, whereas two more PTFs (Wösten et al.,
1999; Schaap et al., 2001) also could be used to expand the multimodel, which includes 57 rather than 19 individual models.
Values of weights in Eq. [1] were negative for some models
when the regression-based methods were applied. Bayesianmodel averaging and Akaike criteria application do not allow
negative weights. The regression-based methods treat weights asﬁtting coeﬃcients that can have any sign. The reason for negative weights is that the model predicts behavior opposite to thatmeasured. That may happen because of prediction errors of someindividual PTFs. Consider a PTF that predicts water retention
that is too low. Then the ﬂow model will predict the much fastermovement of an inﬁltration pulse through the soil proﬁle thanshould happen in reality. In such a case, a water loss will be predicted at some depths for some period of redistribution of theDiscussion
Multimodel simulations rely on the assumption that thephysics of each individual model is correct. Unfortunately, PTFs
are merely estimators of model parameters, and as such do notaﬀect the description of the physical processes involved, in thiscase the Richards equation based on mass conservation in combination with the Darcy–Buckingham law for the ﬂuid ﬂux. Thisﬂow model was applicable in principle since it could be ﬁtted reasonably well to the data (Fig. 1). Our attempts to use multimodel
simulations based on the Richards equation are hence warranted.
The similarity in simulated physical processes caused strong correlations between predictions from individual models (Fig. 2 and3), which is a typical feature of multimodel predictions.
Diﬀerent methods used for the multimodel predictions were
found to perform quite diﬀerently. The superensemble with SVD
www.vadosezonejournal.org · Vol. 8, No. 1, February 20097

inﬁltration pulse while in actuality water will still be accumulating in these horizons. In such a case, the model with this PTF
will indeed predict behavior opposite to that observed, and thevalues of weights from regression-based multimodeling may be
negative for this model. Note that the SVD method, which allows
negative weights, performed the best.
The list of methods to build the multimodel ensemble inthis work is far from exhaustive. For example, principal component analysis (Young, 2002) can be used to replace correlated
predictors Fi in Eq. [1] with a smaller number of uncorrelated
predictors. The number of independent variables in Eq. [1] canbe decreased by using only a subset of models (Kharin and Zwiers,
2002). Various statistics have been proposed to select the “best”
subset of input variables (Neter and Wasserman, 1974; Kharin
and Zwiers, 2002; Regonda et al., 2006). Young (2002) andRegonda et al. (2006) listed several other methods to combine
models. Recent developments include using patterns found inexperimental data to adjust the weights of individual forecasts in
multimodel predictions (Zheng et al., 2004). There is no reason
why a nonlinear combination of individual simulations could not
be used in multimodels. Using data mining techniques, such as
artiﬁcial neural networks, may be beneﬁcial toward that end.
Both the accuracy and, to a lesser extent, the reliability of the
SVD multimodel had temporal variations (Fig. 6). The accuracy
but not the reliability of the SVD multimodel was worse when
the training occurred between Days 135 and 185. Juxtaposing
data in Fig. 1 and Fig. 6 shows that the accuracy decreased fortraining periods that included a period of intensive drying followed by extensive wetting. This was not correctly reﬂected bysimulations based on the Richards ﬂow model. All individualmodels performed somewhat poorly in a qualitative sense during
these time periods, and training was less successful than when the
Richards equation was a satisfactory model. One possible remedy
would be using diﬀerent weighting coeﬃcients by introducing
an indicator variable for prevailing antecedent soil conditions orboundary conditions. The viability of such an approach should
depend on the magnitude of the available data sets. This may
present an interesting avenue for future research.
The relatively poor performance of the superensemble
method was somewhat expected based on previous experience
with climate and hydrologic models (i.e., Regonda et al., 2006).
Using the best model was the second best method in this work.
The superiority of a multimodel prediction to that of a single
model has been observed consistently in meteorological forecasts.
This is attributed to the fact that it is diﬃcult or impossible todeﬁne the single best model (Hagedorn et al., 2005) if various
training and testing periods are considered. Figure 8 shows that
there was no generally superior model in our study, and thatthe model considered the best depends on when the training
period started. Hagedorn et al. (2005) noted that a single model
may perform better in some situations but that the multimodel
approach in the long term will give more reliable simulations.
The premise of superiority of multiple-source simulation
systems is based on the hypothesis that “two or more inaccurate
but independent predictions of the same future events may be
combined in a very speciﬁc way to yield predictions that are on
average more accurate than either of any of them taken individually” (Thompson, 1977). This hypothesis is supported by some
artiﬁcial neural networks studies, speciﬁcally the group method
www.vadosezonejournal.org · Vol. 8, No. 1, February 2009

F . 8. Pedotransfer func ons (PTFs) that provided the most accurate simula ons at four depths using a 30-d training period for
diﬀerent days of the beginning of training.

of data handling (e.g., Pachepsky and Rawls, 1999), which uses
combinations of individually poor predictions to create predictors with excellent performance. The results of this study seem tofurther corroborate this hypothesis.
The weights of modeling results with individual PTFs
depended on the training period (Fig. 7). This is an expected
outcome from the regression analysis when diﬀerent regression
coeﬃcients are obtained if diﬀerent data subsets are used todevelop the regression. Similarly, the training and testing accuracies were also dependent on the training period (Fig. 4). In spiteof diﬀerences in weight values, however, the overall accuracy ofthe multimodel depended much more on the method used to
combine predictions than on the testing data set (Fig. 4). The
selection of the training data set aﬀected the participation ofindividual models in the multimodeling results, but the overall
accuracy and reliability of the multimodel did not vary much. In
particular, the accuracy and reliability of the SVD multimodel
remained relatively high for all training data sets (Fig. 4).
We include here a comment about terminology. The term
ensemble prediction is sometimes also used to refer to the integration of the results of many simulation runs but for cases where
the initial conditions and parameters or models are randomly
perturbed. The term multimodel ensemble prediction is used whenboth diﬀerent models and their input perturbations are used inindividual simulation runs.
In this study, we deliberately made no attempt to discriminate between PTFs or to single out the best PTF or PTFs. First,
the multimodel concept does not presume any such analysis since
it is agreed from the outset that all models will contribute. Second,
based on the data in Fig. 8, the result would depend on the
training period. We note that several attempts have been madeto relate the performance of PTFs to similarity in geographic
regions, to the size of the database used to develop a PTF, to the8Heimovaara, T.J., and W. Bouten. 1990. A computer controlled 36-channel
time domain reﬂectometry system for monitoring soil water contents. Water Resour. Res. 26:2311–2316.
Hill, C.M. 1998. Methods and guidelines for eﬀective model calibration. Water
Resour. Invest. Rep. 98–4005. USGS, Denver, CO.
Jacques, D. 2000. Analysis of water ﬂow and solute transport at the ﬁeld scale.
Ph.D. diss. 454. Faculteit Landbouwkundige en Toegepaste Wetenschappen, Katholieke Univ. Leuven, Leuven, Belgium.
Jacques, D., J. Simunek, A. Timmerman, and J. Feyen. 2001. Calibration of the
Richards’ and convection–dispersion equations to ﬁeld-scale water ﬂow
and solute transport under rainfall conditions. J. Hydrol. 259:15–31.
Jolliﬀe, I.N., and D.B. Stephenson. 2003. Forecast veriﬁcation: A practitioner’s
guide in atmospheric science. John Wiley & Sons, New York.
Kharin, V., and F.W. Zwiers. 2002. Climate predictions with multimodel ensembles. J. Clim. 15:793–799.
Krishnamurti, T.N., C.M. Kishitawal, Z. Zhang, T. Larow, D. Bachiochi, and E.
Williford. 2000. Multimodel ensemble forecasts for weather and seasonal
climate. J. Clim. 13:4196–4216.
Mayr, T., and N.J. Jarvis. 1999. Pedotransfer functions to estimate soil water retention parameters for a modiﬁed Brooks–Corey type model. Geoderma 91:1–9.
McLachlan, G.J., and T. Krishnan. 1997. The EM algorithm and extensions.
John Wiley & Sons, New York.
Neter, J., and W. Wasserman. 1974. Applied linear statistical models: Regression,
analysis of variance, and experimental designs. Richard D. Irwin, Homewood, IL.
Pachepsky, Ya.A., A.K. Guber, M.Th. van Genuchten, T.J. Nicholson, R.E. Cady,
J. Simunek, and M.G. Schaap. 2006. Model abstraction techniques forsoil water ﬂow and transport. NUREG/CR-6884. U.S. Nucl. Regulatory
Commiss., Washington, DC.
Pachepsky, Ya.A., and W.J. Rawls. 1999. Accuracy and reliability of pedotransfer
functions as aﬀected by grouping soils. Soil Sci. Soc. Am. J. 63:1748–1757.
Pachepsky, Ya.A., and W.J. Rawls. 2004. Status of pedotransfer functions. p. vii–
xviii. In Ya.A. Pachepsky and W. J. Rawls (ed.) Development of pedotransfer functions in soil hydrology. Elsevier, Amsterdam.
Palmer, T.N., A. Alessandri, U. Andersen, P. Cantelaube, M. Davey, P. Délécluse,
et al. 2004. Development of a European multimodel ensemble system forseasonal-to-interannual prediction (DEMETER). Bull Am. Meteorol. Soc
85:853–872.
Palmer, T.N., C. Brankovic, and D.S. Richardson. 2000. A probability anddecision-model analysis of PROVOST seasonal multi-model ensemble
integrations. Q. J. R. Meteorol. Soc. 126:2013–2035.
Pavan, V., and F.J. Doblas-Reyes. 2000. Multimodel seasonal hindcasts over theEuro-Atlantic: Skill scores and dynamic features. Clim. Dyn. 16:611–625.
Peng, P., A. Kumar, H. van den Dool, and A.G. Barnston. 2002. An analysis ofmultimodel ensemble predictions for seasonal climate anomalies. J. Geophys. Res. 107(D23):4710, doi:10.1029/2002JD002712.
Petersen, G.W., R.L. Cunningham, and R.P. Matelski. 1968. Moisture characteristics of Pennsylvania soils: I. Moisture retention as related to texture.
Soil Sci. Soc. Am. Proc. 32:271–275.
Poeter, E., and D. Anderson. 2005. Multimodel ranking and inference in groundwater modeling. Ground Water 43:597–605.
Press, W.H., B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling. 1992. Numerical recipes in FORTRAN 77: The art of scientiﬁc computing. Cambridge
Univ. Press, Cambridge, UK.
Raftery, A.E., R. Balabdaoui, T. Gneiting, and M. Polakowski. 2003. Using
Bayesian model averaging to calibrate forecast ensembles. Tech. Rep. 440.
Dep. of Statistics, Univ. of Washington, Seattle.
Rajagopalan, B., U. Lall, and S.E. Zebiak. 2002. Categorical climate forecasts
through regularization and optimal combination of multiple GCM ensembles. Mon. Weather Rev. 130:1792–1811.
Rajkai, K., and Gy. Várallyay. 1992. Estimating soil water retention from simpler
properties by regression techniques. p. 417–426. In M.Th. van Genuchten
et al. (ed.) Proc. Int. Worksh. on Indirect Methods for Estimating the Hydraulic Properties of Unsaturated Soils, Riverside, CA. 11–13 Oct. 1989.
U.S. Salinity Lab., Riverside, CA.
Rawls, W.J., and D.L. Brakensiek. 1985. Prediction of soil water properties forhydrologic modeling. p. 293–299. In E.B. Jones and T.J. Ward (ed.) Proc.
Symp. Watershed Management in the Eighties, Denver, CO. 30 Apr.–1
May 1985. Am. Soc. Civil Eng., New York.
Rawls, W.J., D.L. Brakensiek, and B. Soni. 1983. Agricultural management ef-

homogeneity of this database in terms of measurement methods, or other factors (Pachepsky and Rawls, 2004). All of theseattempts appear to have been inconclusive. Our results by nomeans make a general statement about the quality or reliabilityof an individual PTF because the results depend on the chosendata set and training period.
The similarity in accuracy between the calibrated Richards
model for the layered soil and the multimodel predictions with
PTFs suggest the interesting possibility of using soil moisturemonitoring data to obtain a more accurate predictive soil waterﬂow model. Monitoring soil moisture and ﬁnding weights forthe multimodel prediction may well be a very viable approach tosimulating ﬁeld water ﬂow in the vadose zone.
A
The work has been partially supported through the Interagency
Agreement RES-02-008, “Model Abstraction Techniques for Soil Water
Flow and Transport.”

References
Armstrong, J.S. 2001. Combining forecasts. p. 417–439. In J.S. Armstrong (ed.)
Principles of forecasting: A handbook for researchers and practitioners.
Kluwer Acad. Publ., Dordrecht, the Netherlands.
Barnston, A.G., S. Mason, L. Goddard, D.G. Dewitt, and S.E. Zebiak. 2003.
Multimodel ensembling in seasonal climate forecasting at IRI. Bull. Am.
Meteorol. Soc. 84:1783–1796.
Bates, J.M., and C.W.J. Granger. 1969. The combination of forecasts. Oper. Res.
Q. 20:451–468.
Baumer, O.M. 1992. Predicting unsaturated hydraulic parameters. p. 341–354.
In M.Th. van Genuchten et al. (ed.) Proc. Int. Worksh. on Indirect Methods for Estimating the Hydraulic Properties of Unsaturated Soils, Riverside, CA. 11–13 Oct. 1989. U.S. Salinity Lab., Riverside, CA.
Brooks, R.H., and A.J. Corey. 1964. Hydraulic properties of porous media. Hydrol. Pap. 3. Colorado State Univ., Fort Collins.
Burnham, K.P., and D.R. Anderson. 2002. Model selection and inference: A
practical information-theoretic approach. Springer-Verlag, New York.
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.) Proc. Int. Worksh. on Indirect Methods
for Estimating the Hydraulic Properties of Unsaturated Soils, Riverside,
CA. 11–13 Oct. 1989. U.S. Salinity Lab., Riverside, CA.
Canarache, A. 1993. Physical-technological maps: A possible product of soil survey for direct use in agriculture. Soil Technol. 6:3–16.
Clemen, R.T. 1989. Combining forecasts: A review and annotated bibliography.
Int. J. Forecast. 5:559–583.
Doblas-Reyes, F.J., M. Déqué, and J.-Ph. Piedelievre. 2000. Model and multimodel spread and probabilistic seasonal forecasts in PROVOST. Q. J. R.
Meteorol. Soc. 126:2069–2088.
Doblas-Reyes, F.J., R. Hagedorn, and T.N. Palmer. 2005. The rationale behindthe success of multi-model ensembles in seasonal forecasting. II: Calibration and combination. Tellus Ser. A 57:234–252.
FAO. 1975. Soil map of the world at 1:5,000,000. Vol. 1. Europe. FAO, Rome.
Graham, R.J., A.D.L. Evans, K.R. Mylne, M.S.J. Harrison, and K.B. Robertson.
2000. An assessment of seasonal predictability using atmospheric generalcirculation models. Q. J. R. Meteorol. Soc. 126:2211–2240.
Guber, A.K., Ya.A. Pachepsky, M.Th. van Genuchten, W.J. Rawls, D. Jacques,
J. Simunek, R.E. Cady, and T.J. Nicholson. 2006. Field-scale water ﬂow
simulations using ensembles of pedotransfer functions for soil water retention. Vadose Zone J. 5:234–247.
Gupta, S.C., and W.E. Larson. 1979. Estimating soil water retention characteristics from particle-size distribution, organic matter percent, and bulk
density. Water Resour. Res. 15:1633–1635.
Hagedorn, R., F.J. Doblas-Reyes, and T.N. Palmer. 2005. The rationale behindthe success of multi-model ensembles in seasonal forecasting: I. Basic concept. Tellus Ser. A 57:219–233.
Hall, D.G.M., M.J. Reeve, A.I. Thomasson, and V.F. Wright. 1977. Water retention, porosity and density of ﬁeld soils. Soil Surv. Tech. Monogr. 9.
Rothamsted Exp. Stn., Harpenden, UK.

www.vadosezonejournal.org · Vol. 8, No. 1, February 20099

fects on soil water processes. I: Soil water retention and Green–Ampt parameters. Trans. ASAE 26:1747–1752.
Rawls, W.J., D.L. Brakensiek, and K.E. Saxton. 1982. Estimation of soil water
properties. Trans. ASAE 25:1316–1320.
Rawls, W.J., D. Giménez, and R. Grossman. 1998. Use of soil texture, bulk density and slope of the water retention curve to predict saturated hydraulicconductivity. Trans. ASAE 41:983–988.
Regonda, S.K., B. Rajagopalan, M. Clark, and E. Zagona. 2006. A multimodel ensemble forecast framework: Application to spring seasonal
ﬂows in the Gunnison River Basin. Water Resour. Res. 42:W09404,
doi:10.1029/2005WR004653.
Saxton, K.E., W.J. Rawls, J.S. Romberger, and R.I. Papendick. 1986. Estimating generalized soil-water characteristics from texture. Soil Sci. Soc. Am. J.
50:1031–1036.
Schaap, M.G. 2004. Graphic user interfaces for pedotransfer functions. p. 349–
356. In Ya.A. Pachepsky and W.J. Rawls (ed.) Development of pedotransfer functions in soil hydrology. Elsevier, Amsterdam.
Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. Rosetta: A computer
program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol. 251:163–176.
Shukla, J., J. Anderson, D. Baumhefner, C. Brankovic, Y. Chang, E. Kalnay,
et al. 2000. Dynamical seasonal prediction. Bull. Am. Meteorol. Soc.
81:2593–2606.
Simunek, J., M. Sejna, and M.Th. van Genuchten. 1998. The HYDRUS-1D
software package for simulating the one-dimensional movement of water,
heat, and multiple solutes in variably-saturated media. Version 2.0. IGWMC- TPS- 70. Int. Ground Water Modeling Ctr., Colorado School ofMines, Golden.
Soil Survey Staﬀ. 1999. Soil Taxonomy: A basic system of soil classiﬁcation formaking and interpreting soil surveys. 2nd ed. U.S. Gov. Print. Oﬃce,
Washington, DC.
Thompson, P.D. 1977. How to improve accuracy by combining independent
forecasts. Mon. Weather Rev. 105:228–239.
Tomasella, J., and M.G. Hodnett. 1998. Estimating soil water retention characteristics from limited data in Brazilian Amazonia. Soil Sci. 163:190–202.
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.
Varallyay G., K. Rajkai, Ya.A. Pachepsky, and R.A. Shcherbakov. 1982. Mathematical description of soil water retention curve. (In Russian.) Pochvovedenie 4:77–89.
Vereecken, H., J. Maes, J. Feyen, and P. Darius. 1989. Estimating the soil moisture retention characteristics from texture, bulk density and carbon content. Soil Sci. 148:389–403.
Williams, J., P. Ross, and K. Bristow. 1992. Prediction of the Campbell water
retention function from texture, structure, and organic matter. p. 427–442.
In M.Th. van Genuchten et al. (ed.) Proc. Int. Worksh. on Indirect Methods for Estimating the Hydraulic Properties of Unsaturated Soils, Riverside, CA. 11–13 Oct. 1989. U.S. Salinity Lab., Riverside, CA.
Wösten, J.H.M., A. Lilly, A. Nemes, and C. Le Bas. 1999. Development anduse of a database of hydraulic properties of European soils. Geoderma
90:169–185.
Ye, M., S.P. Neuman, and P.D. Meyer. 2004. Maximum likelihood Bayesian
averaging of spatial variability models in unsaturated fractured tuﬀ. Water
Resour. Res. 40:W05113, doi:10.1029/2003WR002557.
Young, G. 2002. Combining forecasts for superior prediction. p. 107–111. In
Conf. on Probability and Statistics in the Atmos. Sci., 16th, Orlando, FL.
12–17 Jan. 2002. Am. Meteorol. Soc., Boston, MA.
Yun, W.T., L. Stefanova, and T.N. Krishnamurti. 2003. Improvement of themultimodel superensemble technique for seasonal forecasts. J. Clim.
16:3834–3840.
Zheng, X., M. Sugi, and C.S. Frederiksen. 2004. Interannual variability andpredictability in an ensemble of climate simulations with the MRI-JMA
AGCM. J. Meteorol. Soc. Jpn. 82:1–18.

www.vadosezonejournal.org · Vol. 8, No. 1, February 200910