WATER RESOURCES RESEARCH, VOL. 35, NO. 5, PAGES 1329 –1345, MAY 1999

Estimating hysteresis in the soil water retention function
from cone permeameter experiments
Jiřı́ Šimůnek
U.S. Salinity Laboratory, ARS, USDA, Riverside, California

Radka Kodešová and Molly M. Gribb
Department of Civil and Environmental Engineering, University of South Carolina, Columbia

Martinus T. van Genuchten
U.S. Salinity Laboratory, ARS, USDA, Riverside, California

Abstract. Data obtained from modified cone penetrometer experiments were used toestimate the hysteretic soil hydraulic properties with a parameter estimation techniquewhich combined a numerical solution of the Richards equation with Marquardt-Levenberg
optimization. The modified cone penetrometer was designed to inject water into a soil
through a cylindrical screen, measure the infiltration rate with time, and track the
movement of the wetting front using two tensiometer rings positioned above the screen.
After reaching relatively stable tensiometer readings during the experiments, the source ofwater was cut off and pressure head readings measured while water in the soil profile
redistributed. Cumulative inflow and pressure head readings for two experiments withdifferent supply pressures were analyzed to obtain estimates of the soil water retention
and hydraulic conductivity functions. Analysis of flow responses obtained during theinfiltration period, and of those obtained during the combined infiltration andredistribution phases, demonstrated the importance of hysteresis of the soil hydraulicfunctions. We found that the redistribution phase could not be described accurately whenhysteresis was neglected. Hysteresis in the soil hydraulic functions was modeled using a
relatively simple empirical model in which wetting scanning curves are scaled from themain wetting curve and drying scanning curves are scaled from the main drying curve.
This model was deemed adequate for our examples. Optimization results for various
combinations of unknown soil hydraulic parameters were compared to results of standard
laboratory and in situ methods. Estimates of the saturated hydraulic conductivity were wellwithin the range of in situ measurements. The estimated main hysteretic loops of the soil
water retention curve were for the most part situated between the wetting and drying
curves obtained with standard methods.

1.

Introduction

Hysteresis can significantly influence water flow and solutetransport in variably saturated porous media [e.g., Vachaud
and Thony, 1971; Royer and Vachaud, 1975; Gillham et al.,
1979; Kaluarachchi and Parker, 1987; Jaynes, 1992]. The phenomenon has been observed under both laboratory [e.g., Topp,
1969; Gillham et al., 1979] and field [e.g., Vachaud and Thony,
1971; Royer and Vachaud, 1975; Watson et al., 1975] conditions.
The significance of hysteresis has also been demonstrated inseveral numerical studies [e.g., Gillham et al., 1979; Kaluarachchi and Parker, 1987; Russo et al., 1989]. Although longrecognized as being important, hysteresis is usually neglectedin water flow studies because of a lack of good data.
Hysteresis refers to the nonunique relationship between thepressure head, h, and the water content, u, in the soil water
retention function u (h) (Figure 1). This relationship displays
considerable variations in u for the same h depending upon the
Copyright 1999 by the American Geophysical Union.
Paper number 1998WR900110.
0043-1397/99/1998WR900110$09.00

history of soil wetting and drying. During infiltration when thewater content and the pressure head are monotonically increasing, the retention curve can be described by a uniquefunction. Similarly, during evaporation or gravity drainage
when the water content and the pressure head are monotonically decreasing, the retention curve can also be described by aunique but different function. When drying reverses into awetting process and vice versa, u (h) is no longer unique andhysteresis must be introduced [Russo et al., 1989]. When thesoil is wetted from the residual water content or drained fromthe saturated water content, the main wetting or drying curves
are followed, respectively. When a wetting or drying process isreversed while following the main hysteresis curve, the retention curve follows a primary hysteresis curve. Secondary andhigher-order scanning curves are a result of additional reversals [Jaynes, 1992; Kool and Parker, 1987]. The hysteresis process is described in detail by Mualem [1974], Kool and Parker[1987], Luckner et al. [1989], and Jaynes [1992].
Numerous models for describing hysteresis in the soil water
retention curve have been developed. Two major groups ofmodels can be distinguished: physically based models and empirical models. Domain models are the most common theoret-

1329

1330

ŠIMŮNEK ET AL.: ESTIMATING HYSTERESIS

Figure 1. Idealized hysteretic soil water retention function
showing main, primary, and secondary wetting and drying
curves.

ical models. These models divide the soil matrix into a set ofpore domains characterized by wetting pore radii that control
the filling of pores by water and drying pore radii that control
their emptying. A similarity assumption for wetting and drying
pore radii [Mualem, 1973, 1974] is often invoked to avoid therequirement for a large number of scanning curves for model
calibration. With the similarity assumption, only the main wetting and main drying curves are needed to describe the entire
hysteretic phenomenon. Depending on whether pore drainage
is assumed to be independent or dependent upon the state ofneighboring pores, independent-domain models [Mualem,
1973, 1974] and dependent-domain models [Mualem, 1984a;
Mualem and Dagan, 1975] have been developed.
Empirical analytical models assume that the primary, secondary, and higher-order scanning curves can be scaled fromthe main hysteresis curve [Klute and Heerman, 1974; Hoa et al.,
1977; Scott et al., 1983; Kool and Parker, 1987]. Different analytical expressions describing the soil water characteristic curve
were used in these studies. Kool and Parker [1987] coupled theanalytical model of van Genuchten [1980] with the simplified
scaling approached used by Scott et al. [1983] to describe thescanning curves. Scott et al. [1983] assumed that the shapeparameters for all drying scanning curves are the same as thosefor the main drying curve and similarly that the shape parameters for all wetting scanning curves are the same as those forthe main wetting curve. Scanning curves are then calculated byvarying the residual and saturated water contents for the wetting and drying scanning curves, respectively. Kool and Parker[1987] also assumed that one of the two shape parameters isthe same for both wetting and drying, thus decreasing thenumber of required parameters. An advantage of the model by
Scott et al. [1983] is that one can obtain the shape parametersof the drying and wetting curves from knowledge of any main(or primary), secondary, or higher-order scanning drying andwetting curves, respectively.
Most methods for measuring the soil hydraulic properties
remain relatively time consuming and costly, are often limited
to a relatively narrow range of water contents, and usuallypertain to relatively restrictive initial and boundary conditions.
Measurement of soil water hysteresis is even more time con-suming since both water contents and pressure heads must be
measured for at least the main drying and wetting curves, andpreferably also for additional primary, secondary, and higherorder scanning curves. Generation of a data set suitable forcalibrating the domain model of hysteresis can require several
months or even years, depending upon whether knowledge of
only the main hysteresis loop is needed, or also additional
primary and secondary scanning curves. Therefore attempts
have been made to develop domain models that require only
one branch of the retention curve for calibration of hysteresis.
These models, however, have not been found to be very successful [Mualem and Morel-Seytoux, 1978].
A more expedient way for obtaining information about hysteresis in the soil hydraulic properties is to design appropriate
transient experiments involving hysteresis and to analyze thedata with a suitable inverse procedure. The effort of determining hydraulic properties can thus be shifted from experimentation to computation [Kool and Parker, 1988]. A transient
experiment is first carried out with known initial and boundary
conditions and subsequently simulated using the governing
water flow equation and selected analytical functions representing the soil hydraulic properties. Although the analytical
functions must be selected in advance, their parameters are
generally unknown. Starting with the initial parameter estimates, the parameters are further adjusted in repeated solutions until deviations between measured and computed waterflow attributes (such as water contents, pressure heads, or
fluxes) are minimized. The parameter updates can be done
manually (trial-and-error calibration) or by using an automated minimization algorithm.
Inverse methods in subsurface hydrology were initially usedalmost exclusively for saturated flow problems (see the review
by Yeh [1986]); they are now also often utilized to analyze
unsaturated zone experiments (see the review by Kool et al.
[1987]). Parameter estimation techniques were first used withlaboratory experiments, such as one-step-type or multisteptype outflow methods [Kool et al., 1985; van Dam, 1992, 1994;
Eching and Hopmans, 1993] and evaporation experiments [Ciollaro and Romano, 1995; Santini et al., 1995; Šimůnek et al.,1998], for which the initial and boundary conditions and homogeneity of the soil sample can be readily controlled. The
application of inverse methods to field experiments, such as theinstantaneous profile method [Dane and Hruska, 1983], disc
infiltrometers [Šimůnek and van Genuchten, 1996, 1997], or the
extraction method [Inoue et al., 1998], is more complex becauseof soil heterogeneity and uncertainty associated with theboundary conditions. In addition, field methods have been
generally limited to the soil surface or the near surface environment.
A new cone permeameter method currently under development [Gribb, 1996; Gribb et al., 1998; Kodešová et al., 1998a, b]
for estimating soil hydraulic properties has the potential foruse up to depths of 30 m or more. Cone penetrometers wereoriginally used to obtain soil strength characteristics by measuring the tip resistance and sleeve friction during penetration
at a constant rate. To obtain the hydraulic properties, a modified cone penetrometer, instrumented with a porous filter
close to the penetrometer tip and two tensiometer rings above
the filter, is used. The device is pushed into a soil to the desired
depth, and a constant head is applied to the filter. The volume
of water imbibed into the soil is monitored during infiltration.
The tensiometer rings register movement of the wetting front

ŠIMŮNEK ET AL.: ESTIMATING HYSTERESIS

during infiltration, as well as the subsequent redistribution
process after the source of water has been turned off.
Gribb [1996] gave a detailed numerical analysis of this experiment, including a study of the identifiability of the soil
hydraulic parameters. She showed that the inverse solution was
most sensitive to the saturated hydraulic conductivity K s andthe shape parameter a and least sensitive to the saturated
water content u s and the parameter n in van Genuchten’s[1980] retention function (see equation (10)). The method wasrecently used to estimate K s and the wetting branch of theretention curve of a sandy soil in a laboratory aquifer system[Gribb et al., 1998; Kodešová et al., 1998a, b]. Gribb et al. [1998]
also showed that the cone permeameter can be used in fully
saturated media to estimate K s . Kodešová et al. [1998b] evaluated the effects of cone placement in a soil on the optimized
soil hydraulic properties. One set of tests was performed after
the cone permeameter was buried in the soil and a second set
after the permeameter was pushed to the desired depth. They
concluded that for the sandy soil tested, the method of placement had only little impact on the optimized hydraulic parameter values. They also showed that inclusion of the final watercontent data in the optimization problem improved estimatesof u s . Only the wetting parts of the infiltration experiments
were considered in the above mentioned studies.
In this paper we analyze data from both the infiltration andredistribution parts of the experiments utilized by Kodešová etal. [1998b] so as to estimate the hydraulic parameters of both
the wetting and drying branches of the soil water retention
characteristic. We use only data from two tests performed withthe prototype instrument after it was buried (not pushed) inthe soil; the tests involved different applied pressure heads but
similar initial conditions. The two tests were selected in order
to avoid any disturbance of the soil profile caused by conepushing with concomitant changes in the estimated soil hydraulic parameters. Kodešová et al. [1998b] showed that theobserved flow responses could be fit very well with physicallyreasonable estimates of all optimized hydraulic parameters,with the exception of u s which was too high for the two examples. For cases in which u s was fixed at the laboratory-derived
value or was reasonably optimized using additional water content data, poorer fits of the measured flow responses wereobtained. To improve the solution, we first investigate here theeffects of including additional optimized parameters, such as
the anisotropy coefficient k A and the pore connectivity parameter l (see equations (3) and (11), respectively), while considering only the infiltration parts of the experiments. The sameprocedure is then applied to the infiltration and redistribution
parts of the experiments with and without consideration ofhysteresis. The resulting soil hydraulic characteristics are compared with independently measured laboratory values.

2.

Theory2.1. Cone Permeameter Experimental Setup: Problem
Definition
A cone penetrometer was modified for hydraulic testing.
The prototype was discussed in detail by Gribb et al. [1998].
The device consists of four parts: a shaft, a 5-cm screen, andtwo porous ceramic rings serving as tensiometers located 5 and9 cm above the screened section (Figure 2). After the cone
permeameter is installed in a soil, a constant head is applied tothe screen using a microprocessor-controlled solenoid valve
assembly. Cumulative flow volume infiltration into the soil

1331

Figure 2. Schematic of the cone permeameter.
through the screen is determined from scale readings of themass of water removed from a source bottle. The tensiometerrings are connected to pressure transducers to continuously
record the transient pore water pressure heads in the soil.
Injection of water continues until relatively stable tensiometerreadings are reached. The source of water is then cut off, andpressure heads are continuously measured during the ensuing
redistribution phase.
Cumulative inflow and pressure head readings are subsequently analyzed to obtain estimates of the soil water retention, u (h), and hydraulic conductivity, K(h), functions. The
unknown parameters of K(h) and u (h) are obtained by minimizing an objective function expressing the difference between measured flow variables and predicted system responses. The system response is modeled by numericalsolution of the appropriate governing variably saturated flow
equation, augmented with parameterized hydraulic functions
and suitable initial and boundary conditions. Initial estimatesof the optimized system hydraulic parameters are iteratively
improved during the minimization process until a desired degree of precision is obtained.
2.2.

Governing Variably Saturated Flow Equation

The governing equation for radially symmetric isothermal
Darcian flow in a variably saturated anisotropic rigid porous
medium is given by the following modified form of Richards’
equation:

S F
GD

­u
1 ­~rv r! ­v z 1 ­
52
2
5
­t
r ­r
­z
r ­r
1

­
­z

SF

K K Azj­h
1 K Azz­ xj

rK K Arj­h
1 K Arz
­ xj

GD
(1)

where r is the radial coordinate [L]; z is the vertical coordinate
positive upward [L]; t is time [T]; u is the volumetric watercontent [L 3 L 23 ]; v r and v z are volumetric fluxes in radial andvertical directions [L T 21 ]], respectively; h is the pressure
head [L]; K A
ij are components of a dimensionless anisotropy
tensor KA ; and K is the unsaturated hydraulic conductivity
function [L T 21 ] given by
K~h! 5 K sK r~h!

(2)where K r is the relative hydraulic conductivity [–] and K s thesaturated hydraulic conductivity in the vertical direction [L
T 21 ]. The anisotropy tensor K A
ij in (1) may be used to accountfor an anisotropic medium. The diagonal entries of K A
ij areequal to 1 and the off-diagonal entries are equal to zero for an

1332

ŠIMŮNEK ET AL.: ESTIMATING HYSTERESIS

voked unstructured triangular finite element mesh, is shown in
Figure 3.
The solution of (1) for conditions involving hysteresis in soil
hydraulic properties requires knowledge of the initial distribution of the pressure head within the flow domain,h~r, z, t! 5 h i~r, z!

t50

(6)

where h i is the initial pressure head [L], as well as additional
information characterizing the history status, that is, information whether the initial condition is on the main drying curve,
the main wetting curve, or on a particular scaling curve. The
initial water content must also be specified when the initialcondition is not on the main wetting or drying branch of theretention curve but inside of the main hysteresis loop.
The boundary condition along the screen section of the coneduring infiltration is as follows (Figure 3):h~r, z, t! 5 h 0~t! 2 ~ z 2 z 0!

Figure 3. Computational domain and finite element mesh.

r 5 r 0, z 0 , z , z 0 1 L(7)where h 0 is the supply pressure head imposed at the bottom of
the screen [L], z 0 is the coordinate of the bottom of the screen[L], L is the length of the screen [L], and r 0 is the radius of
the screen [L]. Notice that the supply pressure head, h 0 , canbe variable with time; h 0 for our experiment was kept constantduring the initial (first) stage of the infiltration process. The
cone, however, was still completely filled with water at the
moment when the source bottle was disconnected from thepermeameter. This water infiltrated during the second stage ofthe infiltration process, which required us to calculate thesupply pressure head from the remaining volume of water inthe cone, the geometric characteristics of the interior of thecone body, and the actual infiltration rate. During this second
stage of infiltration, the supply pressure head was calculated
with

S E D
t

h 0~t! 5 f~V~t!! 5 f V 0 2
isotropic medium. Since in the following we assume that the
principal directions of the anisotropy tensor KA coincide with
global coordinate system (r, z), the off-diagonal entries areequal to 0. The degree of anisotropy can then be represented
by an anisotropy coefficient, k A , which relates the first andsecond principal components of the anisotropy tensor KA as
follows [Mualem, 1984b]:k A 5 K Arr/K Azz(3)

Assuming that K A
zz 5 1, the anisotropy tensor can be writtenas
KA 5

U k0 10 U

Equation (1) simplifies to­u 1 ­
5
­t
r ­r

S

­h
rKk
­r
A

A

D

­
1
­z

(4)

S

­h
K
1K
­z

D

Notice that we assume that anisotropy is unaffected by the
degree of saturation; that is, (3) holds for both saturated andunsaturated soils.
Equation (5) subject to specified initial and boundary conditions is solved numerically using the finite element method
for space discretization and finite differences with a fully implicit scheme for temporal discretization [Šimůnek et al., 1996].
A schematic of the transport domain, together with the in-(8)

where V and V 0 are the volumes of water at the current timeand at the beginning of the second infiltration stage inside of
the cone [L 3 ], respectively; t 0 represents the end of the first
infiltration stage and start of the second stage [T]; and I is theactual infiltration rate [L 3 T 21 ]. During each time step the
supply pressure head at the screen was kept constant and theinfiltration rate corresponding with this Dirichlet-type boundary condition calculated. The volume of water in the cone
penetrometer was then correspondingly decreased, thus resulting in the supply pressure head for the new time step. Once allwater from the cone had infiltrated, the redistribution stage
begins. The following boundary condition then applies along
the screen:v r~r, z, t! 5 0(5)

I~t! dtt0

r 5 r 0, z 0 , z , z 0 1 L(9)Other no-flow boundaries are implemented far enough awayfrom the source such that they do not affect water flow during
the calculations.
2.3.

The Unsaturated Soil Hydraulic Properties

We implement the soil hydraulic functions of van Genuchten[1980], who used the statistical pore-size distribution model of
Mualem [1976] to obtain a predictive equation for the unsaturated hydraulic conductivity function in terms of soil water

ŠIMŮNEK ET AL.: ESTIMATING HYSTERESIS

retention parameters. The expressions of van Genuchten[1980] are given by

u ~h! 5H

ur 1

us 2 ur@1 1 u a hu n# m
us

h,0
h$0

m 2
K~h! 5 K sS le@1 2 ~1 2 S 1/m
e ! #

(10)
(11)

where S e is the effective water content [–]; K s is the saturated
hydraulic conductivity [L T 21 ]; u r and u s denote the residual
and saturated water contents [L 3 L 23 ], respectively; l is apore-connectivity parameter [–]; and a [L 21 ], n [–], and m(5 1 2 1/n) [–] are empirical parameters. The predictive
K(h) model is based on the capillary model of Mualem [1976]
in conjunction with (10). The above equations contain six independent parameters: u r , u s , a , n, l, and K s . The poreconnectivity parameter l in the hydraulic conductivity function
is usually assumed to be 0.5 as an average for many soils[Mualem, 1976].
We incorporated hysteresis by using the empirical model
introduced by Scott et al. [1983], which assumes that drying
scanning curves are scaled from the main drying curve and thatwetting scanning curves are scaled from the main wetting
curve. This model was also previously employed by Kool andParker [1987], who modified the formulation to account for air
entrapment, and by Vogel et al. [1996], who also consideredhysteresis in the hydraulic conductivity function. The adopted
procedure for modeling hysteresis in the retention function
requires that both the main drying and main wetting curves beknown. These two curves are described with (10) using thew
w
w
parameter vectors ( u dr , u ds , a d , n d ) and ( u w
r , u s , a , n ),where the superscripts d and w indicate drying and wetting.
The following restrictions are expected to hold in most practical applications:u r 5 u dr 5 u wr

ad # aw(12)ting cycles result in a considerably lower water content thanwould be predicted by the domain model of Mualem [1984a].
This pumping process is contrary to observed behavior [Jaynes,
1992]. However, “pumping” can become apparent only after
several reversals in the wetting and drying cycles and hence
should not be important in this study since our experiments
involve only one reversal from wetting to drying.
2.4.

Formulation of the Inverse Problem

The objective function F to be minimized during the parameter estimation process is formulated using cumulative infiltration data in combination with transient pressure head readings
from two tensiometers. The objective function for multiple
measurement sets is defined as [Šimůnek and van Genuchten,
1996]

d
s

w
s

d

(13)
w

so that the parameters u , u , a , and a are the only independent parameters describing hysteresis in the retention function.
If the main hysteresis loop is not closed at saturation, thewater content at saturation for a particular wetting scanning
curve can be evaluated using the empirical relationship of Aziz
and Settari [1979]:

u s 5 u ds 2
R5

u ds 2 u D
1 1 R~ u ds 2 u D!

11
2 d
u 2 u ws
u s 2 u dr

(14)

d
s

where uD represents the reversal point. In our analysis, however, we will assume that u s 5 u ds 5 u w
s , so that the hysteresis
retention model is characterized with five parameters ( u r , u s ,a d , a w , n). With additional three parameters characterizing
the conductivity model (l, k A , and K s ), the entire homogeneous system is fully described with eight independent parameters.
The above model of hysteresis may lead to a so called“pumping effect,” which occurs when several drying and wet-

OS Om

F~b, q 1, · · · , q m! 5

nj

vj

j51

w i, j@q *j~t i! 2 q j~t i, b!# 2

i51

D

(15)
where m represents the number of different sets of measurements (cumulative infiltration data and/or tensiometer readings), n j is the number of measurements in a particular set,q *j (t i ) are specific measurements at time t i for the jth measurement set, b is the vector of optimized parameters (e.g., u r ,u s , a d , a w , n, l, k A , and K s ), q j (t i , b) are the corresponding
model predictions for parameter vector b, and v j and w i, j areweights associated with a particular measurement set or point,respectively. We assume for now that the weighting coefficients
w i, j in (15) are equal to 1, that is, the variances of the errors
inside a particular measurement set are assumed the same.
The weighting coefficients v j , which minimize differences inweighting between different data types because of different
absolute values and numbers of data involved, are given by[Clausnitzer and Hopmans, 1995]

We also invoke the often assumed restriction [Kool and Parker,
1987]
n 5 nd 5 nw

1333

vj 5

1
n js 2j(16)

thus defining the objective function as the average weighted
squared deviation normalized by measurement variances s 2j .
Minimization of the objective function F is accomplished byusing the Levenberg-Marquardt nonlinear minimization
method [Marquardt, 1963]. This method was found to be very
effective and has become a standard in nonlinear least squares
fitting among soil scientists and hydrologists [van Genuchten,
1981; Kool et al., 1985, 1987; Šimůnek and van Genuchten,
1996].

3.

Methods

3.1.

The Laboratory Aquifer

Cone tests were carried out in a laboratory aquifer measuring 4.7 3 4.7 3 2.6 m. The aquifer material was a sandy soil
with occasional kaolin pockets, underlain by 20 cm of gravel.
The bulk density of undisturbed soil samples ranged from 1.65
to 1.69 g cm23. The soil hydraulic properties were determinedusing several standard methods. The wetting branch of theretention curve was obtained by capillary rise tests [Lambe,
1951] and a computer-automated extraction/sorption laboratory testing method [Znidarčič et al., 1991; Ray and Morris,
1994]. The drainage branch was measured with pressure platetests. The pressure plate and the computer-automated tests

1334

ŠIMŮNEK ET AL.: ESTIMATING HYSTERESISTable 1. Hydraulic Parameters Obtained From Various Laboratory Tests
Testa w , cm21

a d , cm21n [–]u r [–]

u s [–]

K s , cm s21

Capillary riseComputer-automatedPressure plate method
Constant head laboratory
Slug
Guelph permeameter
Cone permeameter(saturated conditions)

0.0864
0.139
zzz
zzz
zzz
zzz
zzz

zzz
zzz
0.045
zzz
zzz
zzz
zzz

3.60
2.17
1.61
zzz
zzz
zzz
zzz

0.008
0.042
0.008
zzz
zzz
zzz
zzz

0.329
0.345
0.35
zzz
zzz
zzz
zzz

zzz
zzz
zzz
0.003850.00725
0.00349
0.0134

were performed on undisturbed samples, while the capillaryrise tests were carried out on repacked samples. The soil hydraulic parameters for both branches of the retention curve
were obtained by fitting (10) to u (h) data using the RETC
code [van Genuchten et al., 1991]. Soil hydraulic parameters
determined with the different laboratory methods are presented in Table 1.
A series of drive tube samples was taken from the testingarea and subjected to laboratory constant head permeability
tests. Besides the laboratory tests, slug tests and Guelph permeameter tests were performed in the laboratory aquifer inother studies [Scaturo, 1993; Singleton, 1997]. Slug test data
were analyzed using the Bower and Rice [1976] equation, while
Guelph permeameter data were analyzed according themethod of Reynolds [1993]. The saturated hydraulic conductivity was also estimated using Hvorslev’s [1951] analyticalequation adopted for infiltration from a cone penetrometer
into a saturated soil [Gribb et al., 1998]. Table 1 summarizes K svalues obtained with the various methods.
3.2.

Cone Penetrometer Tests

The cone penetrometer was placed in a hole in the laboratory aquifer and the soil around the cone was backfilled suchthat the center of the screened section was 65 cm below the soil
surface. This method of placement ensured relative homogeneity of the soil surrounding the prototype without possible soil
disturbance due to pushing the cone to the testing depth. The
water table in the laboratory aquifer was raised to groundsurface and then lowered first to 26 and then to 48 cm below
the surface of the aquifer. Under these conditions the cone wasused as a piezometer to obtain saturated hydraulic conductivity
values [Gribb et al., 1998]. The water table was subsequently
lowered to approximately 190 cm below the soil surface toestablish unsaturated conditions in the aquifer. This depth
ensured that tests would not be influenced by the presence ofa water table. The laboratory aquifer was left undisturbed for
2 weeks before running the first infiltration test. The second
infiltration test was performed about 3 weeks after the first one
using the same experimental setup. We selected these two
infiltration tests for further analysis.
During experiment 1 a supply pressure (h 0 ) of 52.5 cm wasapplied at the bottom of the screen for a period of 400 s (0 ,t , 400 s), while a pressure of 32.5 cm was applied for 445 s
during the experiment 2. The total volumes of water infiltrated
during experiments 1 and 2 were 11.2 and 6.95 l, respectively.
Tensiometer readings were obtained for 840 s (0 , t , 840s). Initial pressure heads of the upper tensiometer were approximately equal to 250 cm for both tests. This relatively low
pressure head corresponds to less than 10% of the effective
saturation value for this particular soil; the experiment hence
covered more than 90% of the range of saturation. Pressure

heads of the upper tensiometer at the end of the infiltration
period were 216 and 218 cm and dropped to about 234 and232 cm at the end of the experiments. The data collected forboth tests are shown in Figure 4. Although measurements wererecorded every second, only data collected at 5-s intervals were
included in objective function (15).
3.3.

Numerical Inversions

Table 2 summarizes the different numerical inversions carried out on both data sets. The first five inversions (set 1, runs
1–5) analyzed only the first infiltration stage for both experiments to obtain estimates of the wetting branches of soil hydraulic properties. Hysteresis was not considered in these runs.
In the first five inversions we optimized the parameters a w , n,
and K s alone or together with different combinations of theparameters u s , l, and k A . Using numerically generated data,
Gribb [1996] showed that the inversion procedure is very insensitive to the saturated water content u s , which is the reason
why we fixed this parameter in several optimizations at a valueof 0.35. This value was selected on the basis of calculations
from the bulk density and results of laboratory tests (Table 1).
Unless optimized, the parameter l was assumed to be equal to
0.5 [Mualem, 1976], while the anisotropy coefficient k A was set
equal to 1.0. The residual water content was always equated to
the value obtained from capillary rise experiments ( u r 5
0.008), since we did not expect this parameter to be identifiable from the near-saturated experiments. The first five inversions (set 1) were repeated using measured data from theentire experiments, that is, using both the infiltration and redistribution data. Hysteresis again was neglected in these optimizations.
The next eight optimizations (set 2, runs 6 –13) involved data
entirely from experiments 1 and 2 and considered hysteresis in
the soil hydraulic properties. These data sets contained information about both wetting and drying and about the corresponding branches of the soil hydraulic properties. Again, theparameters a w , n, and K s were optimized in all eight inversions. The parameter a d was either optimized independently
or assumed coupled with a w by the relationshipaw 5 2ad(17)

which is often used as a reasonable approximation [Kool andParker, 1987; Nielsen and Luckner, 1992]. Again, the parameters u s , l, and k A , with or without a combination of u s and lwere optimized in the numerical inversions. Parameter constraints had to be specified to avoid numerical instabilities forhigh values of n and l (n [ [1.1, 7.5], l [ [220, 120]).
Constraints were also imposed on u s ( u s # 0.6) to avoid
physically unrealistic values. The initial estimates of the soil
hydraulic parameters in most inversions were as follows: u s 5

ŠIMŮNEK ET AL.: ESTIMATING HYSTERESIS

1335

Figure 4. Experimental data for (a) experiment 1 and (b) experiment 2.
0.35, a d 5 0.3, a w 5 0.5, n 5 4.0 (6.0), l 5 0.5, k A 5
1.0, and K s 5 0.01. Variances of all measured values ofpressure heads and cumulative infiltration volumes were used
as measurement variances in (16).

Depending upon the wetting or drying history of a soil, awide range of water contents can be associated with a particular measured pressure head. If no information about both theinitial water contents and the initial pressure heads is available,

Table 2. Summary of Optimization Runs
Run

Hysteresis

us

aw

1
2
3
4
5

no
no
no
no
no

0.35
opt
0.35
0.35
opt

opt
opt
opt
opt
opt

6
7
8
9
10
11
12
13

yes
yes
yes
yes
yes
yes
yes
yes

opt
opt
0.35
0.35
0.35
0.35
opt
opt

opt
opt
opt
opt
opt
opt
opt
opt

Numbers represent fixed parameters. Opt, optimized.
aOptimized as a w 5 2 a d .

Set 1

Set 2

n

Ks

l

kA

ad

opt
opt
opt
opt
opt

opt
opt
opt
opt
opt

0.5
0.5
0.5
opt
opt

1.0
1.0
opt
1.0
1.0

NA
NA
NA
NA
NA

opt
opt
opt
opt
opt
opt
opt
opt

opt
opt
opt
opt
opt
opt
opt
opt

0.5
0.5
0.5
0.5
opt
opt
opt
opt

1.0
1.0
opt
opt
1.0
1.0
1.0
1.0

opt
opta
opt
opta
opt
opta
opt
opta

1336

ŠIMŮNEK ET AL.: ESTIMATING HYSTERESISTable 3. Results of Numerical Inversions for Experiment 1 Without Considering HysteresisRunF

u s [–]

a w , cm21l [–]k A [–]

a d , cm21

Rank

1
2
3
4
5

0.06040.008770.00879
0.00956
0.00839

zzz
0.459
zzz
zzz
0.405

0.045
0.05200.05270.03570.0442

zzz
zzz
zzz
14.2
2.28

zzz
zzz
1.37
zzz
zzz

zzz
zzz
zzz
zzz
zzz

5
2
3
4
1

1*
2*
3*
4*
5*

0.09860.07770.07080.08200.0475zzz
0.600azzz
zzz
0.600aFitted to Infiltration and Redistribution Data
0.03494.02
0.0115zzz
0.02972.34
0.0122zzz
0.03292.42
0.00764zzz
0.0288
5.41
0.013
10.6
0.02683.22
0.0118
4.65

zzz
zzz
1.96
zzz
zzz

zzz
zzz
zzz
zzz
zzz

5
3
2
4
1

n [–]

K s , cm s21

Fitted to Only Infiltration Data
5.76
0.01146.77
0.01196.11
0.00965
6.05
0.01195.42
0.0119aUpper parameter constraint.

then this possible variation in water content must be considered in the optimization process. Therefore, set 2 inversions
were carried out twice. We first assumed that the initial condition was derived from the main wetting branch of the retention curve and next from the main drying branch. In general, asoil will be at a scanning curve in between the main wetting andmain drying branches. In this way, both extremes of the initialconditions are hence considered.

4.

Results and Discussion

4.1.

Analysis of the Infiltration Process

Results of the numerical inversions for experiments 1 and 2,
without considering hysteresis, are summarized in Tables 3 and4, respectively. The first five runs in each table used only
infiltration data, while the second set of runs (runs 1* through5*) used both infiltration and redistribution data. In addition
the optimized parameters, Tables 3 and 4 give also the final
values of the objective function and a ranking for each run (1being the best). The rankings presented here are based solely
on the final value of the objective function and do not consider
the number of optimized parameters, the accuracy of theirprediction, or any other criteria involving uniqueness of the
optimization results. Hence inversions with a higher number of
optimized parameters are ranked higher than those with fewer
optimized parameters. We emphasize that the rankings in Tables 3– 8 serve only to simplify interpretation of the results.
Parameters not given for a particular numerical inversion (e.g.,u r , u s , a d , l, k A ) were fixed and are given in Table 2.

Measured and optimized cumulative infiltration volumes
and pressure heads resulting from the first five numerical inversions are plotted in Figure 5. Measured pressure heads andcumulative infiltration volumes were not well described whenu s was fixed and only three parameters ( a w , n, and K s ) were
optimized (run 1). The model predicted an earlier arrival of
the wetting front at both tensiometers than actually observed
and also underestimated the cumulative infiltration volume forboth experiments. Improved fits could be obtained by optimizing u s , the anisotropy factor k A , and/or the tortuosity factor l(runs 2–5 in Tables 3 and 4). Results show that optimization ofu s , k A , and l, or u s and l, along with K s , a , and n yieldedalmost identical flow responses (Figure 5). The very similar
final values of the objective function (Tables 3 and 4) for runs
2–5 make it impossible to clearly distinguish between the effects of the different parameters. The estimated value of u s(run 2) was greater than 0.35 (the value estimated from laboratory soil bulk density and water content measurements) andprovided increased infiltration capacity. The optimized valueof the anisotropy coefficient k A (run 3) was greater than 1 (thevalue associated with isotropic flow), leading to more flow inthe radial direction and, indirectly, also more infiltration. The
fitted value of the pore-connectivity parameter l (run 4) was
greater than the 0.5 and resulted in a faster decrease of thehydraulic conductivity with decreasing water content. This feature yielded a much sharper wetting front. Finally, the simultaneous optimization of u s and l (run 5) generated values thatwere greater than 0.35 and 0.5, respectively, but lower thanvalues resulting from runs 2 and 4, respectively.

Table 4. Results of Numerical Inversions for Experiment 2 Without Considering HysteresisFu s [–]

a w , cm21

1
2
3
4
5

0.04600.005880.00551
0.00456
0.00456

zzz
0.430
zzz
zzz
0.350

0.05010.05700.05770.0486
0.0486

1*
2*
3*
4*
5*

0.10070.07540.07770.04470.0260

zzz
0.497
zzz
zzz
0.848

Run

a

Upper parameter constraint.

l [–]k A [–]

a d , cm21

Rank

zzz
zzz
zzz
2.71
2.69

zzz
zzz
1.27
zzz
zzz

zzz
zzz
zzz
zzz
zzz

5
4
3
1
2

Fitted to Infiltration and Redistribution Data
0.04064.74
0.00889zzz
0.03963.00
0.00966zzz
0.04103.01
0.00741zzz
0.02885.00
0.0091219.5
0.02082.75
0.00976
20.0a

zzz
zzz
1.48
zzz
zzz

zzz
zzz
zzz
zzz
zzz

5
3
4
2
1

n [–]

K s , cm s21

Fitted to Only Infiltration Data
6.36
0.00900
7.50a
0.00951
7.50a
0.008106.71
0.00952
6.72
0.00952

ŠIMŮNEK ET AL.: ESTIMATING HYSTERESIS

1337

Figure 5. Measured and optimized pressure heads and cumulative infiltrations for (a) experiment 1 and (b)experiment 2 during the first infiltration stage.
Despite differences in the sets of optimized parameters, as
discussed above, the other simultaneously fitted parameters( a w , n, and K s ) showed remarkable consistency among the
optimization runs. The optimized K s value ranged from0.00965 to 0.0119 cm s21 and from 0.00810 to 0.00952 cm s21
for experiments 1 and 2, respectively, while values of a w
ranged from 0.0357 to 0.0527 cm21 and from 0.0484 to 0.0577
cm21, and those of n between 5.42 and 6.77 and between 6.36 and
7.50, for experiments 1 and 2, respectively. On the basis of anumerical ranking of the value of the objective function andconsidering the values of optimized parameters, the best solution
seems to be run 5 with us and l as additional fitting parameters.
The parameters optimized for the infiltration process werenext used to predict observed pressure heads for the redistribution phase of the experiment. Hysteresis was either neglected ( a d 5 a w ) or considered ( a d 5 a w / 2) using thesimplifying assumption of (17). A comparison of calculated
and measured pressure heads at lower tensiometer for experiment 1 is shown in Figure 6a. This figure clearly shows that
although the infiltration part of the experiment was fitted very
successfully, the resulting optimized parameters could not be
used to describe the redistribution stage. Measured pressure

heads decreased much faster and dropped to lower values thanthose simulated when hysteresis was neglected. On the other
hand, considering hysteresis with the simplifying assumption of(17) resulted in calculated pressure heads that decreased fasterthan the measured values. Measured pressure heads were always in between those predicted with and without hysteresis.
Similar results were also obtained for experiment 2 (Figure6b). Values for the objective function F for optimizations using
infiltration data only but which applied to the entire experiment, were very high for all cases (between 0.2 and 0.9), as
could be expected from inspection of Figure 6.
4.2. Simultaneous Analysis of Infiltration andRedistribution Data
The soil hydraulic parameters were next estimated using
data from both the infiltration and redistribution phases forcases in which (1) hysteresis was neglected, (2) hysteresis wasconsidered with the parameters a d and a w being fitted independently, and (3) hysteresis was considered but restricted tothe simplifying assumption (17). Two different runs were executed for each case with hysteresis: one set assuming that thewater content corresponding to the initial pressure heads could be

1338

ŠIMŮNEK ET AL.: ESTIMATING HYSTERESIS

Figure 6. Measured and optimized pressure heads at the lower tensiometer for (a) experiment 1 and (b)experiment 2 during both infiltration and redistribution. The input parameters were optimized using infiltration data only. The redistribution stage was simulated with ( a w 5 2 a d ) or without ( a w 5 a d ) hysteresis.
taken from the main wetting branch and a second set assumingthat the initial condition was given by the main drying branch.
Results for cases when hysteresis was not considered aresummarized in Tables 3 and 4, and the measured and simulated flow responses are shown in Figure 7. A relatively good
agreement between measured and calculated cumulative infiltration volumes was obtained for both experiments. The overall fit of the measured pressure head was better than thatshown in Figure 6, where the redistribution phase was predicted using wetting branch hydraulic data. Calculated pressure heads during the final stages of wetting were generally
lower, and those during redistribution generally higher, thanthe measured values, while the wetting front reached the tensiometers earlier than observed. Although values of the objective function were much lower than those for simulations presented in Figure 6, this was often achieved with physicallyunreasonable parameter values (see, for example, u s and l forruns 4* and 5*). Although still relatively small, the optimized
K s now also shows a broader range of values (0.00764 – 0.013
cm s21 and 0.00741– 0.00976 cm s21 for experiments 1 and 2,
respectively) as compared to those derived from infiltration
only. The values of a w (0.0265– 0.0349 cm21 and 0.0208 –0.0410 cm21) and n (2.34 –5.41 and 2.78 –5.00) estimated fromthe combined infiltration and redistribution process are now
considerably lower than those obtained from the infiltration
data only. The lower n values reflect the inability of the numerical model to predict the sharp moisture fronts at bothtensiometers. The need to simulate a sharp wetting front consequently resulted in extremely high values of l for runs 4* and 5*.
Results of the optimizations using the hysteretic hydraulicproperty model with initial conditions defined by the wetting
branch of the retention curve are presented in Figure 8 andTables 5 and 6. As before, only values of the optimized parameters are shown in the tables. All optimizations allowing
for an independent fit of a d and a w (runs 6, 8, 10, and 12)
resulted in excellent agreement between measured and calculated pressure heads and cumulative infiltration volumes forboth experiments. Both the arrival and the shape of the wetting
front, as well as the decline of pressure heads during redistribution, were now captured with remarkable accuracy. Small
discrepancies between the measured and optimized pressure
heads were present during the second stage of infiltration after
the cone was disconnected from the source bottle and only
water remaining in the cone body continued to infiltrate
through the screen (before the redistribution stage started).
This second infiltration stage, lasting between 35 and 40 s, may

ŠIMŮNEK ET AL.: ESTIMATING HYSTERESIS

1339

Figure 7. Measured and optimized pressure heads and cumulative infiltration for experiment 1 (top) andexperiment 2 (bottom) during both infiltration and redistribution. Hysteresis was not considered.

not have been properly characterized in the numerical solution. Some additional water, perhaps from the water supply
tube, may not have been accurately accounted for by boundary
condition (8). When optimized, u s values again were overpredicted, perhaps suggesting some anisotropy of the soil or a
sharper decrease in the hydraulic conductivity function withdecreasing pressure head as reflected by a higher value of l.
Even with these uncertainties related to u s , k A , and l, theother soil hydraulic parameters ( a , n, K s ) remained exceptionally consistent among the various optimization runs. For
example, K s was predicted to be in the range of 0.00971–
0.0120 cm s21 for experiment 1 and 0.0834 – 0.00969 for experiment 2. The parameters a d and a w for runs 6, 8, and 12
ranged from 0.0325 to 0.0328 cm21 and from 0.0499 to 0.0515
for experiment 1 and from 0.0359 to 0.0368 and from 0.0500 to
0.0559 for experiment 2, and slightly lower for run 10. The
parameter n was predicted to be within the narrow intervals of
5.27–5.61 and 4.75–5.59 for the two experiments. All of theseintervals are remarkably narrow considering the relatively highsensitivity of the parameter estimation technique.
The numerical solutions obtained with optimizations constrained by (17) ( a w 5 2 a d ) (runs 7, 9, 11, and 13) generally
underpredicted the cumulative infiltration volumes as well as
the pressure heads during infiltration. However, the redistribution phase was described equally well as when a d and a w
were both optimized. The final values of the objective functions

when both ad and aw were optimized were, on average, about halfof those obtained when constraint (17) was invoked. Although aw
and Ks were again found to be within relatively narrow intervals(0.0583–0.0674 cm21 and 0.00923–0.0123 cm s21 for experiment
1 and 0.0674–0.0753 and 0.00838–0.0102 for experiment 2), n
became then unreasonably high, as it often reached the upper
constraint of 7.50 imposed to avoid numerical instabilities.
Results of the optimizations using the hysteretic model andthe initial condition defined by the drying branch of the retention curve are presented in Figure 9 and Tables 7 and 8. These
inverse solutions now returned higher initial water contents
than was the case with the initial condition given by the wetting
branch. Since all other water flow attributes were similar, thehigher initial water contents resulted in higher optimized u svalues as well. However, when u s was fixed at 0.350, the numerical inversions yielded (as expected) higher values of theanisotropy coefficient k A and exponent l. Values of the otheroptimized parameters ( a w , a d , n, and K s ) remained almost
the same as those obtained when the wetting branch was used
for the initial condition. Again, we obtained excellent correspondence between measured and calculated values of thepressure head and the cumulative infiltration volume when a w
and a d were fitted independently. When a w and a d wereconstrained by (17), the predicted pressure heads during theinfiltration phase were again underpredicted, except for run 11of experiment 1.

1340

ŠIMŮNEK ET AL.: ESTIMATING HYSTERESIS

Figure 8. Measured and optimized pressure heads and cumulative infiltration for experiment 1 (top) andexperiment 2 (bottom) during both infiltration and redistribution. Hysteresis was considered; the initialcondition was assumed to be on a main wetting branch of the retention curve.

4.3. Comparison of Optimized and Independently
Measured Soil Hydraulic Properties
Figure 10 compares the estimated and independently measured soil water retention functions, as well as estimated hydraulic conductivity functions, for experiment 1 (optimizations
from infiltration data only). Retention curves (Figure 10a)obtained by parameter estimation resembled fairly closely theexperimental laboratory wetting curves determined with thecapillary rise (CR) and computer-automated (CA) methods.
The main deviations involve sharp decreases in estimatedwater contents away from saturation at low suction values.
Desaturation occurred at somehow higher tensions than
expected from the laboratory test results. Notice that theinverse solutions for different combinations of optimized
parameters all yielded quite similar retention curves, except
near saturation. Differences near to saturation are primarily
due to uncertainty in the estimated us values (see also Table
3).
Differences between the various optimized hydraulic conductivity functions (Figure 10b) obtained in this work are also

Table 5. Results of Numerical Inversions for Experiment 1 Considering Hysteresis With Initial Conditions Assumed From
the Wetting CurveRunF

u s [–]

a w , cm21n [–]K s , cm s21l [–]k A [–]

a d , cm21

Rank

6
7
8
9
10
11
12
13

0.01250.03310.01220.0291
0.02250.08630.01250.0169

0.453
0.498
zzz
zzz
zzz
zzz
0.453
0.663

0.04990.06400.05150.06520.03890.05830.0501
0.0674

5.35
7.50a
5.27
7.50a
5.61
7.50a
5.41
7.50a

0.01200.01230.00971
0.009230.01180.01220.01200.0119

zzz
zzz
zzz
zzz
6.25
5.01
0.480
0.710

zzz
zzz
1.37
1.52
zzz
zzz
zzz
zzz

0.0325zzz
0.0328zzz
0.0287zzz
0.0325
zzz

2–3
7
1
6
5
8
2–3
4

a

Upper parameter constraint.

ŠIMŮNEK ET AL.: ESTIMATING HYSTERESIS

1341

Table 6. Results of Numerical Inversions for Experiment 2 Considering Hysteresis With Initial Conditions Assumed From
the Wetting CurveRunF

u s [–]

a w , cm21n [–]K s , cm s21l [–]k A [–]

a d , cm21

Rank

6
7
8
9
10
11
12
13

0.01280.03840.01240.03710.01300.06540.01220.0229

0.423
0.443
zzz
zzz
zzz
zzz
0.389
0.615

0.05530.07120.05590.07220.04380.06740.05000.0753

5.59
7.45
5.31
7.50a
4.75
6.72
4.94
7.50a

0.00966
0.01010.00834
0.00838
0.009670.01020.00969
0.00981

zzz
zzz
zzz
zzz
3.58
2.08
1.35
20.577

zzz
zzz
1.25
1.33
zzz
zzz
zzz
zzz

0.0364zzz
0.0368zzz
0.0329zzz
0.0359
zzz

3
7
2
6
4
8
1
5

a

Upper parameter constraint.

relatively small and would be acceptable for most practicalapplications. Hydraulic conductivities at or near saturationvary only by about 20%. The largest deviations occurred for the
optimization run considering anisotropy when the saturated
hydraulic conductivity in the vertical direction, K vs , was smaller
than those obtained with the other runs. On the other hand,the estimated saturated hydraulic conductivity in the radial
direction, K rs , was somewhat higher than those determinedwith other methods. The optimized K s values were on average
only about 15% and 40% smaller for experiments 1 and 2,
respectively, than K s values determined from a slug test per-formed with the cone penetrometer at the same location (Table 1). Numerically estimated values of K s were approximately
50% and 25% higher than those obtained with the slug test,
respectively. Finally, inverse solutions of the cone penetrometer test gave hydraulic conductivities which were approximately
three times larger than those determined with the constanthead laboratory and Guelph permeameter tests.
Figures 11 and 12 present for experiments 1 and 2, respectively, the estimated and measured soil water retention andestimated hydraulic conductivity functions obtained using datafrom the entire cone penetrometer test (infiltration plus redis-

Figure 9. Measured and optimized pressure heads and cumulative infiltrations for experiment 1 (top) andexperiment 2 (bottom) during both infiltration and redistribution. Hysteresis was considered; the initialcondition was assumed to be on a main draining branch of the retention curve.

1342

ŠIMŮNEK ET AL.: ESTIMATING HYSTERESISTable 7. Results of Numerical Inversions for Experiment 1 Considering Hysteresis With Initial Conditions Assumed From
the Drying CurveRunF

u s [–]

a w , cm21n [–]K s , cm s21l [–]k A [–]

a d , cm21

Rank

6
7
8
9
10
11
12
13

0.01290.02990.01250.02570.03730.10330.01270.0278

0.522
0.580
zzz
zzz
zzz
zzz
0.520
0.568

0.05010.05930.05200.06070.0406
0.06150.04690.0623

5.37
5.36
5.01
4.85
7.50a
6.74
4.91
5.39

0.01200.01240.008630.008200.01160.01120.01200.0122

zzz
zzz
zzz
zzz
7.32
10.4
0.982
0.0561

zzz
zzz
1.63
1.83
zzz
zzz
zzz
zzz

0.0315zzz
0.0319zzz
0.0290zzz
0.0306
zzz

3
6
1
4
7
8
2
5

a

Upper parameter constraint.

tribution), including hysteresis, and with the initial condition
defined by the drying branch of the retention curve. Only the
curves obtained with an independent fit for a w and a d are
given. The retention curves (Figures 11a and 12a), again, show
some uncertainty with respect to u s when this parameter was
optimized. The correspondence between wetting branch retention curves for the different inverse solutions is now slightly
better than that shown in Figure 10. The wetting branches ofthe optimized curves are reasonably close to the independently
measured curves for the wetting branch. Although the drying
retention curves are located higher (with respect to the pressure head axis) than the wetting curves, their shapes are similar. This is due to the assumption that the parameter n is thesame for both wetting and drying [Kool and Parker, 1987].
Given the excellent fit of the pressure head data, and in view of
the uncertainty in the description of the second infiltration
phase of the experiment, we did not expect that the available
data would allow for an independent fit of n for the wetting
and drying branches. Optimized curves were situated between
the laboratory wetting and drying branches of the retention
curve for most of the pressure head range. The restriction
imposed on n probably caused a relatively poor description of
the drying branch of the retention curve as compared to results
obtained with the pressure plate method. The pressure plate
method resulted in n values that differed significantly fromthose obtained with the capillary rise or computer-automatedmethods. Since the redistribution phase covered only a relatively limited part of the drying branch, with no information onwater contents and water fluxes available, n should be morereflective of the wetting than for the drying branch if indeed
Kool and Parker’s [1987] model is appropriate for this particular soil.
The hydraulic conductivity K(h) functions (Figures 11b and12b) exhibit distinct hysteretic loops. Except for run 8, theoptimized K(h) functions again show better correspondence

between particular branches (wetting or drying) for the different optimizations than those shown in Figure 10b. Optimization run 8 was obtained with the assumption of soil anisotropy.
While the vertical saturated conductivity for this curve wasestimated to be about 20% lower than for the other parameter
combinations (as shown in Figures 11 and 12), the estimatedhorizontal conductivity was about 10% higher.
In our experiments we obtained, similar to Kodešová et al.
[1998b], somewhat larger estimates of u s and K s for experiment 1 having a higher applied pressure head as compared toexperiment 2. Still, we obtained good agreement between thehydraulic properties obtained from the first and second cone
penetrometer experiments. This agreement is especially goodfor the estimated properties when hysteresis was considered.

5.

Summary and Conclusions

We have presented a parameter estimation analysis of cone
permeameter data for simultaneous determination of the soil
water retention and hydraulic conductivity functions characterizing the hydraulic properties of unsaturated soils. We
showed that one can estimate not only the wetting branches ofthe soil hydraulic characteristics from the infiltration part ofthe experiment but also simultaneously the wetting and drying
curves via analysis of both the infiltration and redistribution
parts of the test. Hysteresis was described using a relatively
simple empirical model, which scales the wetting scanning
curves from the main wetting curve and the drying scanning
curves from the main drying curve. Drying and wetting curves
were characterized with the same parameters, n, u r , u s , andK s , and with different parameters a d and a w which we either
assumed to be independent, or constrained by the relationshipa w 5 2 a d . Fixing u r and u s at reasonable values and optimizing the other parameters resulted in a less satisfactory fit of themeasured flow response. To improve the solution, u s or addi-

Table 8. Results of Numerical Inversions for Experiment 2 Considering Hysteresis With Initial Conditions on a Drying
Curve
Run

F

u s [–]

a w , cm21n [–]K s , cm s21l [–]k A [–]

a d , cm21

Rank

6
7
8
9
10
11
12
13

0.01298
0.035780.012430.032250.017640.07158
0.01194
0.02195

0.450
0.544
zzz
zzz
zzz
zzz
0.453
0.592

0.0551
0.0647
0.05690.06850.04450.06750.04020.0727

5.69
4.22
5.65
5.42
5.99
7.50a
3.97
7.19

0.009640.01040.007900.007930.009470.01010.009800.00980

zzz
zzz
zzz
zzz
4.25
2.49
4.35
20.470

zzz
zzz
1.36
1.47
zzz
zzz
zzz
zzz

0.0359zzz
0.0363zzz
0.0333zzz
0.0306
zzz

3
7
2
6
4
8
1
5

a

Upper parameter constraint.

ŠIMŮNEK ET AL.: ESTIMATING HYSTERESIS

1343

tional parameters such as the anisotropy coefficient, k A , andthe pore-connectivity parameter, l, were optimized. Interestingly, soil anisotropy had the same effect (in terms of themeasured or predicted flow response) as increasing u s or l. It
seems that our results indicate that knowledge of the initialand final water contents of the soil around the screen is necessary to clearly distinguish between the effects of these threeparameters on the resulting water flow field. Water content
data were not measured in the experimental setup as presented
here (or by Gribb [1996] and Gribb et al. [1998]); incorporationof time domain reflectometry (TDR) elements in the cone
permeameter could provide the required additional information.
Although it was not possible from the data collected with thecurrent measurements (cumulative infiltration volumes andpressure heads at two positions above the water source) todistinguish between the specific effects of u s , l, and k A on thewater flow field, the parameters a w , a d , n, and K s were all

Figure 11. (a) Measured and estimated soil water retentionand (b) estimated hydraulic conductivity curves using combined infiltration and redistribution data from experiment 1
and considering hysteresis (CR, capillary rise (wetting curve);
CA, computer-automated test (wetting curve); PP, pressure
plate (drying curve); d, drying; w, wetting).

Figure 10. (a) Measured and estimated soil water retentionand (b) estimated hydraulic conductivity curves using infiltration data of experiment 1 (CR, capillary rise (wetting curve);
CA, computer-automated test (wetting curve)).

uniquely estimated and showed a high degree of consistency.
We found that the simple hysteretic model was satisfactory fordescribing differences in wetting and drying behavior observed
during this cone experiment. Simultaneous analysis of both
infiltration and redistribution, while considering hysteresis in
the soil hydraulic properties, increased the identifiability of the
optimized parameters with respect to the solution based solely
on wetting process data. The estimated soil hydraulic properties estimated from cone permeameter experiments 1 and 2
were similar to each other. The estimated wetting branches of
the retention curve compared fairly well with retention curves
fitted to results obtained with the capillary rise and computerautomated methods. Failure to similarly describe closely thedrying branch obtained using the pressure plate method was
attributed to limitations in the implemented hysteresis model
which assumes the same n values for both the wetting anddrying branches of the retention curve. The estimated soil
saturated hydraulic conductivity compared well with independent estimates.

1344

ŠIMŮNEK ET AL.: ESTIMATING HYSTERESIS

Figure 12. (a) Measured and estimated soil water retentionand (b) estimated hydraulic conductivity curves using combined infiltration and redistribution data from experiment 2
and considering hysteresis (CR, capillary rise (wetting curve);
CA, computer-automated test (wetting curve); PP, pressure
plate (drying curve); d, drying; w, wetting).

In this study we analyzed data collected after the cone permeameter was carefully buried in the soil profile to limit theeffects of soil disturbance on the optimized soil hydraulic properties. Although Kodešová et al. [1998b] investigated the effectsof cone placement on optimized soil hydraulic properties andconcluded that disturbance had a little impact on the values ofthe estimated parameters for this particular sandy soil, we
believe that additional studies of placement and soil disturbance must be carried out for different soils.
Acknowledgment. Thanks are due to reviewers (Jos van Dam andWolfgang Dürner) for their many useful suggestions and comments.

References
Aziz, K., and A. Settari, Petroleum Reservoir Simulation, pp. 395– 401,
Applied Science, Barking, U. K., 1979.
Bard, Y., Nonlinear Parameter Estimation, 341 pp., Academic, San
Diego, Calif., 1974.

Bower, H., and R. C. Rice, A slug test for determining hydraulicconductivity of unconfined aquifers with completely or partially penetrating wells, Water Resour. Res., 12(3), 432– 428, 1976.
Ciollaro, G., and N. Romano, Spatial variability of the soil hydraulicproperties of a volcanic soil, Geoderma, 65, 263–282, 1995.
Clausnitzer, V., and J. W. Hopmans, Non-linear parameter estimation:
LM_OPT, general-purpose optimization code based on the Levenberg-Marquardt algorithm, Land, Air and Water Resour. Pap.
100032, Univ. of Calif., Davis, 1995.
Dane, J. H., and S. Hruska, In-situ determination of soil hydraulicproperties during drainage, Soil Sci. Soc. Am. J., 47(3), 619 – 624,
1983.
Eching, S. O., and J. W. Hopmans, Optimization of hydraulic functions
from transient outflow and soil water pressure data, Soil Sci. Soc.
Am. J., 57, 1167–1175, 1993.
Gillham, R. W., A. Klute, and D. F. Heermann, Measurement andnumerical simulation of hysteretic flow in a heterogeneous porous
medium, Soil Sci. Soc. Am. J., 43(6), 1061–1067, 1979.
Gribb, M. M., Parameter estimation for determining hydraulic properties of a fine sand from transient flow measurements, Water Resour. Res., 32(7), 1965–1974, 1996.
Gribb, M. M., J. Šimůnek, and M. F. Leonard, Development of a cone
penetrometer method to determine soil hydraulic properties, J.
Geotech. Geoenviron. Eng., 124(9), 820 – 829, 1998.
Hoa, N. T., R. Gaudu, and C. Thirriot, Influence of the hysteresis
effect on transient flows in saturated-unsaturated porous media,
Water Resour. Res., 13, 992–996, 1977.
Hvorslev, M. J., Time lag and soil permeability in ground-water observations, WES Bull. 36, U.S. Army Corps of Eng., Vicksburg,
Miss., 1951.
Inoue, M., J. Šimůnek, J. W. Hopmans, and V. Clausnitzer, In-situ
estimation of soil hydraulic functions using a multistep soil water
extraction technique, Water Resour. Res., 34(5), 1035–1050, 1998.
Jaynes, D. B., Estimating hysteresis in the soil water retention function,
in Proceedings of the International Workshop on Indirect Methods for
Estimating the Hydraulic Properties of Unsaturated Soils, edited by
M. T. van Genuchten et al., pp. 219 –232, Univ. of Calif., Riverside,
1992.
Kaluarachchi, J. J., and J. C. Parker, Effects of hysteresis with air
entrapment on water flow in the unsaturated zone, Water Resour.
Res., 23(10), 1967–1976, 1987.
Klute, A., and D. F. Heerman, Soil water profile development under aperiodic boundary condition, Soil Sci., 117, 265–271, 1974.
Kodešová, R., M. M. Gribb, and J. Šimůnek, A new CPT method for
estimating soil hydraulic properties, in Proceedings of the First International Conference on Site Characterization, vol. 2, edited by P. K.
Robertson and P. W. Mayne, pp. 1421–1425, A. A. Balkema, Rotterdam, Netherlands, 1998a.
Kodešová, R., M. M. Gribb, and J. Šimůnek, Estimating soil hydraulicproperties from transient cone permeameter data, Soil Sci., 163(6),436 – 453, 1998b.
Kool, J. B., and J. C. Parker, Development and evaluation of closedform expressions for hysteretic soil hydraulic properties, Water Resour. Res., 23(1), 105–114, 1987.
Kool, J. B., and J. C. Parker, Analysis of the inverse problem fortransient unsaturated flow, Water Resour. Res., 24(6), 817– 830, 1988.
Kool, J. B., J. C. Parker, and M. T. van Genuchten, Determining soil
hydraulic properties from one-step outflow experiments by parameter estimation, I, Theory and numerical studies, Soil Sci. Soc.
Am. J., 49, 1348 –1354, 1985.
Kool, J. B., J. C. Parker, and M. T. van Genuchten, Parameter estimation for unsaturated flow and transport models—A review, J.
Hydrol., 91, 255–293, 1987.
Lambe, W. T., Capillary phenomena in cohesionless soils, Trans. Am.
Soc. Civ. Eng., 116, 401– 423, 1951.
Luckner, L., M. T. van Genuchten, and D. R. Nielsen, A consistent setof parametric models for the two-phase flow of immiscible fluids inthe subsurface, Water Resour. Res., 25(10), 2187–2193, 1989.
Marquardt, D. W., An algorithm for least-squares estimation of nonlinear parameters, SIAM J. Appl. Math., 11, 431– 441, 1963.
Mualem, Y., Modified approach to capillary hysteresis based on asimilarity hypothesis, Water Resour. Res., 9(5), 1324 –1331, 1973.
Mualem, Y., A conceptual model of hysteresis, Water Resour. Res.,
10(3), 514 –520, 1974.
Mualem, Y., A new model for predicting the hydraulic conductivity ofunsaturated porous media, Water Resour. Res., 12(3), 513–522, 1976.

ŠIMŮNEK ET AL.: ESTIMATING HYSTERESIS
Mualem, Y., A modified dependent-domain theory of hysteresis, Soil
Sci., 137, 283–294, 1984a.
Mualem, Y., Anisotropy of unsaturated soils, Soil Sci. Soc. Am. J., 48,
505–509, 1984b.
Mualem, Y., and G. Dagan, A dependent domain model of capillary
hysteresis, Water Resour. Res., 11(3), 452– 460, 1975.
Mualem, Y., and H. J. Morel-Seytoux, Analysis of a capillary hysteresis
model based one a one-variable distribution function, Water Resour.
Res., 14(4), 605– 610, 1978.
Nielsen, D. R., and L. M. Luckner, Theoretical aspects to estimate
reasonable initial parameters and range limits in indentification procedures for soil hydraulic properties, in Proceedings of the International Workshop on Indirect Methods for Estimating the Hydraulic
Properties of Unsaturated Soils, edited by M. T. van Genuchten et al.,pp. 147–160, Univ. of Calif., Riverside, 1992.
Ray, R. P., and K. B. Morris, Automated laboratory testing for soil/
water characteristic curves, in Proceedings of 1st International Conference on Unsaturated Soils, vol. 2, edited by E. Alonso and P.
Delage, pp. 236 –241, A. A. Balkema, Rotterdam, Netherlands,
1994.
Reynolds, W. D., Unsaturated hydraulic conductivity: Field measurement, in Soil Sampling and Methods of Analysis, edited by M. R.
Carter, pp. 633– 644, Lewis, New York, 1993.
Royer, J. M., and G. Vachaud, Field determination of hysteresis insoil-water characteristics, Soil Sci. Soc. Am. J., 39(2), 221–223, 1975.
Russo, D., W. A. Jury, and G. L. Butters, Numerical analysis of solutetransport during transient irrigation, 1, The effect of hysteresis andprofile heterogeneity, Water Resour. Res., 25(10), 2109 –2118, 1989.
Santini, A., N. Romano, G. Ciollaro, and V. Comegna, Evaluation ofa laboratory inverse method for determining unsaturated hydraulicproperties of a soil under different tillage practices, Soil Sci., 160,
340 –351, 1995.
Scaturo, D. M., Evaluation of multi-level direct push sampling forhydraulic conductivity analysis, M. S. thesis, Univ. of S. C., Columbia, 1993.
Scott, P. S., G. J. Farquhar, and N. Kouwen, Hysteresis effects on net
infiltration, in Advances in Infiltration, Publ. 11-83, pp. 163–170, Am.
Soc. of Agric. Eng., St. Joseph, Mich., 1983.
Šimůnek, J., and M. T. van Genuchten, Parameter estimation of soil
hydraulic properties from the tension disc infiltrometer experimentby numerical inversion, Water Resour. Res., 32(9), 2683–2696, 1996.
Šimůnek, J., and M. T. van Genuchten, Estimating unsaturated soil
hydraulic properties from multiple tension disc infiltrometer data,
Soil Sci., 162(6), 383–398, 1997.
Šimůnek, J., M. Šejna, M., and M. T. van Genuchten, The HYDRUS-2D software package for simulating water flow and solutetransport in two-dimensional variably saturated media, version 1.0,
IGWMC - TPS 53, Int. Ground Water Modeling Cent., Colo. Schoolof Mines, Golden, 1996.
Šimůnek, J., O. Wendroth, and M. T. van Genuchten, A parameter

1345

estimation analysis of the evaporation method for determining soil
hydraulic properties, Soil Sci. Soc. Am. J., 62(4), 894 –905, 1998.
Singleton, J. E., Hydraulic characteristics of a laboratory aquifer, M. S.
thesis, Univ. of S. C., Columbia, 1997.
Topp, G. C., Soil-water hysteresis measured in a sandy loam andcompared with the hysteretic domain model, Soil Sci. Soc. Am. Proc.,
33, 645– 651, 1969.
Vachaud, G., and J.-L. Thony, Hysteresis during infiltration and redistribution in a soil column at different initial water contents, Water
Resour. Res., 7(1), 111–127, 1971.
van Dam, J. C., J. N. M. Stricker, and P. Droogers, Inverse method for
determining soil hydraulic functions from one-step outflow experiment, Soil Sci. Soc. Am. Proc., 56, 1042–1050, 1992.
van Dam, J. C., J. N. M. Stricker, and P. Droogers, Inverse method todetermine soil hydraulic functions from multistep outflow experiment, Soil Sci. Soc. Am. Proc., 58, 647– 652, 1994.
van Genuchten, M. T., A closed-form equation for predicting thehydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44,892– 898, 1980.
van Genuchten, M. T., Non-equilibrium transport parameters frommiscible displacement experiments, Res. Rep. 119, U.S. Salinity Lab.,
U.S. Dep. of Agric., Riverside, Calif., 1981.
van Genuchten, M. T., F. J. Leij, and S. R. Yates, The RETC code forquantifying the hydraulic functions of unsaturated soils, EPA/600/291-065, US EPA, Off. of Res. and Dev., U.S. Environ. Prot. Agency,
Washington, D. C., 1991.
Vogel, T., K. Huang, R. Zhang, and M. T. van Genuchten, The HYDRUS code for simulating one-dimensional water flow, solute transport, and heat movement in variably-saturated media, version 5.0,
Res. Rep. No 140, U.S. Salinity Lab., U.S. Dep. of Agric., Riverside,
Calif., 1996.
Watson, K. K., R. J. Reginato, and R. D. Jackson, Soil water hysteresis
in a field soil, Soil Sci. Soc. Am. Proc., 39, 242–246, 1975.
Yeh, W. W.-G., Review of parameter identification procedures ingroundwater hydrology: The inverse problem, Water Resour. Res.,
22(2), 95–108, 1986.
Znidarčič, D., T. Illangasekare, and M. Manna, Laboratory testing andparameter estimation for two-phase flow problems, in GeotechnicalEngineering Congress, edited by F. G. McClean et al., Geotech. Spec.
Publ. 27(2), pp. 1089 –1099, 1991.
M. M. Gribb and R. Kodešová, Department of Civil and Environmental Engineering, University of South Carolina, Columbia, SC
29208.
J. Šimůnek and M. T. van Genuchten, U.S. Salinity Laboratory,
ARS, USDA, 450 W. Big Springs Rd., Riverside, CA 92507.
(jsimunek@ussl.ars.uda.gov)(Received March 6, 1998; revised December 3, 1998;
accepted December 7, 1998.)

1346