SOIL SCIENCE SOCIETY OF AMERICAJOURNAL
JULY-AUGUST 1998

VOL. 62

No. 4

DIVISION S-l—SOIL PHYSICS
Neural Network Analysis for Hierarchical Prediction of Soil Hydraulic Properties
Marcel G. Schaap,* Feike J. Leij, and Martinus Th. van Genuchten
ABSTRACT

is often time consuming and expensive, while the results
may not be accurate. An alternative is the use of pedotransfer functions (PTFs), which estimate the hydraulicproperties through correlation with more easily measured or widely available soil parameters (Bouma andvan Lanen, 1987; van Genuchten and Leij, 1992).
A variety of PTFs with different mathematical concepts, predicted properties, and input data requirements
have been developed. Quasi-physical methods by Aryaand Paris (1981), Haverkamp and Parlange (1986), andTyler and Wheatcraft (1989) use the concept of shapesimilarity between pore- and particle-size distributions.
The vast majority of PTFs, however, are empirically
based on relatively simple linear regression equations.
Although considerable differences exist among PTFs interms of the required input data, all of them use at least
some information about the particle-size distribution.
When only the textural classification is known, simple"class" PTFs can be used to provide average hydraulicproperties for each soil textural class (Carsel and Parrish, 1988; Wosten et al., 1995). When the actual particlesize distribution is known, PTFs that predict continuously changing hydraulic properties across the textural
triangle can be used.
Pedotransfer function predictions may be improvedby extending the input data through addition of basic
soil properties like bulk density, porosity, or organic
matter content (Rawls and Brakensiek, 1985; Vereecken et al., 1989). Additional improvements may beachieved by including one or more water retention data
points (Rawls et al., 1992; Williams et al., 1992). Ahuja
et al. (1989) and Messing (1989) similarly improved predictions of saturated hydraulic conductivity, Ks, by using
effective porosity data, which they defined as the total
porosity minus the water content at 10 or 33 kPa pres-

The solution of many field-scale flow and transport problems requires estimates of unsaturated soil hydraulic properties. The objective of this study was to calibrate neural network models for predictionof water retention parameters and saturated hydraulic conductivity,
A's, from basic soil properties. Twelve neural network models weredeveloped to predict water retention parameters using a data set of1209 samples containing sand, silt, and clay contents, bulk density,
porosity, gravel content, and soil horizon as well as water retention
data. A subset of 620 samples was used to develop 19 neural networkmodels to predict K,. Prediction of water retention parameters andA', generally improved if more input data were used. In a more detailed
investigation, four models with the following levels of input data wereselected: (!) soil textural class, (ii) sand, silt, and clay contents, (iii)sand, silt, and clay contents and bulk density, and (iv) the previous
variables and water content at a pressure head of 33 kPa. For water
retention, the root mean square residuals decreased from 0.107 forthe first to 0.060 m3 nr' for the fourth model while the root meansquare residual K, decreased from 0.627 to 0.451 log(cm d '). The
neural network models performed better on our data set than four
published pedotransfer functions for water retention (by ~ 0.01-0.05
m3 m ') and better than six published functions for K, (by -0.1-0.9
order of magnitude). Use of the developed hierarchical neural networkmodels is attractive because of improved accuracy and because it
permits a considerable degree of flexibility toward available input
data.

C

ONCERN ABOUT THE QUALITY of soil and water resources has motivated the development of increasingly sophisticated models for predicting water flow andsolute transport in unsaturated soils. These models generally require knowledge of the soil water retention,0(/z), and unsaturated hydraulic conductivity, K(h),where 0 is the volumetric water content and h is thepressure head. Direct measurement of these properties

Abbreviations: PTF, pedotransfer function; pb, bulk density; EP10,
effective porosity at 10 kPa; EP33, effective porosity at 33 kPa; EPlsgt>,
effective porosity at 1500 kPa; HOR, horizon; POR, porosity; PTF,
pedotransfer function; RMSR, root mean square residual, Eq. [4] and[5]; SSC, sand, silt, clay content; TXT, textural class.

USDA-ARS, U.S. Salinity Lab., 450 West Big Springs Rd., Riverside,
CA 92507-4617. Received 2 June 1997. *Corresponding author(mschaap@ussl.ars.usda.gov).
Published in Soil Sci. Soc. Am. J. 62:847-855 (1998).
847

848

SOIL SCI. SOC. AM. J., VOL. 62, JULY-AUGUST 1998

sure head. Other researchers have predicted soil hydraulic properties using more limited or extended sets of
input variables (Rawls et al., 1992; Schaap and Bouten,
1996; Vereecken et al., 1989, 1990). Such hierarchicalapproaches are of great practical use since they permit
more flexibility toward the required input data when
predicting the hydraulic properties.
Neural-network-based PTFs were recently used byPachepsky et al. (1996), Schaap and Bouten (1996), andTamari et al. (1996). The feed-forward backpropagation
or radial basis functions used by these researchers areable to approximate any continuous (nonlinear) function to any desired degree of accuracy (Hecht-Nielsen,
1990; Haykin, 1994). An advantage of neural networks,compared with traditional PTFs, is that neural networks
require no a priori model concept. The optimal relations
that link input data (basic soil properties) to output data
(hydraulic parameters) are obtained and implemented
in an iterative calibration procedure. For a data setof 204 sandy soil samples, Schaap and Bouten (1996)
showed that neural networks made predictions with significantly smaller errors than more traditional linear
regression approaches. Pachepsky et al. (1996) used adata set of 230 soil samples. They found that neuralnetworks predicted water retention points better thanmultilinear regression, but that the two methods produced comparable results when water retention parameters were predicted. Tamari et al. (1996) used synthetic
K(h) data sets and found that neural networks were notbetter than multilinear regression models if the uncertainty in the data was large. The neural networks performed better than regression when real soil data
were used.
It was the objective of this study to predict van Genuchten (1980) water retention parameters and saturated
hydraulic conductivity for a hierarchical input structure.
The corresponding sequence of neural network PTFs
yields a collection of models of which the most suitable
one can be selected depending on the available basic

soil properties. We also used neural networks to establish which basic soil properties are the most relevant
for predicting the hydraulic properties. Contrary to previous work on neural network PTFs, a data set was
employed that covers most of the textural triangle, thus
ensuring broad applicability of the hierarchical approach. Conventional error criteria and independent
data were used to evaluate the performance of the hierarchical neural network models, and to test performanceof previously published PTFs. The uncertainty in thepredicted soil hydraulic properties was evaluated by
combining the neural network approach with the bootstrap method (Efron and Tibshirani, 1993).
MATERIALS AND METHODS
Data Set
The data set for calibration and testing was extracted froma database consisting of 4515 laboratory samples taken fromabout 30 sources in the USA (W.J. Rawls, USDA-ARS Hydrology Lab., Beltsville, MD, 1996, personal communication).
We selected 1209 samples for water retention (containing six
to 13 points with different pressure head ranges), soil texture,and bulk density; a subset of 620 samples had Ks measurements. Figures la and Ib show the distributions of the samples
in the USDA textural triangle for the entire 1209-sample data
set and the /Cs subset. For subsequent analysis, the samples inthe sandy clay, clay, and silty clay classes were grouped into
one class ("clays"); silt and silt loam classes were grouped
into "silts"; 405 samples were from A horizons, 588 from B
horizons, and 216 from C horizons.
Water retention data for each of the 1209 samples werefitted to the van Genuchten (1980) equation:e. - er

100 A 0

80

100
20

80

40

60

Clay [%]

Silt [%]

Clay [%]

60

100

80

60

[1]

[1 + (ah)"]"where 6S and 6r are1 saturated and residual water contents
respectively; a (cm" ) and n are curve shape parameters, andfollowing van Genuchten (1980), m = 1 — lln. Fitting was
carried out with the simplex or amoeba algorithm (Nelder

40

Silt [%]
40

100

100

80

60

40

20

100

Sand[%]
Sand[%]
Fig. 1. Textural distribution of the (a) water retention data set and the (b) saturated hydraulic conductivity (Ks) subset. S: sand, C: clay, Si: silt,

L: loam, s: sandy, c: clayey, si: silty, 1: loamy.

SCHAAP ET AL.: NEURAL NETWORK ANALYSIS FOR HYDRAULIC PROPERTY PREDICTION

and Mead, 1965; Press et al., 1988) with the following constraints: 0.0 < 6r < 0.3 m3 m~3, 0.6(j>< 6S < <(> m3 m~ 3 (where
4> is the total porosity), 0.0001 < a < 1.000 cirr1, and 1.001< « < 10. The parameters a and n were then log-transformed
to obtain approximately normal distributions. The same transformation was also carried out for Ks measurements.

Published Pedotransfer Functions
In this study we used four previously published PTFs topredict water retention parameters (Rl through R4 in Table
1) and six published PTFs to predict saturated hydraulic conductivity (Kl through K6 in Table 2). Tables 1 and 2 also
show the required input data for each PTF. All PTFs consistof the same coefficients and variables as in the original publications. We did not recalibrate the PTFs for our data set becauseof their empirical nature. Recalibration would not only require
adaption of the coefficients but also an evaluation of whether
the input variables or expressions used are actually appropriate for the current data set. Recalibration might have resulted in completely different PTFs than were published.
The Rl and R3 functions predicted parameters of the
Brooks and Corey (1964) water retention function, given by:

6(70 = er,BC + (es,BC - er,BC)[2]

where hb is the air-entry pressure and X is an empirical parameter. The function Rl (Rawls and Brakensiek, 1985) assumed
that the saturated water content, 6sBC, is equal to the porosity.
The function R3 (Cosby et al., 1984) predicted Brooks-Coreyparameters according to Campbell (1974), i.e., with 6 as theindependent variable and 6rjBC = 0.
The R2 and R4 functions predicted parameters of Eq. [1].
In R2, Rawls and Brakensiek (1985) transformed the BrooksCorey parameters as predicted by Rl into the van Genuchten(1980) parameters using the expressions in Table 1. For R4,
m in Eq. [1] is set to 1 (Vereecken et al., 1989).
Functions Kl, K2, and K4 (Table 2) provide relatively
straightforward expressions between K,, and the listed soil
properties. Function K3 of Saxton et al. (1986) predicts a K(Q)
Table 1. Input and output data of four pedotransfer functions(PTF) that predict selected soil water retention parameters.
PTF
Input
OutputEq. [2]Rl Rawls andBrakensiek (1985)
Sand, clay, <(>0,3C

Sand, clay, <)>Sand, clay, <t>
R2 Rawls andBrakensiek (1985)

Eq. [1], m = 1 - 1/ne, = e,Bc
a = I/Atn = X+ 1

R3 Cosby et al.
(1984)
Sand
Sand
Clay
R4 Vereecken
et al. (1989)

NA
NA
NA
NA

e,
e,

NA
NA
NA
NA

O^BC
h\,
\
Eq. [1], m = 1

0.771

a
n
Eq. [2], e,,BC = 0

Clay, pb
6,
Clay, Organic C Os
Sand, clay, Pb, aOrganic C
Sand, clay
n

0.809
0.966
0.703
0.848
0.680
0.560

t Published coefficient of determination for calibration; NA = notavailable.

849

relationship using sand and clay contents. Substitution of thetotal porosity, 4>, for 9 yields an estimate for Ks. The coefficientof determination of this PTF is very high (R2 = 0.95, Saxton
et al., 1986) since it was fitted to predictions of a PTF reported
earlier by Rawls et al. (1982). Function K5 by Ahuja et al.
(1989) expresses Ks as a power function of effective porosity
parameter (EP33, calculated from the total porosity minus thewater content, 633, at 33 kPa). As parameter values in thispower function we used 1015 for the multiplication factor and
4.0 for the exponent as suggested by Rawls et al. (1992).
Finally, K6 is a quasi-physical PTF by Mishra and Parker(1989, 1992) that calculates Ks based on closed-form expressions by Mualem (1976), Brutsaert (1968), and van Genuchten (1980).

Neural Network Models
Since several textbooks on neural networks are available(e.g., Hecht-Nielsen, 1990; Haykin, 1994), we provide only abrief summary of the neural network approach. In this study
we used a back-propagation neural network with one hidden
layer. This type of neural network is a nonlinear data transformation structure consisting of input and output nodes connected to a number of hidden nodes by adaptable coefficients.
The number of input and output nodes corresponds to thenumber of input and output variables. The number of hidden
nodes depends on the complexity of the underlying problem
and is determined empirically by calibrating neural networks
with different numbers of hidden nodes. After Schaap andBouten (1996), we used six hidden nodes. Both the hidden
and output nodes contain sigmoidal transfer functions thatprovide the neural network with nonlinear capabilities. The
coefficients were obtained in an iterative calibration procedure based on, in our case, the Levenberg-Marquardt algorithm (Marquardt, 1963). This algorithm minimizes the objective function:0(t,t') = 2

(t,,-

[3]

i=l 7=1

where t and t' are the measured and predicted output variables,
NZ is the total number of calibration samples, JV, is the total
number of hydraulic parameters (one for Ks and four for theretention parameters in Eq. [1]).
A neural network model obtained in the calibration procedure should always be tested on independent data. Usually
this is done by calibrating the neural network on one part ofa data set and subsequently testing the network on the otherpart. We followed this approach, and also estimated uncertainty in neural network predictions by combining multiple
Table 2. Input and output data of six pedotransfer functions(PTFs) that predict the saturated hydraulic conductivity, Ks.
Input data
OutputPTF
0.839
Kl Cosby et al.
K,
Sand(1984)NAK,
K2 Brakensiek et al.
Sand, clay, <h(1984)
K3 Saxton et al.
Sand, clay, 00.95(1986)K4 Vereecken et al.
Sand, clay, pi,. Organic C
K,
0.20(1990)K5 Ahuja et al.
Effective porosity at 33 kPa K,
0.33-0.71§(1989)K6 Mishra and Parker ()„ 6,, aK,NA(1991)____________________________________t Published coefficient of determination for calibration; NA = notavailable.
± See text for explanation.
§ Ahuja et al., 1989.

850

SOIL SCI. SOC. AM. J., VOL. 62, JULY-AUGUST 1998

calibration and testing runs with the bootstrap method (Efron
and Tibshirani, 1993).
The bootstrap method is a nonparametric technique forcalculating the probability distribution of any statistic, in ourcase hydraulic parameters as predicted by neural networks.
The bootstrap does not make assumptions about the shapeof the distribution of the underlying population, nor does itmatter how complex is the calculation of the statistic (Efron
and Tibshirani, 1993). It is assumed that the data set involved
is a good representation of the original population with theunderstanding that this data set is only one particular realization of that population. Calibration of a model on anotherrealization would therefore always lead to slightly different
predictions of 9r, 0S, a, n, and Ks and thus cause uncertainty
about the true value of these parameters. Unfortunately, we
often have only one realization of the population. Bootstrapping simulates different realizations by repeated random resampling with replacement of the original data set of size N
to yield B bootstrap data sets of size N. Each bootstrap data
set contains somewhat different data, which results in B neuralnetwork models, all of which may differ slightly. Efron andTibshirani (1993) suggested B to be between 50 and 200; weused B = 60.
Because resampling is done with replacement, each sample
has a chance of 1 - [(TV - 1)1N]N (approximately 63%) to be
selected once or multiple times for a bootstrap data set. Eachbootstrap data set thus contains one or more copies of 63%
of the samples. Although formally not part of the bootstrap
method, the remaining 37% of the samples can be used tocarry out an independent test (cf. Efron and Tibshirani, 1993,
p. 252-254). By calculating averages and standard deviationsof B testing results, one obtains robust values of the predicted
hydraulic parameters and associated uncertainty estimates forindependent data.
The combined neural network-bootstrapping analysis was
carried out with a slightly adapted TRAINLM routine of theneural network toolbox (version 2.0) of the MATLAB1 package (version 4.0, MathWorks Inc., Natick, MA) and default
optimization parameters. By plotting the error evolution vs.
the number of iterations, we found that 50 iterations weregenerally enough to find the minimum in O(t,t'). In a limited
number of calibrations, a local minimum was found [O(t,t')
was larger than twice the average minimum value]. In thosecases the neural network calibration procedure was restarted
with a different set of random initial weights. Depending on
the neural network model, about six to 24 h on a P5-166 MHz
PC were required for calibration. Testing the calibrated neuralnetwork models with independent data needed only a few
seconds.
Hierarchical Approach
Two series of 12 water retention and 19 K, neural networkmodels were developed in a hierarchical approach using input
data listed in the first column of Table 4. The first 12 Ks models
use exactly the same input data as the retention models; thelast seven Ks models use input data that are not considered
for the retention models.
First, we carried out analyses using minimal input data, i.e.,
only textural class information (TXT in Table 4) similar tothe class RTFs of Carsel and Parrish (1988) and Wosten et al.
(1995). Next, neural networks were developed that used sand,silt, and clay fractions as input (SSC). Although one of thesefractions always contains redundant information, we used all1
Trade names are provided for the benefit of the reader and do
not imply endorsement by the USDA.

three as input to allow the neural network to select the most
relevant variable. In a subsequent model, bulk density (pb)
was added to the input variables (SSCpb). Next, we substituted
porosity (FOR) for pb (SSCPOR). Since bulk density andporosity are closely related, the change in input variable was
not expected to greatly affect the results. Information wasadded about the gravel content (particle size >2 mm) to yield
the SSCpt+Gravel model or the horizon information to give
SSCpb + HOR model. Horizon information was encoded as
one variable (HOR) with A horizons as 0, B horizons as 0.5,
and C horizons as 1.
Subsequent neural network models also considered one orseveral input variables derived from retention data. Although
it was the purpose of this study to predict retention parametersand Ks, limited retention measurements are sometimes available (cf. Soil Survey Staff, 1995). Assimilation of available
retention data in PTFs may improve the prediction of hydraulic properties (e.g., Rawls et al., 1992). Seven models used thesame input variables as the SSCpb model, but with additionof one or two retention points (6i0, 633, and 61500 at h = 10, 33,
and 1500 kPa, respectively; see Table 4). The water contents
were obtained by evaluating the fitted retention parametersof Eq. [1] at the appropriate pressure heads.
Four additional neural network models were developed for
Ks using the same input as the SSCpb model and effective
porosity (i.e., SSCpb + EP10, SSCpb + EP33, and SSCpb + EPi500
models) or estimated retention parameters (SSCpb + 0r0san).
Effective porosity was calculated by subtracting the watercontent at the appropriate pressure head (10, 33, or 1500 kPa)
from the total porosity. Following Ahuja et al. (1989) andMessing (1989), we also tried to predict K, from EP10 or EP33,
exclusively. Finally, we attempted to predict Ks using only
retention parameters (0r, Os, a, n ) so as to allow a comparison
with the PTF of Mishra and Parker (1989, 1992).
We evaluated TXT, SSC, SSCpb, and the SSCpb+033 models
by textural class because they appeared to be the most usefulmodels for practical applications.
Evaluation Criteria
The neural network predictions were evaluated in terms ofthe coefficient of determination (R2) between predicted andfitted or measured hydraulic parameters (6r, 6S, a, n, and Ks).
We did not evaluate the published water retention PTFs inthis way since most predict parameters for functions otherthan Eq. [1]. Additionally, we compared the predictions of
the neural networks and the published PTFs with measured
data as root mean square residuals (RMSRwrc and RMSRKs).
The RMSRwrc was obtained by converting the predicted parameters to water contents at the appropriate pressure heads
and calculating:- fl'.V

RMSRwrc = J-J-2[4]~ «parwhere TV is the number of samples (1209), L the number ofmeasured retention points for each sample, 0,, and 0,'; aremeasured and predicted water retention points, and «par is thenumber of parameters (four). The root mean square residualfor log(^s) was similarly calculated as:[5]where M is the number of samples (620), while log(/Q andlog(.Ks') denote measured and predicted values, respectively.

851

SCHAAP ET AL.: NEURAL NETWORK ANALYSIS FOR HYDRAULIC PROPERTY PREDICTION

For all models in the hierarchical approach listed in Table
4, we present average R2 and RMSR values of the 60 bootstrap
models. The variability of the bootstrap predictions will bequantified by standard deviations (the predicted parameters
had approximately normal distributions). The RMSR resultsof published PTFs were calculated for the entire data set of1209 water retention samples or the 620-sample Ks subset
without bootstrapping. Applying the bootstrap to evaluate thepublished PTFs is not useful because, contrary to the neuralnetworks, these models stay the same for all bootstrap data
sets.

RESULTS AND DISCUSSION
The results of the direct fit of Eq. [1] to retention data
for each textural class (Table 3) are generally similar to
those reported by Carsel and Parrish (1988) with thefollowing exceptions. The average 0S values were lower
for loam (0.356 m3 m~ 3 for our data vs. 0.43 m3 m~ 3 forCarsel and Parrish [1988]) and higher for clays (0.441
vs. 0.38 m3 m^3). Further, the values for log(a) weregenerally lower because we calculated the average oflog(a), but Carsel and Parrish (1988) computed the average of a. For example, for the loam class we foundlog(a) to be -2.11 log(cm~]) while Carsel and Parrish(1988) found -1.29 log^m"1). We also found higher
values for Ks for the silty clay loam class [0.87 vs. 0.22log(cm d"1)] and the clays [1.10 vs. 0.68 log(cm day"1)].
The average RMSRwrc of the direct fit varies from 0.010
to 0.025 m3 m~ 3 for sandy clay loam and silts, respectively; for the entire data set it was 0.020 m3 m~3.
The results of the hierarchical approach (Table 4)show an increase in the R2 and a decrease in RMSRwrc
and RMSR/& if more input variables are used. Thereis a substantial variation among the different output
variables. The behavior of R2 is discussed first.
All neural network models show a low R2 for 0r; thisindicates that 0r could not be predicted well from theavailable input data. A reason for this poor performanceis probably the asymptotic nature of 0r, which is an
extrapolated water content at infinite suction while thehighest suction available in most cases was not morethan 1500 kPa. The mean and standard deviation of thefitted 0r in Table 3 suggest a large coefficient of variationin each soil textural class. Table 4 shows that when
Oisoo was included in the input, the correlation remained
poor: the maximum R2 was 0.41 (SSCpb+0330i5oo)- Schaap

and Bouten (1996) also found poor correlations between
predicted and fitted 0r values. Vereecken et al. (1989)
found much higher ,R2 values (0.70) for 0r for a rangeof Belgian soils while the predicted 0r values tended tobe higher than our class averages because of a somewhat
different retention function by setting m = 1 in Eq. [1].
A comparison of R2 for 0S of the TXT and SSC models
with those of the SSCpb and SSCPOR models showsthat information about the soil bulk density is important
for accurately predicting 0S. If only basic soil properties
were used to predict 0S, the maximum R2 was 0.63
(SSCpb + Gravel). If water contents were used as additional input, R2 increased to 0.68 for the SSCpb + 0io
neural network model, and to 0.77 for the SSCpb +0io033 model.
The results for log(a) and log(n) were rather similar.
A slight increase in R2 is noticeable for the SSC model,
compared with the TXT model. Correlations were between 0.35 and 0.40 if only basic soil properties wereused. The addition of horizon information or gravel tothe input did not improve the prediction, nor did thereplacement of pb with porosity (POR). Addition of 010
or 033 to the input improved the correlations to 0.58 forboth log(a) and log(rc) (SSCpb + 033 model) and to 0.65
for log(ot) and 0.70 for log(n) if both 0 i0 and 033 wereused (SSCpb + 010033 model). The 0i500 variable did notprovide useful information for predicting log(a) orlog(n).
The correlation for log(/Q improved if the SSC model
was used instead of TXT; further improvement wasachieved with the SSCpb model (Table 4). Addition ofgravel or horizon information or replacing pb with POR
did not significantly improve the correlation. The maximum correlation that could be obtained with basic soil
properties alone was 0.58. Including information derived from the water retention data improved the correlation to as much as 0.72 (SSCpb + EPio). Information
contained in 0150o was apparently not as useful as adding
0io and 033 or EP10 and EP33. Notice that the results for
the SSCpb + 0io and SSCpb + 033 models were very
similar to those for the SSCpb + EPi0 and SSCpb + EP33
models. This finding is not surprising since the input
data were essentially the same (porosity can be derived
from bulk density while effective porosity can be calculated from porosity and water content). Using only EPi0
or EP33 as input, as was done by Ahuja et al. (1989)

Table 3. Average values and standard deviations (<r, in parentheses) of fitted hydraulic parameters, root mean square residuals (RMSR),
and the number of samples for nine textural classes.
Class
Sand
Loamy sand
Sandy loam
Loam
Silts*
Sandy clay loam
Clay loam
Silty clay loam
CIays§

log(a)

A'.rct

97
135
337
137
220
104
77
47
55

3
———————— m3 r-111—————————————————
0.044 (0.019)
0.413 (0.057)
0.039 (0.037)
0.395 (0.072)
0.389 (0.094)
0.031 (0.049)
0.054 (0.067)
0.356 (0.082)
0.065 (0.062)
0.441 (0.103)
0.076 (0.074)
0.379 (0.066)
0.091 (0.067)
0.439 (0.077)
0.111 (0.062)
0.460 (0.056)
0.081 (0.088)
0.441 (0.068)

t Total 1209 samples for water retention and 620 for A,.
t Silt and silty loam.
§ Clay, silty clay, and sandy day.

log(cm ')-1.57 (0.21)-1.49 (0.53)-1.57 (0.58)-2.11 (0.82)-2.51 (0.49)-1.80 (0.66)-1.95 (0.60)-2.36 (0.38)-1.89 (0.55)log(n)

RMSR,.

0.462 (0.200)
0.194 (0.108)
0.150 (0.094)
0.195 (0.140)
0.260 (0.131)
0.132 (0.100)
0.188 (0.128)
0.240 (0.110)
0.107 (0.059)

m3m 3
0.019
0.018
0.021
0.017
0.025
0.010
0.016
0.015
0.014

97
117
199
32

62
80
6
10
17

log (cm d ')2.71 (0.51)
1.92 (0.61)
1.53 (0.65)
0.99 (0.63)
1.04 (0.54)
1.29 (0.70)
0.67 (0.58)
0.87 (0.55)
1.10 (0.43)

852

SOIL SCI. SOC. AM. J., VOL. 62, JULY-AUGUST 1998

Table 4. Average values and standard deviations (a, in parentheses) of coefficients of determination (R2) and root mean square residuals(RMSR,,TC for water retention and RMSRKs for Ks) for 12 models that predict retention parameters and 19 models that predict Ks.
These results are evaluations of the models on independent data.
R2

Model inputt

o.

0,log(a)

log(n)

iog<*yRMSR,™
3

TXT

ssc

SSCpb
SSCPOR
SSCpb + Gravel
SSCpb + HOR

sscpb + 0,,,
SSCpb + 633
sscp,, + e™
sscpb + eMe33

sscpb + e,0e15oo
SSCpb + 6336™
SSCpb + EPio
SSCpb + EPjj
SSCpb + EP15M
SSCpb + 0,«,an
EPIO
EP,,
ere,an0.12 (0.03)
0.16 (0.05)
0.18 (0.04)
0.19 (0.04)
0.21 (0.05)
0.17 (0.06)
0.23 (0.05)
0.26 (0.06)
0.35 (0.07)
0.27 (0.05)
0.38 (0.06)
0.41 (0.07)-

0.11 (0.02)
0.10 (0.02)
0.54 (0.04)
0.53 (0.04)
0.63 (0.03)
0.53 (0.04)
0.68 (0.03)
0.56 (0.03)
0.52 (0.04)
0.77 (0.02)
0.68 (0.03)
0.55 (0.03)-

0.29 (0.03)
0.35 (0.04)
0.38 (0.04)
0.37 (0.04)
0.39 (0.04)
0.36 (0.03)
0.55 (0.03)
0.58 (0.03)
0.37 (0.04)
0.65 (0.03)
0.57 (0.04)
0.63 (0.03)-

0.32 (0.05)
0.37 (0.05)
0.38 (0.05)
0.39 (0.04)
0.40 (0.04)
0.38 (0.04)
0.55 (0.04)
0.58 (0.04)
0.39 (0.05)
0.70 (0.04)
0.57 (0.04)
0.64 (0.03)-

0.42 (0.04)
0.47 (0.04)
0.57 (0.04)
0.57 (0.04)
0.58 (0.04)
0.57 (0.04)
0.70 (0.03)
0.70 (0.03)
0.59 (0.04)
0.70 (0.03)
0.70 (0.03)
0.68 (0.04)
0.72 (0.03)
0.70 (0.03)
0.60 (0.03)
0.68 (0.03)
0.67 (0.03)
0.65 (0.04)
0.62 (0.03)

rn'rn
0.107 (0.004)
0.104 (0.006)
0.087 (0.005)
0.087 (0.004)
0.087 (0.006)
0.087 (0.004)
0.065 (0.003)
0.060 (0.005)
0.081 (0.004)
0.058 (0.006)
0.064 (0.005)
0.061 (0.004)-

RMSRKslog(cm d~*)0.627 (0.025)
0.602 (0.026)
0.533 (0.024)
0.536 (0.025)
0.539 (0.024)
0.537 (0.028)
0.448 (0.023)
0.451 (0.019)
0.529 (0.024)
0.448 (0.019)
0.447 (0.022)
0.463 (0.027)
0.435 (0.019)
0.452 (0.019)
0.523 (0.024)
0.467 (0.024)
0.479 (0.021)
0.494 (0.025)
0.504 (0.022)

t TXT = textural class; SSC = sand, silt, clay contents; pb = bnlk density; HOR = horizon; H,,,, 033, 0IHH, = water content at 10 kPa, 33 kPa, and 1500
kPa; EPio, EP33, EP1500 = effective porosity at 10 kPa, 33 kPa, and 1500 kPa; 0,, 0,, a, n = retention parameters in Eq. [2].

and Messing (1989), reduced the effectiveness of thepredictions somewhat compared with the SSCpb + EPi0
or SSCpb + EP33 Ks models. The same was true if only
van Genuchten (1980) parameters were used (cf. 6r0sa«
and SSCpb + 0r6san models), similar to the work ofMishra and Parker (1989,1992). These results show that
SSCpb and effective porosity data or fitted van Genuchten (1980) parameters do not convey exactly the same
information. Combining basic soil properties with retention data will lead to more accurate predictions.
The RMSR results for the retention models andlog(/Q models show similar but decreasing trends to thecorrelations for the individual parameters. The RMSRwrc
values were relatively high for the TXT input data; individual water retention data points were predicted witha standard deviation of 0.107 m3 m~ 3 while log(Ks) was
predicted with a RMSR& of 0.627 log(cm d~'). As observed by Williams et al. (1992) and Rawls et al. (1992),
inclusion of a retention point improved the prediction.
The use of two retention points did not improve theprediction much (RMSRwrc was 0.060 m3 m~ 3 for
SSCpb+e33 compared with 0.058 m3 m~ 3 for SSCpb +610653). Although the RMSRwrc values decreased whenmore input data were used, they remained high at 0.058
m3 m~ 3 (SSCpb + 610633), which is well above the RMSRwrcof the direct fit of Eq. [1] to the retention data (0.020
m3 m~3). The RMSR& values decreased to 0.447 log(cmd"1) for the SSCpb + 6i0615oo model. Notice that R2 and
RMSR values do not necessarily provide similar information: a comparison of the SSCpb + 633 and SSCpb +610633 models shows that the R2 values of most retention
parameters increased substantially while the RMSRwrc
decreased only modestly. The latter error measure is
probably of most direct interest.
The standard deviations in Table 4 provide information about the variability of the R2 and RMSR among
the predictions of the 60 bootstrap models. The standard
deviations were in most cases about 5 to 10% of theaverage R2 and RMSR values. This implies that models
with different performances are obtained when slightly
different calibration data sets are used. A single calibration, typical for most or all previously published PTFs,
would have had a chance of about 95% to be within
two standard deviations of the reported average R2 and
RMSR values. This result indicates that there is apparently a substantial uncertainty in the predicted hydraulic properties.
Evaluation of Selected Models by Textural ClassThe RMSRwrc values for the TXT, SSC, SSCpb andthe SSCpb + 633 models were generally the lowest forthe sandy loam and sandy clay loam classes (Table 5).
Notice that the RMSRwrc for the clay loam class showeda significant improvement from the SSCpb to the
SSCpb + 633 model, whereas the RMSRwrc for the sandclass decreased little even when 633 was included. For
this class much of the information about the particle-size
distribution is lost in the sand fraction. The prediction ofretention parameters for sandy soils probably could be
improved if more fractions within the sand particle-sizeclass were available (Vereecken et al., 1989; Schaap andBouten, 1996).
A direct comparison of the hierarchical models andpublished PTFs is not easy. First, the published PTFs
were calibrated on different data sets with possibly different distributions of data. Formally, we can only evaluate how well they perform on our, independent, data
set. But if we assume that our data set is a fairly good
representation of soils in the USA (cf., Fig. la and Ib),
some degree of comparison is still possible, especiallybecause we used the bootstrap to evaluate our neuralnetworks on independent data. A second reason thatcomplicates comparison is that neural networks andpublished PTFs use somewhat different input data. By
approximation, the input data for function R3 (Cosby

853

SCHAAP ET AL.: NEURAL NETWORK ANALYSIS FOR HYDRAULIC PROPERTY PREDICTION
Table 5. Root mean square residuals for water retention data (RMSR wrc ) obtained on .
four published pedotransfer functions (PTFs) for nine textural classes.

, samples with four hierarchical models andHierarchical models
Class
Sand
Loamy sand
Sandy loam
Loam
Siltsi
Sandy clay loam
Clay loam
Silty clay loam
Clays§
Data set average

TXT

SSC

0.096

————— nr m
0.092

SSCp,,

0.107

0.102
0.102
0.105
0.123

0.095

0.100
0.124
0.096
0.101
0.107
0.110
0.104

0.095

0.112
0.105
0.105
0.107

Published PTFsf

Rl

SSCp^s

R3

R4

0.107
0.126
0.115

0.113
0.112

0.097
0.105
0.195
0.145

0.104
0.134

0.087
0.091
0089
0.086

R2
————— m' m

36
49
124
52
82
38
28

0.080
0.054
0.050
0.057
0.071
0.052
0.050
0.065
0.073
0.060

0.092
0.100
0.072
0.081
0.096
0.078
0.096

0.106
0.094
0.087

17
20

0.104
0.114

0.091
0.086
0.114
0.084
0.109
0.099
0.109
0.101

0.097

0.159
0.126

135
218
106
78
47
56

0.111
0.071
0.130
0.112
0.136

0.090
0.120
0.126
0.112
0.111

0.132

99
136
335

0.098

t R1-R4: published water retention PTFs, see Table 1 for references. PTF R3 uses approximately the same input as the SSC neural network model while
PTFs Rl, R2, and R4 are comparable with the SSCpb model. Underlined values denote cases where the published PTFs provide better predictions thancorresponding neural networks.
I Silt and silt) loam.
§ Clay, silty clay, and sandy clay.

et al., 1984) model is comparable to the SSC model
while the input data for Rl and R2 (Rawls and Brakensiek, 1985) models and R4 (Vereecken et al., 1989)
are similar to the SSCpb model. The TXT and SSCpb +633 models do not have published PTFs with corresponding input data.
Based on average RMSRwrc values (bottom of Table5), the hierarchical neural network models provided
better estimates than the published PTFs by about 0.01
to 0.02 m3 m~ 3 for the SSC and SSCpb models. The resultsin Table 5 also show that, because of an additionaltransformation step in its derivation (cf., Table 3), function R2 of Rawls and Brakensiek (1985) did not predictvan Genuchten (1980) parameters as well as function
Rl predicted the Brooks and Corey (1964) parameters.
Coincidentally, function R2 was used by Carsel andParrish (1988) to calculate class-average van Genuchten(1980) parameters. A few (underlined) cases are shownin Table 5 where the previously published PTFs provided better estimates than corresponding hierarchicalneural network models. For example, function Rl predicted water retention slightly better for the silty clay
loam class than the comparable SSCpb neural networkmodel. Function R4 was somewhat better for the sand,loamy sand, and sandy clay loam classes. Following

Tietje and Tapkenhinrichs (1993) and Kern (1995), wenote that the function R4 (Vereecken et al., 1989) performed relatively well in comparison with other published PTFs. However, the neural network models yield
better results for our data set, both on average and formost textural classes.
Table 6 shows similar results for log(/Q predictedwith four hierarchical neural network models and six
published PTFs (K1-K6; Table 1). With respect to input
data use, function Kl (Cosby et al., 1984) is comparable
to the SSC model, while functions K2 (Brakensiek etal., 1984), K3 (Saxton et al., 1986), and K4 (Vereecken
et al., 1990) are similar to the SSCpb model. The PTFs
K5 (Ahuja et al., 1989) and K6 (Mishra and Parker,
1992) are comparable to the models using EP33 and the
9r6san, respectively (Table 4).
The results in Table 6 indicate that, for most texturalclasses, the RMSRXs decreased by about 0.2 log(cm d"1)
from the TXT models to the SSCpb + 633 model. The
RMSR of the sandy clay loam, clays, and clay loam
classes showed less pronounced decreases. All previously published PTFs had higher RMSRs and lower
R2 than the neural network models. In terms of average
RMSRfo, the best published PTF for predicting Ks was
Kl (Cosby et al., 1984). In terms of R2, the best published

Table 6. Root mean square residuals for data including saturated hydraulic conductivity (RMSRKs) obtained for Ntcst samples for four
hierarchical models and six published saturated hydraulic conductivity pedotransfer functions (PTFs).
Hierarchical models
Class

TXT

SSC
•

SSCp.
,

•«5 I*-1

Sand
Loamy sand
Sandy loam
Loam
SiltslSandy clay loam
Clay loam
Silty clay loam
Clays§
Dataset average
R1

0.510
0.622
0.652
0.643
0.546
0.724
0.684
0.637
0.491
0.627
0.42

0.433
0.620
0.638
0.667
0.487
0.711
0.642
0.468
0.453
0.602
0.47

Published PTFsf

J

Kl

K

0.408
0.604
0.521
0.488
0.369
0.668
0.573
0.423
0.506
0.533
0.57

0.319
0.416
0.466
0.470
0.372
0.598
0.586
0.363
0.459
0.451
0.70

36
45
74
12
23
30
2

3
6

0.613
0.666
0.795
0.860
0.550
0.957
0.980
0.624
0.545
0.746
0.30

K2

K5

K6

——— log (c- J ~ n ——————————————
0.421
0.519
0.883
0.362
0.670
1.090
0.877
0.511
0.625
0.815
0.627
0.777
0.587
0.847
0.530
0.611
0.862
1.153
0.569
0.462
0.868
1.344
1.178
1.155
0.624
0.790
0.630
0.678
1.003
0.526
0.960
1.852
1.441
0.914
1.704
0.674
0.791
0.761
0.934
0.822
0.43
0.42
0.22
0.54

K3

0.525
1.208
1.615
1.317
0.891
1.660
0.726
0.749
1.707
1.332
0.21

K4

97
117

199
32
62
80
6
10
17

f K1-K6: published saturated hydraulic conductivity PTFs, see Table 2 for references. PTF Kl is approximately comparable to the SSC neural networkmodel. PTFs K2, K3, and K4 are comparable to the SSCpb model. K5 and K6 are comparable to the EP33 and 6resan models in Table 4, respectively.
$ Silt and silty loam.
§ Clay, silty clay, and sandy clay.

854

SOIL SCI. SOC. AM. J., VOL. 62, JULY-AUGUST 1998

PTF was K5 (Ahuja et al., 1989). The overall R2 for thisPTF was 0.54, which is lower than the R2 values of 0.65
and 0.67 obtained for the EP10 and EP33 neural networkmodels (Table 4). The low R2 in Table 6 for K4 wasconsistent with the low calibration R2 listed in Table 2.
The results for the PTFs of Cosby et al. (1984) andSaxton et al. (1986) were not as good as those reportedby the researchers (0.839 and 0.95, respectively; cf., Table 2). Finally, we note that K6 performed relatively
poorly, with an RMSR/fs that exceeded one order ofmagnitude.

CONCLUSIONSIn this study we used a hierarchical neural network
approach to predict hydraulic parameters. The basis of
the analysis was a data set of 1209 samples, of which,horizon, sand, silt, and clay content, bulk density, porosity, and six to 13 retention points were known. A totalof 620 samples also had measured Ks data. We tested
12 different configurations of input data to predict water
retention parameters of van Genuchten (1980), and we
tried 19 different input configurations to predict Ks. The
results show that 6r was difficult to predict: R2 neverexceeded 0.41, even with the most relevant input. The
correlation for the other parameters [6S, log(a), log(«),
and log(Ks)] increased from relatively low levels for an
input configuration that used only textural classes, toabout 0.70 for configurations that used sand, silt, andclay contents, bulk density, and one or two water content
points. Root mean square residual water contents decreased from 0.107 m3 m~ 3 for a model that used only
information about the soil textural class, to about 0.058
m3 m~ 3 for the best neural network model. However,even with the addition of a measured water content to
the input, the RMSRwrc always stayed well above the
RMSRwrc obtained for direct fits (=0.020 m3 m~ 3 ). Measured water retention curves appear to contain information that cannot be predicted from macroscopic variables like sand, silt, and clay contents, bulk density, andwater retention. The RMSRfc results decreased fromabout 0.627 log(cm d~') for a model that used only
information about the soil textural class, to about 0.450log(cm d"1) for the best neural network models.
Uncertainty estimates for the predicted hydraulicproperties were generated by combining the hierarchical models with the bootstrap method. Results of thisapproach showed that there can be a considerable uncertainty in predicted hydraulic properties. These uncertainties can be used in conjunction with the Richards
equation to generate uncertainty estimates in simulated
unsaturated flow processes.
Four selected neural network models were found tobe better than four previously published PTFs for water
retention parameters. Only the PTF of Vereecken et al.
(1989) performed nearly as good as the neural networkmodel with similar input requirements. None of theinvestigated existing PTFs for Ks were found to be better
than our neural network models.
We think that using the models from our hierarchicalapproach is not only attractive because of improved

accuracy, but also because a neural network model canbe chosen to match a particular availability of data.
Because all neural network models were calibrated on
the same data set, predictions of different models will beconsistent. A final note may be that the neural networkmodels are probably somewhat more difficult to implement than more traditional regression models. However, because the neural networks used in this study arebasically simple matrix-vector operations, these limitations can be overcome. Currently, work is underway toimplement some of the neural network models in userfriendly software.
ACKNOWLEDGMENTS

This research was partially supported by NASA, Project(95) 805/60110. We would like to thank WJ. Rawls at the
USDA-ARS Hydrology Lab., Beltsville, MD, for providing
his soil moisture database.

LEU & TORIDE: SOLUTE TRANSPORT IN FINITE SOIL COLUMNS

855