Geoderma 86 Ž1998. 61–81

Using an inverse method to estimate the hydraulicproperties of crusted soils from tension-disc
infiltrometer data
a,)
ˇ ˚
Jirı
, Rafael Angulo-Jaramillo b, Marcel G. Schaap a ,ˇ´ Simunek
Jean-Pierre Vandervaere b, Martinus Th. van Genuchten aaU.S. Salinity Laboratory, USDA, ARS, 450 West Big Springs Road, RiÕerside,
CA 92507-4617, USAb
Laboratoire d’etudedes Transferts en Hydrologie et EnÕironment (LTHE, CNRS URA 1512,´
INPG, UJF) BP 53, 38041 Grenoble Cedex 9, France
Received 15 July 1997; accepted 24 March 1998

Abstract
An inverse procedure was used to estimate the soil hydraulic characteristics of a two-layeredsoil system—soil surface crust and subsoil—from data obtained during a tension-disc infiltration
experiment. The inverse procedure combined the Levenberg–Marquardt nonlinear parameter
optimization method with a numerical solution of the axisymmetric variably-saturated flow
equation. The objective function was defined in terms of the cumulative infiltration curve and thefinal water content measured directly below the tension-disc infiltrometer at the end of theexperiment; this final water content was assumed to correspond to the final supply pressure head.
We analyzed two infiltration experiments carried out with a 25-cm diameter tension-disc infiltrometer. One experiment was carried out on a two-layered system, and a second after removal of thesurface crust covering the sandy subsoil. Both experiments were performed with six consecutive
supply tensions. We first analyzed the infiltration experiment for the subsoil only, thus yielding its
hydraulic characteristics. Subsequent analysis of the infiltration experiment for the two-layered
system with known hydraulic properties of the subsoil provided estimates of the hydraulicproperties of the surface crust. We further compared the estimated hydraulic parameters of the
subsoil with those obtained using Wooding’s analytical method wWooding, R.A., 1968. Steady
infiltration from a shallow circular pond. Water Resour. Res. 4, 1259–1273x and predictions basedon a neural network model requiring textural input information. All three methods generated
roughly the same results. The numerical inversion technique proved to be a convenient tool for)

Corresponding author. Tel: q1-909-369-4865; Fax: q1-909-342-4964; E-mail:jsimunek@ussl.ars.usda.gov
0016-7061r98r$ - see front matter q 1998 Elsevier Science B.V. All rights reserved.
PII: S 0 0 1 6 - 7 0 6 1 Ž 9 8 . 0 0 0 3 5 - 462ˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

estimating the soil hydraulic properties of both the surface crust and the subsoil. q 1998 Elsevier
Science B.V. All rights reserved.
Keywords: parameter estimation; soil hydraulic properties; tension-disc infiltrometer; neuralnetworks

1. Introduction
A large number of measurements were carried out during the soil program
part of the HAPEX ŽHydrologic Atmospheric Pilot Experiment.-Sahel Žclimatic
region bordering the southern Sahara from Senegal to Somalia. regional-scale
experiment Ž Goutorbe et al., 1994; Cuenca et al., 1997. . Soil-related measurements included neutron probe and TDR measurements of the profile soil watercontent, tensiometer measurements of the soil water potential, particle-size
analyses, and infiltration rate measurements. The type of data collected and theexperimental set-ups are documented by Cuenca et al. Ž 1997. and Vandervaere
et al. Ž1997. , among others. The overall objective of the experiment was toimprove the understanding of the physical processes of the Earth’s ecosystem,
the hydrologic system, and the atmospheric system ŽCuenca et al., 1997..
An important feature affecting the water balance in many tropical andsubtropical regions is the development of surface crusts in loamy and sandy soilsŽCasenave and Valentin, 1992; Valentin and Bresson, 1992. . These crusts
control the infiltration process and play a critical role in separating runoff andinfiltrating water from the incoming rainfall. Valentin and Bresson Ž 1992.
classified surface crusts according to their formation, leading to three majorclasses: structural crusts, erosion crusts and depositional crusts, each withseveral subclasses. While structural crusts are formed by the impact of water
drops, depositional crusts are formed by the translocation of fine particles andtheir deposition at a certain distance from their original location Ž Valentin andBresson, 1992.. Valentin and Bresson characterized erosion crusts as being
rigid, thin and smooth surface layers enriched in fine particles, with the coarser
particles removed by wind or overland flow. The surface crusts have usually ahigher bulk density and smaller pores, and consequently a lower saturated
conductivity than the underlying soil. Although the thickness of a surface crust
is generally only a few millimeters, the small saturated hydraulic conductivity
can significantly decrease the infiltration rate.
Several quasi-analytical solutions ŽHillel and Gardner, 1969, 1970; Ahuja,
1974; Smilles et al., 1982. have been developed and used to calculate fluxes
during steady and transient infiltration into a crusted soil profile from knowledgeof the hydraulic properties of the crust and the underlying soil. These solutions
are all for one-dimensional water flow. Because surface crusts are typically very
thin, few field methods exist to estimate their hydraulic properties ŽVandervaere
et al., 1997.. Vandervaere et al. Ž1997. recently proposed a method which usesˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

63

tension-disc ŽPerroux and White, 1988; Ankeny et al., 1991; Reynolds andElrick, 1991; Logsdon et al., 1993, among many others. infiltration data atseveral water supply potentials, together with information from pre-installedmini-tensiometers below the soil crust to estimate the hydraulic conductivity,
matric flux potential and sorptivity. Recent developments in transient flow
theory ŽHaverkamp et al., 1994. were used for the analysis.
Another promising method for describing the soil crustrsubsoil multi-layered
system, and as an alternative to the traditional analysis of tension-disc infiltration experiments based on the analytical solution of Wooding Ž1968. , assuming
a homogeneous profile, is the use of parameter estimation using numericalˇ ˚
inversion. Simunek
and van Genuchten Ž1996, 1997. recently suggested a
numerical inversion method to estimate the hydraulic properties from cumulative
infiltration data obtained with a tension-disc permeameter at several consecutive
tensions. Information required for a successful inversion, in addition to the
cumulative infiltration with time, included information about the initial and finalˇ ˚
water contents before and after the infiltration experiment Ž Simunek
and van
Genuchten, 1997. . Using numerically simulated data, they showed that themethod can provide information on not only the hydraulic conductivity as is
usually the case when invoking quasi-analytical methods, but also the soil water
retention curve.
Soil hydraulic properties can also be estimated indirectly from correlations
with basic soil data such as the sand, silt and clay fractions, bulk density andrororganic matter content. Such indirect methods, also called pedotransfer functionsŽPTFs., are often used to generate soil hydraulic properties in situations where
measurements are too expensive, too cumbersome, or too difficult to carry out.
Several PTFs predicting water retention and the saturated hydraulic conductivity, K s , were tested on independent data by Tietje and Tapkenhinrichs Ž1993.,
Kern Ž1995. , Williams et al. Ž 1992. , and Tietje and Hennings Ž 1996. . Thesestudies show that many different types of PTFs exist in terms of input data, thepredicted properties, mathematical structure and accuracy. Recently Pachepsky
et al. Ž1996., Schaap and Bouten Ž 1996. and Tamari et al. Ž 1996. used neuralnetwork analyses to improve the PTF predictions. Schaap et al. Ž 1998. showedthat a hierarchical neural network system with increasing predictive capability
could be obtained when more input variables were used for the estimation of theretention parameters of van Genuchten Ž 1980. and K s . They also showed thatneural networks provided significantly better predictions than 10 previously
published PTFs.
In this paper we will analyze two tension-disc infiltration experiments inorder to estimate the soil hydraulic properties of a two-layered system involving
a thin surface crust and the subsoil. We will first analyze the infiltration
experiment carried out on the subsoil to obtain its hydraulic properties. This
information is subsequently used to estimate the hydraulic properties of thesurface crust by analyzing infiltration into the two-layered system. First, we will

64

ˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunekˇ ˚
use the method of Simunek
and van Genuchten Ž1997. based on a numericalanalysis of the field experiment. Results for the subsoil will be compared withthose obtained with the classical steady-state analysis of Wooding Ž 1968. .
Results will also be compared with soil hydraulic functions predicted fromtextural information using neural network analysis, as well as those given by
Carsel and Parrish Ž 1988. for a particular textural class. The neural network
analysis will be applied only to the subsoil since no soil textural information wasavailable for the very thin soil crust.

2. Experimental data
A large number of infiltration measurements with tension infiltrometersŽReynolds and Elrick, 1991; Ankeny et al., 1991; Thony et al., 1991. were taken
as part of soil program of the HAPEX-Sahel regional-scale experiment during
the summers of 1992 and 1993. The measurements were carried out at theCentral Super Site East ŽCuenca et al., 1997. using tension infiltrometers having
disc diameters of 25, 8 and 4.85 cm. In this study we will analyze two
experiments conducted on a fallow site with an erosion type crust using a 25-cm
diameter disc. A thin Žf 0.3 cm. layer of sand was placed between the discmembrane and the soil surface to improve hydraulic contact. Readings of thewater supply tube were done visually. The soil was characterized according to
the textural triangle of U.S. Department of Agriculture as a loamy sand with aparticle size distribution given in Table 1. In one experiment Ž denoted as T2. thedisc permeameter was placed directly on top of the soil surface Ž on the surface
crust., while for a second experiment Ž denoted as T1. the crust was removed andthe infiltrometer placed on the subsoil. The left-hand parts of Tables 2 and 3summarize for the two experiments the applied supply pressure heads, h 0 , thetime intervals during which h 0 was applied Ž t init is the starting time, and t final
Table 1
Particle size distribution of the loamy sandParticle diameter range Žmm.

Percent by mass

0.000–0.002
0.002–0.020
0.020–0.050
0.050–0.100
0.100–0.200
0.200–0.500
0.500–1.000
1.000–2.000

11.2
1.9
2.5
12.0
26.6
37.7
7.3
0.8

ˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

65

Table 2
Wooding’s analysis of experiment T1 Žno crust.
h0
wcmx

t initwsx

t finalwsx

QŽ h 0 .
wml sy1 xhawcmx

K Ž ha .
wcm sy1 xa)wcmy1 xKs
wcm sy1 x

y11.5
y9.0
y6.0
y3.0
y1.0
y0.1

0
2430
3510
4470
5295
6015

2250
333
4290
5115
5835
6450

1.275
1.564
1.833
2.224
2.608
3.233

y10.25
y7.5
y4.5
y2.0
y0.55

0.001280.001180.00159
0.002160.004160.08180.05300.06440.07980.23840.002960.00176
0.002130.002530.00473

the final time. and the final Ž ‘steady-state’. infiltration rates, QŽ h 0 ., for aparticular supply pressure head. Tables 2 and 3 also give results of Wooding’s
analysis Žto be discussed later. .
Disturbed gravimetric samples were taken at two depths before and after the
infiltration experiments to determine the initial and final water contents. Undisturbed samples were also collected to obtain the dry bulk density needed toconvert gravimetric water contents into volumetric values. The bulk density was
found to be highest for the top soil Ž 1.794 g cmy3 for the soil layer to 7.5 cm.
and gradually decreased with depth Ž1.632 g cmy3 at a depth of 17.5 cm. . The
initial water contents in the surface layer for experiments T1 and T2 were 0.149
and 0.0369 cm3 cmy3, respectively, and the final water contents below thetension-disc infiltrometer at the end of the infiltration experiments 0.275 and
0.232, respectively.

3. Data analysis
The tension infiltrometer data will be analyzed in two ways. One method usesˇ ˚
a numerical solution of the variably-saturated water flow equation ŽSimunek
etTable 3
Wooding’s analysis of experiment T2 Žwith crust.
h0
wcmx

t initwsx

t finalwsx

QŽ h 0 .
wml sy1 xhawcmx

K Ž ha .
wcm sy1 xa)wcm y1 x

Ks
wcm sy1 x

y10.0
y7.0
y5.0
y3.0
y1.0
y0.5

0
2040
3240
4390
5520
7320

1740
2940
4080
5220
7020
8115

0.06390.06430.06520.06940.0676
0.0677

y8.5
y6.0
y4.0
y2.0
y0.75

0.223=10y5
0.912=10y5
0.316=10y4
y0.204=10y4
0.236=10y5

0.00177
0.00757
0.03054y0.0130
0.01770.226=10y5
0.955=10y5
0.357=10y4
y0.199=10y4
0.236=10y5

ˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

66al., 1996. coupled with the minimization method of Marquardt–LevenbergŽMarquardt, 1963.. The objective function being minimized will be defined interms of cumulative infiltration data and the final water content at the end of theinfiltration experiment, while the initial condition is given in terms of the waterˇ ˚
content Ž Simunek
and van Genuchten, 1997. . The final water content below thetension-disc infiltrometer, u f , is assumed to correspond to the final supply
pressure head at the end of the experiment. The data will be analysed also using
the analytical solution of Wooding Ž1968. to the final Ž ‘steady-state’. infiltration
fluxes at particular supply pressure heads. The soil hydraulic parameters of the
subsoil will also be estimated by neural network prediction from textural
information ŽSchaap et al., 1998. .
3.1. Numerical model
The governing flow equation for radially symmetric isothermal Darcian flow
in a variably-saturated isotropic rigid porous medium is given by the following
modified form of the Richards’ equation ŽWarrick, 1992. :

Eu
Et

1 E
s

r Er

ž

rK

Eh
Er

E

/

q

Eh

ž /
K

Ez

Ez

EK
q

Ez

Ž1.

where u is the volumetric water content wL3 Ly3 x, h is the pressure head wLx, K
is the hydraulic conductivity wL Ty1 x, r is a radial coordinate wLx, z is a verticalcoordinate wLx positive upward, and t is time wTx. Note that we assume in Eq. Ž1.
that the porous medium is isotropic, although not necessary homogeneous.
Initial and boundary equations applicable to a disc-tension infiltrometer experiment ŽWarrick, 1992. are as follows:u Ž r , z ,t . s u i Ž z .

ts0

Ž2.

h Ž r , z ,t . s h 0 Ž t .

0 - r - r0 , z s 0

Ž3.

r ) r0 , z s 0

Ž4.

y

E h Ž r , z ,t .
Ez

y1s0

u Ž r , z ,t . s u i Ž z .

r2qz2™`Ž5.

where u i is the initial water content wL3 Ly3 x, h 0 Ž t . is the time-variable supply
pressure head imposed by the tension-disc infiltrometer wLx and r 0 is the discradius wLx. Note that the initial condition is given in terms of the water content.
ˇ ˚
Simunekand van Genuchten Ž 1997. showed that having this initial condition, as
compared to using the pressure head, insures a more stable and unique solutionof the inverse problem. Soil surface boundary conditions below the discpermeameter and the remaining soil surface are represented by Eqs. Ž 3. and Ž 4. ,respectively. Eq. Ž5. assumes that all subsurface boundaries are so far removedˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

67

from the supply source that they are not affected during the infiltration experiment. Eq. Ž 1., subject to the above initial and boundary conditions, was solved
numerically using a quasi-three-dimensional Žaxisymmetric. finite element codeˇ ˚
HYDRUS-2D as documented by Simunek
et al. Ž1996..
3.2. Soil hydraulic properties
A model of the unsaturated soil hydraulic properties must be selected prior toapplication of the numerical solution of the Richards’ equation. Most analyticalsolutions for infiltration are based on the exponential function of Gardner
Ž1958., which describes the dependence of the unsaturated hydraulic conductivity on the pressure head. This function, however, was found to be the least
successful ŽLeij et al., 1997. of five different functions studied in terms
describing observed hydraulic conductivity data taken from the UNSODA
unsaturated soil hydraulic database Ž Leij et al., 1996. . In our analysis we will
use the unsaturated soil hydraulic functions of van Genuchten Ž1980. which
performed better. These functions are given by
Se Ž h . s

u Ž h . y ur
us y ur1s

Ž1 q < a h < n .

u Ž h . s us
K Ž u . s K s Sel 1 y Ž 1 y Se1r m .
K Žu . s Ks

m

h-0

Ž6.

hG0
m 2

h-0

Ž7.

hG0

where Se is the effective water content wyx, K s is the saturated hydraulicconductivity wL Ty1 x, ur and us denote the residual and saturated water contents
wL3 Ly3 x, respectively; l is a pore-connectivity parameter wyx, and a wLy1 x andn wyx and m Žs 1 y 1rn. wyx are empirical parameters. The pore-connectivity
parameter l in the hydraulic conductivity function was estimated by Mualem
Ž1976. to be 0.5 as an average for many soils. Assuming l s 0.5, the hydrauliccharacteristics defined by Eqs. Ž 6. and Ž 7. contain five unknown parameters: ur ,us , a , n, and K s . Since saturation will never be reached during a tension
infiltration experiment, K s and us in this study are interpreted as being
extrapolated, empirical parameters outside the range of the disc experimentˇ ˚ŽSimunekand van Genuchten, 1996. . Also, tension-disc infiltration in general isa wetting process Ž assuming that one can neglect internal drainage at the initialpressure head and during short periods when the infiltration process is interrupted in order to refill the water supply tube. ; this means that the hydraulicparameters in Eqs. Ž6. and Ž7. represent wetting branches of the unsaturated
hydraulic properties.

ˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

68

3.3. InÕerse solution
The objective function F to be minimized during the parameter estimation
process can be formulated using either cumulative infiltration data only, or
cumulative infiltration data in combination with additional information such as
the observed water content corresponding to the final supply pressure head, ortransient pressure head andror water content measurements within the soil
profile. The objective function for multiple measurement sets is defined asˇ ˚ŽSimunekand van Genuchten, 1996.
nj

m

ž

F Ž b ,q1 , . . . ,qm . s Ý Õj Ý wi j q j ) Ž t i . y q j Ž t i , b .
js1

js1

2

/

Ž8.

where m represents the different sets of measurements such as the cumulative
infiltration data Ž I . and the final water content Ž u f . , n j is the number ofmeasurements in a particular set, q j )Ž t i . are specific measurements at time t ifor the jth measurement set, b is the vector of optimized parameters Ž e.g., ur ,us , a , n, and K s ., q j Ž t i , b . are the corresponding model predictions for theparameter vector b , and Õj and wi , j are weights associated with a particularmeasurement set or point, respectively. We assume for now that the weighting
coefficients wi , j in Eq. Ž8. are equal to one, that is, the variances of the errors
inside a particular measurement set are assumed the same. The weighting
coefficients Õj , which minimize differences in weighting between different datatypes because of different absolute values and numbers of data involved, aregiven by ŽClausnitzer and Hopmans, 1995.
1
Õj s

n j sj 2Ž9.

thus defining the objective function as the average weighted squared deviation
normalized by measurement variances sj 2. Since the final water content is only
one number and its variance can not be defined, the weight for this data point isassumed to be one.
Minimization of the objective function F is accomplished by using theLevenberg–Marquardt nonlinear minimization method Ž Marquardt, 1963. . Thismethod was found to be very effective and has become a standard in nonlinear
least-square fitting among soil scientists and hydrologists Ž van Genuchten, 1981;
Kool et al., 1985, 1987. .
3.4. Wooding’s analytical solution
The variably-saturated water flow equation can be solved analytically for
infiltration from a circular source with a constant pressure head at the soilˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

69

surface and with the unsaturated hydraulic conductivity described by the exponential model of Gardner Ž1958.
K Ž h . s K s exp Ž a ) h .

Ž 10.

The analytical solution was derived by Wooding Ž1968. and is given byQ Ž h 0 . s p r 02 q

ž

4 r0

a)

/ Ž .
K h0

Ž 11 .

where Q is the steady-state infiltration rate wL3 Ty1 x, r 0 is the radius of the discwLx, h 0 is the wetting pressure head wLx, K Ž h 0 . is the unsaturated hydraulicconductivity wL Ty1 x at pressure head h 0 , and a ) is the sorptive numberŽWhite and Sully, 1987. . The first term on the right represents the effect ofgravitational forces and the second term the effect of capillary forces. Wooding’s
analytical solution has two unknown variables, i.e., the hydraulic conductivity
K Ž h 0 . Žor the saturated hydraulic conductivity K s after using Gardner’s equation. and the sorptive number a ). The traditional analysis of tension-discexperiments using the solution of Wooding Ž 1968. hence requires two steadystate fluxes obtained with the same disc infiltrometer at different tensionsŽAnkeny et al., 1991. , or with two infiltrometers of different diameters at thesame tension Ž Smettem and Clothier, 1989. . We note here that Wooding’s
analytical solution assumes that the soil is homogeneous, isotropic, and has auniform initial water content. When these restrictive assumptions are notfulfilled, Wooding’s analysis often leads to negative values of K s ŽHussen and
Warrick, 1993; Logsdon and Jaynes, 1993; Vandervaere et al., 1997. .
The theory for obtaining the unsaturated hydraulic conductivity in the middleof an interval between two successively applied pressure heads was described
previously by Ankeny et al. Ž 1991. , Reynolds and Elrick Ž 1991. and Jarvis andMessing Ž1995., among others. The approach assumes that the sorptive numbera ) in Gardner’s exponential model is constant over the interval between two
adjacent supply pressure heads such thatln)
a iq1r2
s

Qi
Q iq1

h i y h iq1

i s 1, . . . ,n y 1

Ž 12 .

where n is the number of infiltration tensions used, and where the subscript 1r2
on a ) indicates estimation in the middle of two adjacent supply pressure heads:
h iq1r2 s Ž h i q h iq1 .r2. The unsaturated hydraulic conductivity at pressure head
h iq1r2 is then calculated as
K iq1r2 s

Q iq1r2
4 r0
p r 02 q )a iq1r2

i s 1, . . . ,n y 1

Ž 13.

70

ˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

in which the estimated infiltration rate Q iq1r2 in the middle Ž h iq1r2 . between
two adjacent supply pressure heads is calculated as a geometric mean of theactual infiltration rates Q i and Q iq1:
ln Q i q ln Q iq1
Q iq1r2 s exp
i s 1, . . . ,n y 1
Ž 14.
2
The saturated hydraulic conductivity K s can be calculated from Eq. Ž 10. using)
known values of h iq1r2 , K iq1r2 and a iq1r2
as follows
K iq1r2
Ks s
Ž 15.
)
exp Ž a iq1r2
h iq1r2 .

4. Results and discussion
We first analyze experiment T1 for the case without a surface crust. The flow
domain for the numerical solution was discretized into 722 triangular elements
and 200 nodes, with a relatively fine grid close to the soil surface and below thedisc, and larger elements farther away from the disc. Limiting the number of
elements is important to avoid excessive CPU times since the direct numericalsolution has to be carried out many times during the parameter optimization
process. The adopted discretization was thought to be a good compromise
between solution precision and required computer time. The relatively coarse
discretization scheme resulted in mass balance errors which were always smaller
than 0.1%.
Fig. 1 shows the measured and calculated cumulative infiltration curves, aswell as their differences, for experiment T1 Ž no crust. when the objectivefunction, F Ž b , I, u f ., was defined in terms of the cumulative infiltration, I, andthe water content, u f , corresponding to the final pressure head. The latter
information was included as a point on the retention curve, i.e., u Ž h. s
0.275Žy0.1 cm.. The soil profile was considered homogenous with an initialwater content of 0.149. Notice the small breaks in the cumulative infiltration
curve ŽFig. 1a. caused by momentary removal of the infiltrometer from the soil
surface while resupplying it with water and adjusting the tension for a new timeinterval Ž t final , t init . Žsee Tables 2 and 3.. The interruptions in infiltration took
about 3 and 5 min for experiments T1 and T2, respectively. More than 10.5 l ofwater infiltrated during experiment T1. The supply tubes could not hold thismuch water and water hence had to be resupplied during the runs. Theinterruptions, however, did not cause any problem for the numerical solution
since during these times the Dirichlet boundary condition below the discpermeameter could be simply replaced by a zero-flux Neumann type boundary
condition. Notice the very close agreement between the measured and optimized
cumulative infiltration curves ŽFig. 1b. ; the largest deviations were generally
less than 60 ml, being only about 0.5% of the total infiltration volume.

ˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

71

Fig. 1. Experiment T1: measured and optimized cumulative infiltration curves Ža. and theirdifferences Žb. for the soil profile without the surface crust.

Since nearly all optimizations initially gave ur estimates equal to or very
close to zero, we decided to fix ur to zero for all cases. Van Genuchten’shydraulic parameters obtained by numerical inversion, together with the final
value of the objective function F f , are given in Table 4. Notice the relatively
high precision of the final parameter estimates as judged from the normalized
standard deviation of the parameter estimate, NSD w%x given bys
NSD s 100
Ž 16.
bj
where s is the standard deviation of the parameter estimate, calculated fromknowledge of the objective function, the number of observations, the number of
unknown fitted parameters, and an inverse matrix Ž van Genuchten et al., 1991. .
The narrow confidence interval for us reflects the large weight given to the final

72

Parameters

Numerical inversion

F Ž b , I, u f .
urusa wcmy1 xn
K s wcm sy1 x
Ff

Neural network
prediction

F Ž b , I, up .

Value

NSD w%x

Value

NSD w%x

Value

NSD w%x0.0
0.275
0.03031.91
0.002210.000187

y
0.03
2.70
3.49
4.08

0.0
0.296
0.03062.016
0.00209
0.000235

y
0.04
0.95
2.21
2.9

0.05280.32660.02191.65270.000416

23.8
4.35
26.9
7.42
25.1

Value from Carsel
andParrish Ž1988.
0.057
0.41
0.124
2.28
0.00405ˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

Table 4
Soil hydraulic parameters for loamy sand subsoilˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

73

water content in the optimization. Table 5 shows the correlation matrix for the
optimized parameters. The saturated water content, us , was uncorrelated withany other optimized parameters. However, a high negative correlation exists
between K s and n. This correlation reflects the shape of the objective function
which shows a hyperbolical valley in the parameter plane K s and n Žsee figure 6ˇ ˚
in Simunek
and van Genuchten, 1997. . Nevertheless, even these two parameters
were estimated with a relatively high level of confidence ŽTable 4..
The saturated water content, us , estimated by numerical inversion was equal
to 0.275 which is about 0.045 less than the porosity as calculated from the bulk
density. This result is consistent with observations Ž e.g., Klute, 1986; Baumer,
1992. that field-measured us-values are generally about 10–20% lower than theporosity. Also contributing to the observed difference in this case is likely thehigh permeability of the loamy sand subsoil Ž K s was nearly 200 cm dy1 . . The
high value of K s may have caused already significant drainage and redistribution of water in the soil profile during the approximately 60 s it took after
completion of the infiltration experiment to remove the tension disc and collectthe soil samples for laboratory determination of the water content. Hence, thecollected samples may not have corresponded exactly with the final supply
pressure head.
To further study this redistribution problem we repeated the numericalinversion with a modified objective function F Ž b , I, up . which included thefinal water content not as a fixed point of the retention curve, but as a transient
value of the water content, up , below the disc measured 60 s after the end of theinfiltration experiment. This approach permits the numerical model to calculate
the amount of water draining from the upper soil profile before actual sampling
took place. The results of this numerical inversion are also presented in Table 4.
The objective function has now a somewhat higher, but not significantly
different value. The estimated us was slightly higher Ž0.296. than in theprevious case, being now only about 0.025 lower than the porosity. Values ofthe other optimized parameters are about the same. The correlation matrix Ž notfurther presented here. also had a similar structure.
The measured cumulative infiltration curve was also analyzed using Wooding’s analytical solution as given by Eqs. Ž 10. through Ž 15. . Wooding’s analysis
requires steady-state fluxes; if steady-state is not reached the soil hydraulic
Table 5
Correlation matrix for optimized parameters for experiment T1
Parameter

a

n

us

Ks

a
n
us
Ks

1.0
0.899
0.032
y0.8131.0
0.013
y0.9861.0
y0.0091.0

ˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

74

conductivities will be somewhat under- or over-estimated. The minimum timerequired for the infiltration rate to approach steady-state depends on soil texture,
the size of the infiltration disc and the supply pressure head. The steady-state
infiltration rate is reached relatively fast for coarse-textured soils, as in our case.
The time, tgrav , after which gravitational forces dominate capillary effects can becalculated from knowledge of the sorptivity and the conductivity as followsŽPhilip, 1969.:
tgrav s

S

ž /

2

Ž 17.
K
where S is sorptivity wLTy0.5 x. From values given in Vandervaere et al. Ž1997.
and from our estimates of the hydraulic conductivity, we calculated that tgrav isapproximately equal to 4 and 22 min for supply tensions of 1 and 10 cm,
respectively. The time intervals used for our supply tensions were in most cases
larger than tgrav .
We obtained the final infiltration rates for particular supply pressure heads bydetermining the slope of the cumulative infiltration curve Ž Fig. 1a. using linear
regression analysis on the last six measured values. The results of Wooding’s
analysis are summarized in Table 2. This table gives unsaturated hydraulicconductivities K Ž h a . corresponding to the middle of the pressure head interval,
h a , between two consecutive supply pressure heads, the sorptive number a ) ,and the saturated hydraulic conductivity, K s , as extrapolated from K Ž h a . using
the estimated a ). All calculated K s-values were of the same order, except whenestimated from the highest supply pressure head interval in which case K s was
about twice as high. Also, a ) for this pressure head interval was about threetimes the other values.
Fig. 2 compares the unsaturated hydraulic conductivity function as obtainedwith the numerical parameter estimation method with local values of thehydraulic conductivity obtained for pressure heads in the middle of two supply
tensions using the Wooding’s analysis Ž Table 2. . Both methods gave almost
identical results for pressure heads in the interval between y2 and y10.25 cm.
The value of the hydraulic conductivity for the highest pressure head interval
was overestimated by a factor of two using Wooding’s analysis. Fig. 2 showsthat when the underlying assumptions of the two methods are satisfied Ž i.e.,
approximate homogeneity of the soil profile and steady-state flux for analyticalmethod., both methods provide approximately the same results.
Table 4 gives, in addition to van Genuchten’s parameters obtained bynumerical inversion, the hydraulic parameters obtained with a neural network
prediction from the measured soil texture classes and bulk density, and theaverage parameter values for the loamy sand soil textural group as estimated by
Carsel and Parrish Ž 1988. from an analysis of a large number of soils. We used
two neural network models derived by Schaap et al. Ž1998. to predict theretention parameters of Eq. Ž 6. , as well as K s , from sand, silt and clayˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

75

Fig. 2. Unsaturated hydraulic conductivities of the subsoil calculated using Wooding’s analyticalsolution for particular pressure heads and using the numerical inversion.

percentages Žderived from Table 1. and the subsoil bulk density. The three-layer
back-propagation neural networks were calibrated and validated on a large soil
database consisting of 1209 samples from the USA ŽRawls, personal communication.. The variation of textural, bulk density and hydraulic data in thisdatabase makes the calibrated neural network widely applicable. Validation
results showed that the neural networks performed better than traditional regression equations. Water retention and saturated hydraulic conductivities werepredicted with an average accuracy of, respectively, 0.087 cm3 cmy3 and0.53 logŽcm dy1 . over the entire data set ŽSchaap et al., 1998..
Fig. 3 shows both the retention and the hydraulic conductivity functions forthe parameters in Table 4. Notice the relatively good agreement between thehydraulic properties obtained by numerical inversion and the neural network
prediction. The entire retention curve obtained by numerical inversion using
objective function F Ž b , I, u f . was shifted by about 0.05 cm3 cmy3 downwards
as compared to the neural network prediction, whereas the water content interval
between us and ur was about the same Žequal to 0.27.. The initial condition forexperiment T2 was 0.036 which is lower than ur predicted with the neuralnetwork; this result explains the fact that our numerical inversion gave a smaller
value of ur . The neural network model predicted us to be more or less equal tothe porosity. The values of a predicted by numerical inversion and the neuralnetwork were about the same, as was the case with n. However, K s predictedwith the neural network was about five times lower than K s predicted bynumerical inversion. The value of n obtained by numerical inversion was
similar to the average n of the loamy sand soil textural group as given by Carsel
and Parrish Ž 1988. . The highest difference between parameters obtained withnumerical inversion and average parameters of the loamy sand textural groupwere for a . Parameters for textural groups as derived by Carsel and Parrish

76ˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

Fig. 3. Soil hydraulic functions obtained with numerical inversion for experiment T1, using neuralnetwork prediction for data from Table 1, using Wooding’s analysis, and for loamy sand texturalgroup ŽCarsel and Parrish, 1988..

Ž1988. represent an entire soil class, and hence some differences should be
expected. Neural network predictions of the soil hydraulic parameters, on theother hand, provide parameter values for a specific soil, and hence should yield
better agreement.
In the second step we analyzed experiment T2 when the tension permeameter
was placed directly on the soil surface. The cumulative infiltration data wereagain analyzed by numerical inversion and using Wooding’s method. We stresshere that the assumptions of Wooding’s solution were not fulfilled since the soil
profile was not homogeneous but consisted of two layers, a crust and thesubsoil. This complication, however, is not an issue for the numerical solution
which readily accommodates multi-layered profiles. The soil profile was assumed to consist of a 0.5 cm thick top layer representing the surface crustˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

77Žerosion crusts are in general only about 0.1 cm thick ŽCasenave and Valentin,
1992.; we used 0.5 cm for numerical reasons. , and the underlying subsoil where
hydraulic parameters were assumed known from the earlier analysis involving
the subsoil only Ž experiment T1 in Table 4. . Hence, only the parameters of thesurface layer were now assumed unknown and sought by numerical inversion.
The optimized parameters are given in Table 6. The results indicate much lower
confidence in the estimated parameter values for the surface crust as comparedto parameters for subsoil Ž Table 4. .
The measured and optimized cumulative infiltration curves, and their differences, are shown in Fig. 4. The maximum deviations were slightly higher thanfor experiment T1, about 2.5% of the cumulative infiltration volume. Notice
that, as compared with experiment T1, a much smaller volume of water
infiltrated. The lower infiltrated volume caused the saturated hydraulic conductivity K s of the surface crust to be two orders of magnitude lower than for the
subsoil. Also the water content interval between us and ur was now only about
0.24. This result corresponds well with observations of Valentin and Bresson
Ž1992. and Casenave and Valentin Ž1992. that the porosity of an erosion crust is
generally very low. Casenave and Valentin Ž 1992. estimated the infiltration rate
into erosion-type surface crusts to be between 0.0 to 4.8 cm dy1, with ourestimated value Ž 0.15 cm dy1 . being in this interval. The value of n was morecharacteristic for a clay loam textural group, while a and K s were typical of asilty clay ŽCarsel and Parrish, 1988. .
Table 3 presents results of Wooding’s analysis for experiment T2. Note thatthe ‘steady-state’ infiltration rates increased, as expected, when the supply
tension decreased to 3 cm. After that the infiltration rates actually decreasedslightly. This result likely means that the air-entry value was reached at apressure head of about y3 cm. For infiltration at lower tensions, the hydraulicconductivity remains then equal to K s and only changes in the pressure headgradient can influence the infiltration rate. The resulting values for both K s anda ) fluctuated several orders of magnitude between the different tension intervals, may be due to the underlying assumptions necessary for correct applicationof the Wooding’s solution. In one case K s even reached a negative value,
Table 6
Soil hydraulic parameters for surface crustParameters

Value

NSD w%xur
usa wcmy1 xn
K s wcm sy1 x
Ff

0.0
0.245
0.000411
1.35
0.0000017
0.00227

y
2.51
11.4
8.98
27.1

78

ˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

Fig. 4. Experiment T2: measured and calculated cumulative infiltration curves Ža. and theirdifferences Žb. for the soil profile with the surface crust.

further reflecting that Wooding’s analysis is inappropriate for a two-layered soil
profile.
5. Conclusions
We analyzed cumulative infiltration data from tension-disc infiltrometer byˇ ˚
numerical inversion ŽSimunek
and van Genuchten, 1997. and using the analytical solution of Wooding Ž1968. . We first analyzed the infiltration experimentcarried out on the subsoil only, from which we obtained its soil hydrauliccharacteristics. We subsequently analyzed the infiltration experiment for a
two-layered system—crust and subsoil—with the assumption that the hydrauliccharacteristics of the subsoil are known from the previous analysis, thus leading
to estimates of the hydraulic parameters of the surface crust. The hydraulicˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

79

parameters for the subsoil obtained by parameter estimation and Wooding’s
analytical method compared well with each other, and also with parameters
obtained by means of a neural network prediction ŽSchaap et al., 1998. based on
available soil textural information.
The numerical inversion technique allowed us to obtain the soil hydraulicproperties of both the surface crust and the subsoil. Since direct measurementsof the soil hydraulic properties of the surface crust are extremely difficult toobtain, no information was available for direct comparison with our numericalresults. The saturated hydraulic conductivity K s of the surface crust was foundto be two orders of magnitude lower than for the subsoil.
While many laboratory and field methods exist Ž Klute, 1986. to determine thesoil water retention and unsaturated hydraulic conductivity curves, most methods remain relatively expensive and too cumbersome for applications at relatively larger scales. Parameter estimation techniques provide an attractive alternative to more traditional, time consuming steady state methods which often are
also limited to a relatively narrow range in the water content. Time-savings are
possible since transient flow experiments generally cover a much broader rangeof water contents, while no need exists to reach flow equilibrium. The savings intime, however, are sometimes accompanied by a loss of accuracy since not alltransient experiments will provide enough information to guarantee a uniqueˇ ˚
solution ŽSimunek
and van Genuchten, 1997. . To avoid or limit such a loss of
accuracy, detailed analyses should be carried out before a particular transient
flow experiment is used to estimate the soil hydraulic properties Ž e.g., Toorman
et al., 1992.. An even faster and economical way of estimating the soil hydraulicparameters would be to use indirect methods, such as regression or neuralnetwork analyses, which relate the hydraulic parameters to more easily measured soil taxonomic data such as texture and bulk density ŽSchaap et al., 1998. .
The additional savings in cost and time of indirect methods are generally offset
by an even larger loss of prediction accuracy Ž see Table 4..

References
Ankeny, M.D., Ahmed, M., Kaspar, T.C., Horton, R., 1991. Simple field method for determining
unsaturated hydraulic conductivity. Soil Sci. Soc. Am. J. 55, 467–470.
Ahuja, L.R., 1974. Applicability of the Green–Ampt approach to water infiltration through surface
crust. Soil Sci. 118, 283–288.
Baumer, O.W., 1992. Predicting unsaturated hydraulic parameters. In: van Genuchten, M.Th.,
Leij, F.J., Lund, L.J. ŽEds.., Indirect Methods for Estimating the Hydraulic Properties of
Unsaturated Soils. Proceedings of the International Workshop, Univ. California, Riverside,CA, pp. 341–354.
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water retention
characteristics. Water Resour. Res. 24, 755–769.
Casenave, A., Valentin, C., 1992. A runoff capability classification system based on surface
features criteria in semi-arid areas of West Africa. J. Hydrol. 130, 231–249.

80

ˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

Clausnitzer, V., Hopmans, J.W., 1995. Non-linear parameter estimation: LM_OPT. General-purpose optimization code based on the Levenberg–Marquardt algorithm. Land, Air and Water
Resources Paper No. 100032, University of California, Davis, CA.
Cuenca, R.H., Brouwer, J., Chanzy, A., Droogers, P., Galle, S., Gaze, S.R., Sicot, M., Stricker,
H., Angulo-Jaramillo, R., Boyle, S.A., Bromly, J., Chebhouni, A.G., Cooper, J.D., Dixon, A.J.,
Fies, J.-C., Gandah, M., Gaudu, J.-C., Laguerre, L., Soet, M., Steward, H.J., Vandervaere,
J.-P., Vauclin, M., 1997. Soil measurements during HAPEX-Sahel intensive observation
period. J. Hydrol. 188–189, 224–266.
Gardner, W.R., 1958. Some steady-state solutions of the unsaturated moisture flow equation withapplication to evaporation from a water table. Soil Sci. 85, 228–232.
Goutorbe, J.-P., Lebel, T., Tinga, A., Bessemoulin, P., Brouwer, J., Dolman, A.J., Engman, E.T.,
Gash, J.H.C., Hoepffner, M., Kabat, P., Monteny, B.A., Prince, S., Said, F., Sellers, P.,
Wallace, J.S., 1994. HAPEX-Sahel: a large scale study of land–atmosphere interactions in thesemi-arid tropics. Ann. Geophys. 12, 53–64.
Haverkamp, R., Ross, P.J., Smettem, K.R.J., Parlange, J.Y., 1994. Three-dimensional analysis of
infiltration from the disc infiltrometer: 2. Physically based infiltration equation. Water Resour.
Res. 30, 2931–2935.
Hillel, D., Gardner, W.R., 1969. Steady infiltration into crust-topped profiles. Soil Sci. 108,
137–142.
Hillel, D., Gardner, W.R., 1970. Transient infiltration into crust-topped profiles. Soil Sci. 109,
69–74.
Hussen, A.A., Warrick, A.W., 1993. Alternative analysis of hydraulic data from the disc tension
infiltrometers. Water Resour. Res. 29, 4103–4108.
Jarvis, N.J., Messing, I., 1995. Near-saturated hydraulic conductivity in soils of contrasting texture
measured by tension infiltrometers. Soil Sci. Soc. Am. J. 59, 27–34.
Kern, J.S., 1995. Evaluation of soil water retention models based on basic soil physical properties.
Soil Sci. Soc. Am. J. 59, 1134–1141.
Klute, A., 1986. Water retention: laboratory methods. In: Klute, A. ŽEd.., Methods of Soil
Analysis, Part 1—Physical and Mineralogical Methods, 2nd edn. Agronomy Monograph 9, pp.
635–662.
Kool, J.B., Parker, J.C., van Genuchten, M.Th., 1985. Determining soil hydraulic properties fromone-step outflow experiments by parameter estimation: I. Theory and numerical studies. Soil
Sci. Soc. Am. J. 49, 1348–1354.
Kool, J.B., Parker, J.C., van Genuchten, M.Th., 1987. Parameter estimation for unsaturated flow
and transport models—a review. J. Hydrol. 91, 255–293.
Leij, F.J., Alves, W.J., van Genuchten, M.Th., 1996. The UNSODA unsaturated soil hydraulicdatabase. EPA, Ada, OK.
Leij, F.J., Russell, W.B., Lesch, S.M., 1997. Closed-form expressions for water retention andconductivity data. Ground Water 35 Ž5., 848–858.
Logsdon, S.D., Jaynes, D.B., 1993. Methodology for determining hydraulic conductivity withtension infiltrometers. Soil Sci. Soc. Am. J. 57, 1426–1431.
Logsdon, S.D., McCoy, E.L., Allmaras, R.R., Linden, D.R., 1993. Macropore characterization byindirect methods. Soil Sci. 155, 316–324.
Marquardt, D.W., 1963. An algorithm for least-squares estimation of nonlinear parameters. SIAM
J. Appl. Math. 11, 431–441.
Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated porous
media. Water Resour. Res. 12, 513–522.
Pachepsky, Y.A., Timlin, D., Varallyay, G., 1996. Artificial neural networks to estimate soil water
retention from easily measurable data. Soil Sci. Soc. Am. J. 60, 727–733.
Perroux, K.M., White, I., 1988. Design for disc permeameters. Soil Sci. Soc. Am. J. 52,
1205–1215.

ˇ ˚ et al.r Geoderma 86 (1998) 61–81
J. Simunek

81

Philip, J.R., 1969. Theory of infiltration. Adv. Hydrosci. 5, 215–296.
Reynolds, W.D., Elrick, D.E., 1991. Determination of hydraulic conductivity using a tension
infiltrometer. Soil Sci. Soc. Am. J. 55, 633–639.
Schaap, M.G., Bouten, W., 1996. Modelling water retention curves of sandy soils using neuralnetworks. Water Resour. Res. 32, 3033–3040.
Schaap, M.G., Leij, F.J., van Genuchten, M.Th., 1998. Neural network analysis for hierarchicalprediction of soil water retention and saturated hydraulic conductivity. Soil Sci. Soc. Am. J., inpress.
Šimunek,J., van Genuchten, M.Th., 1996. Estimating unsaturated soil hydraulic properties from˚
tension-disc infiltrometer data by numerical inversion. Water Resour. Res. 32, 2683–2696.
Šimunek,J., van Genuchten, M.Th., 1997. Parameter estimation of soil hydraulic properties from˚
multiple tension disc infiltrometer data. Soil Sci. 162, 383–398.
ˇ ˚
ˇ
Simunek,J., Sejna,M., van Genuchten, M.Th., 1996. The HYDRUS-2D software package forsimulating water flow and solute transport in two-dimensional variably saturated media.
Version 1.0. IGWMC-TPS-53, International Ground Water Modeling Center, Colorado Schoolof Mines, Golden, CO.
Smettem, K.R.J., Clothier, B.E., 1989. Measuring unsaturated sorptivity and hydraulic conductivity using multiple disk permeameters. J. Soil Sci. 40, 563–568.
Smilles, D.E., Knight, J.H., Perroux, K.M., 1982. Absorption of water by soil: the effect of asurface crust. Soil Sci. Soc. Am. J. 46, 476–481.
Tamari, S., Wosten,
J.H.M., Ruiz-Suarez,
J.C., 1996. Testing an artificial neural network for¨´
predicting soil hydraulic conductivity. Soil Sci. Soc. Am. J. 60, 1732–1741.
Thony, J.-L., Vachaud, G., Clothier, B.E., Angulo-Jaramillo, R., 1991. Field measurement of thehydraulic properties of soil. Soil Technol. 4, 11–123.
Tietje, O., Hennings, V., 1996. Accuracy of the saturated hydraulic conductivity prediction bypedo-transfer functions compared to the variability within FAO textural classes. Geoderma 69,
71–84.
Tietje, O., Tapkenhinrichs, M., 1993. Evaluation of pedotransfer functions. Soil Sci. Soc. Am. J.
57, 1088–1095.
Toorman, A.F., Wierenga, P.J., Hills, R.G., 1992. Parameter estimation of hydraulic properties
from one-step outflow data. Water Resour. Res. 28, 3021–3028.
Valentin, C., Bresson, L.-M., 1992. Morphology, genesis and classification of surface crusts inloamy and sandy soils. Geoderma 55, 225–245.
van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity ofunsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.
van Genuchten, M.Th., 1981. Non-equilibrium transport parameters from miscible displacement
experiments. Research Report No. 119. U.S. Salinity Laboratory, USDA, ARS, Riverside, CA.
van Genuchten, M.Th., Leij, F.J., Yates, S.R., 1991. The RETC code for quantifying the hydraulicfunctions of unsaturated soils. EPAr600r2-91-065, U.S. EPA, Office of Research andDevelopment, Washington, DC.
Vandervaere, J.-P., Peugeot, C., Vauclin, M., Angulo-Jaramillo, R., Lebel, T., 1997. Estimating
hydraulic conductivity of crusted soils using disc infiltrometers and minitensiometers. J.
Hydrol. 188–189, 203–223.
Warrick, A.W., 1992. Model for disc infiltrometers. Water Resour. Res. 28, 1319–1327.
White, I., Sully, M.J., 1987. Macroscopic and microscopic capillary length and time scales fromfield infiltration. Water Resour. Res. 23, 1514–1522.
Williams, R.D., Ahuja, L.R., Naney, J.W., 1992. Comparison of methods to estimate soil water
characteristics from limited texture, bulk density, and limited data. Soil Sci. 153, 172–184.
Wooding, R.A., 1968. Steady infiltration from a shallow circular pond. Water Resour. Res. 4,
1259–1273.