WATER RESOURCES RESEARCH, VOL. 34, NO. 5, PAGES 1035–1050, MAY 1998

In situ estimation of soil hydraulic functions using
a multistep soil-water extraction technique
M. Inoue
Arid Land Research Center, Tottori University, Tottori, JapanJ. Šimunek
USDA Salinity Laboratory, Riverside, California

J. W. Hopmans and V. Clausnitzer
Hydrology Program, University of California, Davis

Abstract. Estimation of the retention and unsaturated hydraulic conductivity functions isessential to effectively provide input for water flow and transport simulation andprediction. A parameter optimization procedure is shown as a promising tool to estimate
inversely these hydraulic function parameters from transient soil matric potential andcumulative soil solution extraction measurements. Sensitivity analyses from synthetic data
generated from forward numerical model simulations showed that optimum tensiometer
locations will depend on soil type. Experiments were carried out in both a laboratory
column (Columbia sandy loam) and in the field (Yolo silt loam). In both cases a series of
vacuum extraction pressures was applied to a ceramic soil solution sampler, andcumulative soil solution extraction volume and matric potentials at various positions near
the extraction device were monitored as the soil solution was extracted. In the laboratory
a zero-flux boundary condition was maintained at the bottom of the column, whereas
matric potential measurements were used in the field to define the lower boundary. In
both the field and laboratory experiments, flow at the upper boundary was zero.
Cumulative extraction volume and matric potential data were included in the objectivefunction to be minimized to estimate the hydraulic function parameters. We determinedthat the optimized solution was sensitive to the contact between the ceramic ring and thesurrounding soil. By also optimizing the hydraulic resistance of the ceramic ring of theextraction device, optimization improved the fit between measured and optimized flow
variables. Comparison of the optimized with the independently measured hydraulicfunctions indicated that the in situ estimation using a multistep extraction procedure canprovide accurate soil hydraulic data.

1.

Introduction

As the concern for a safe, clean environment and highgroundwater quality increases, the importance of an accurate
soil physical description of the combined unsaturatedsaturated porous system is increasingly recognized in the fields
of environmental engineering and groundwater hydrology.
Moreover, accurate soil physical data is required for the suiteof increasingly available agrohydrological simulation models
and in the modeling of land surface processes to simulate theexchange of sensible and latent heat between the soil andatmosphere. With this wider interest the spatial scale of interest has shifted to larger dimensions. Soil hydraulic and transport characterization is needed for soil-water systems as large
as a watershed and for depths extending from the rooting zone
to the groundwater. This trend in larger spatial scales brings
with it the need to consider soil heterogeneity within the considered system. Therefore methodologies need to be develCopyright 1998 by the American Geophysical Union.
Paper number 98WR00295.
0043-1397/98/98WR-00295$09.00

oped that allow for a rapid and accurate characterization of thesoil physical properties and its spatial variability.
Currently, many laboratory and field methods exist to determine the highly nonlinear soil hydraulic functions of the soil
water retention and unsaturated hydraulic conductivity curves.
Most methods require restrictive initial and boundary conditions, which make measurements time consuming, range restrictive, and expensive. Excellent reviews have been publishedby Klute and Dirksen [1986], Klute [1986], and Green et al.
[1986]. Dirksen [1991] discusses the application of parameter
estimation by the inverse solution technique as a viable alternative to other traditional laboratory methods for soil hydraulic characterization. The speed at which results for a large
number of soil samples can be obtained via multiplexing of
sensors, and recent developments in improved TDR and pressure sensors make this inverse technique a viable alternative.
The parameter optimization method is especially advantageous since single experiments can yield simultaneous estimates of the soil water retention and unsaturated hydraulicconductivity functions.
Since the early 1970s the inverse problem of parameter
identification for distributed numerical models has been applied in groundwater hydrology and field petroleum engineer-

1035

1036

INOUE ET AL.: IN SITU ESTIMATION

ing. Yeh [1986] presented an excellent review of the inverse
problem as a parameter identification procedure in groundwater hydrology. Its application to the vadose zone started laterand has been limited mostly to parameter estimation of soil
hydraulic properties. Parameter estimation, as defined in thisstudy, involves the indirect estimation of soil hydraulic functions by numerical solution of the governing flow equation andsubsequent comparison of the numerical solution with experimental data. In this procedure, soil hydraulic properties are
assumed to be described by an analytical model, with as yet
unknown parameter values. An experiment is then set up under controlled conditions with prescribed initial and boundary
conditions. During the experiment, one or more flowcontrolled variables are measured. Subsequently, the Richards
equation is solved numerically using the hydraulic functions
with initial estimates for their parameters. These parameters
are then optimized by minimization of an objective function
containing the sums of squared deviations between observed
and predicted flow variables, using repeated numerical simulations of the flow process. This iterative inversion of the flow
equation is in contrast to direct inversion techniques as used inanalytical solutions.
Among the first to suggest the application of computer models to estimate soil hydraulic parameters were Whisler andWatson [1968], who reported on estimation of the saturated
conductivity of a draining soil by matching observed and simulated drainage. The pressure plate outflow method was introduced by Gardner [1956], in which an initially saturated soil was
subjected to a series of step increases in air pressure with thedrainage or outflow measured after each pressure step increase. Given specific assumptions, the analytical solution
yields the soil diffusivity as a function of water content. Doering[1965] improved the outflow method by proposing a one-step
experiment, thereby achieving considerable time savings without loss in accuracy. Other modifications were introduced
more recently by Valiantzas and Kerkides [1990], who extended
the outflow method to simultaneous determination of soil water retention and unsaturated hydraulic conductivity using the
Brooks and Corey [1966] formulation of the soil hydraulic properties.
Although the application of the inverse approach to theoutflow method appeared promising, problems were encountered with the nonuniqueness of the optimized parameters[van Dam, 1990]. Nonuniqueness leads to more than one set ofparameters, yielding minimum values for the objective function determined by local minima, or by the same global minimum at more than one point in the parameter space [Carrera
and Neuman, 1986b]. The study of nonuniqueness problems
has lead to many investigations on the type of experiments, andthe measured flow variables that need to be included in theobjective function.
Kool et al. [1985] and Parker et al. [1985] were the first to
apply the inverse approach by numerical solution of the Richards equation for a transient one-step outflow experiment.
They concluded that uniqueness problems were minimized if
the experiment were designed to cover a wide range in watercontents. Kool and Parker [1988] discussed the advantage ofincluding tensiometric data simultaneously with the outflow
measurements in the inverse approach from a hypotheticalinfiltration and redistribution experiment. In addition, theanalysis of the objective function by Toorman et al. [1992]
indicated that uniqueness problems were minimized if soil
water pressure head data were included in the objective func-tion of a transient one-step outflow experiment. To circumventthe need for additional soil water pressure measurements inthe outflow experiment, van Dam et al. [1994] conducted outflow experiments in which the pneumatic pressure was increased in several smaller steps. Their work for a loam soil
showed that the outflow data of a multistep experiment containsufficient information for unique estimates of the soil hydraulicfunctions. The experimental work by Eching and Hopmans[1993] and Eching et al. [1994] showed how the multistep
method, in combination with automated soil water pressure
measurements during drainage of the soil core, resulted inunique parameter values for the optimized soil hydraulic functions for four different textured soils.
The most recent applications of the parameter optimization
approach for the estimation of soil hydraulic properties include
two-dimensional hypothetical experiments and flow simulations. For example, Šimunek and van Genuchten [1996] demonstrated that tension disc permeameter experiments complemented with soil matric potential data guarantee numericalconvergence and uniqueness of the optimized parameters. In
their subsequent study, Šimunek and van Genuchten [1997]
showed that additional soil matric potential data are notneeded, provided the infiltration experiment is carried out withseveral consecutive tensions and the initial and final watercontents below the tension disc are known. Gribb [1996] usedhypothetical infiltration data from a modified conepenetrometer to indirectly estimate the soil hydraulic functions. The
cone-shaped device included a porous filter near the cone
through which water was injected in the unsaturated soil, andtwo tensiometer rings above the filter. As the volume of infiltrated water is monitored, the soil water matric potential response was measured by the two tensiometer rings.
While a majority of applications of the parameter estimation
technique was only for laboratory studies, the objective of thispaper is mainly to demonstrate the potential application of soil
water extraction method to estimate soil water retention andunsaturated hydraulic conductivity parameters in the field. The
feasibility of the vacuum extraction technique is first demonstrated from numerically generated data. Subsequently, theparameter estimation technique was applied to laboratory andin situ field data.

2.

Materials and Methods

2.1.

Water Flow Theory

The proposed method is based on the premise that the soil’s
hydraulic properties can be estimated from the measurementof extracted soil water volume and soil water potential values
at various locations as a function of time. This is done by
applying a number of vacuum increments to a ceramic soil
water extraction device. Although the experiments will occur
spatially in three dimensions, we can assume axial symmetry
for isotropic soils, which reduces the Richards equation to twodimensions, which for a rigid porous media can be written asC­h 1 ­
5
­t
r ­r

S

rK

­h
­r

D

1

­
­z

S D
K

­h
­z

1

­K
­z

(1)

where C is the water capacity [L 21 ], h is the soil matricpotential [L], K is the unsaturated hydraulic conductivity[L T 21 ], r is the radial coordinate [L], z is the verticalcoordinate (positive upward) [L], and t is time [T]. Equation(1) was solved with the Galerkin finite element method based

INOUE ET AL.: IN SITU ESTIMATION

1037

on the mass conservative iterative scheme proposed by Celia etal. [1990]. An example of a typical flow domain with the applied finite element mesh including geometry notation andlocations of matric potential measurements (T 1 through T 3 ) ispresented in Figure 1. The mesh size of the triangular elements
gradually increases from 0.15 cm for the porous cup to 1.3 cm
in the soil. The smaller size in and near the porous cup isrequired since it is there that the largest potential gradients
occur.
Boundary and initial conditions for which (1) was solved are
dependent on the specific experiment, but they can be generally defined as
h~r, z, t! 5 h i~ z!

t 5 t 0, 0 , r , R(2)

z 5 0, 0 , r , R, t 0 , t , t end(3a)z 5 255 cm, 0 , r , R, t 0 , t , t end(3b)r 5 r i~2.70 cm!, 16.5 , z , 19.5 cm

(4)

Lab
q~r, z, t! 5 0
orField
h~r, z, t! 5 h lb
Lab
h~r, z, t! 5 h exField
r 5 r i~2.70 cm!, 211.5 , z , 28.5 cm

h~r, z, t! 5 h ex
q~r, z, t! 5 0

remaining boundaries, t 0 , t , t end(5)where r 5 R denotes the radius of the flow domain [L], theordinate z 5 0 is placed at the bottom or top of the flow
domain for lab and field experiments (see Figure 1), respectively; t 0 and t end correspond to the beginning and end of theextraction experiment [T], respectively; and r i is the insideradius of the extraction device [L]. Equations (3a) and (3b)specify the bottom boundary condition for the lab and field
experimental setups, respectively; with h lb [L] being the measured pressure head at the bottom of considered transportdomain. Equation (4) represents the applied suction h ex [L] inthe extraction sampler. In solving (1), subject to conditions (2)through (5), the unsaturated hydraulic properties are definedby [van Genuchten, 1980; Mualem, 1976]
Se 5

2.2.

u ~h! 2 u r
1
5
us 2 ur~1 1 u a h nu! m

(6)

and
m 2
K~ u ! 5 K sS 1/e 2@1 2 ~1 2 S 1/m
e ! #K~ u ! 5 K s

h$0

Figure 1. Finite element mesh and geometry notation forflow domain of hypothetical and laboratory experiments.

h,0

(7)
(8)

In expressions (6)–(8), S e is the effective saturation [–]; u r andu s denote the residual and saturated volumetric water contents[–], respectively; a [L 21 ] and n [–] (m 5 1 2 1/n) areempirical parameters; and K s is a fitted saturated hydraulicconductivity [L T 21 ], not necessarily equivalent to an independently measured saturated hydraulic conductivity. Since
the ceramic ring is considered part of the porous system in the
numerical simulations, parameter values were chosen such thatthe ceramic remained fully saturated over the imposed rangeof h ex .

Parameter Optimization

Parameters in (6)–(8) were estimated from maximization of
the log likelihood function [Bard, 1974], which includes differences between the observed and predicted flow variables. Assuming measurement errors to be independent with zeromean, the parameter optimization procedure is equivalent tominimization of a weighted, least squares problem, which is
cast in an objective function, OF(b), with b denoting the vector
containing the optimized parameters

OS Om

OF(b) 5

nj

wj

j51

i51

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

D(9)where j represents the different sets of measurements (cumulative extraction volume, matric potential head at different
locations, or water volume in flow domain), n j is the number ofmeasurements within a particular set, q *j (t i ) are measurementsof type j at time t i , q j (t i , b) are the corresponding model

1038

INOUE ET AL.: IN SITU ESTIMATION

optimized parameters has been given by Kool and Parker[1988], Yeh [1986], and Carrera and Neuman [1986a, b].
The sensitivity coefficients (s i, j ) for the hypothetical experiments described below were calculated from [Šimunek andvan Genuchten, 1996]
s i, j 5 100b j

­q i
­b j

(11)

where s i, j denotes the change of measurement variable q irelative to a 1% change of the parameter b j , and the partial
derivative term is estimated from­q i q i~b 1 Dbe j! 2 q i~b!
5
­b j
Db j

(12)

where e j is the jth unit vector, and Db j 5 0.01b j . Equation(11) allows comparison of sensitivities between parameters,
independent of their unit or absolute value.
The time-averaged sensitivity, S(r, z, b j ), and maximum
sensitivity, S max(r, z, b j ), of the measured soil water potential
values to each of the four parameters b j were calculated for thewhole flow domain from
S~r, z, b j! 5

Figure 2. Layout of laboratory experiment.

1~t 0 2 t end!

E

tend

s h, j~r, z, t! dt

(13a)

t0

and
predictions using the parameters in b, and w j and w i, j areweighting factors associated with data type and data point,respectively. Assuming that measurement errors for all pressure transducers were identical, w i, j was set equal to 1 for allpressure measurements, whereas the water volume measurements were given a weighting factor of 10. Differences inweighting between data types, as caused by differences in magnitudes and their number n j , were minimized by division ofeach data point by the variance of the measurements of data
type j and by n j [Clausnitzer and Hopmans, 1995]
wj 5

1 1
n j s 2j(10)

where s j and n j denote the standard deviation and the numberof data of the j-type measurements, respectively. Thus OF(b) isequal to the weighted average squared deviation between simulated and measured flow variables. An effective method tominimize (9), which uses a combination of the Newton andsteepest descent method, was proposed by Marquardt [1963].
Details of this procedure have been given by Kool et al. [1987]
and Šimunek and van Genuchten [1996]. It suffices to state that
the Levenberg-Marquardt method is a standard method in
nonlinear least squares fitting which, in addition to the sum ofsquared residuals of (9), also provides confidence intervals forthe optimized parameters.
2.3.

Sensitivity Coefficients

An experiment must be designed such that direct information is available for the least sensitive parameters, thereby
eliminating them from the optimized parameter set or providing good initial well-constrained estimates. Choice of data type
and their measurement in space and time should be based ona sensitivity analysis as well, so that their sensitivity to the
optimized parameters is maximum [Šimunek and van Genuchten, 1996]. A review on ill-posedness and error analysis of the

S max ~r, z, b j! 5 Maximum $s h, j~r, z, t, b j!% t 0 , t , t end(13b)
where S h, j is the sensitivity of the soil matric potential toparameter b j calculated for each node of the flow domain using(10). This information allows us to evaluate the importance oftensiometer location in multistep extraction experiments.
2.4.

Numerical Experiments

Soil water extraction accompanied by soil water potential
response was simulated using forward simulations using thetwo-dimensional flow model HYDRUS-2D of Šimunek et al.
[1996]. The dimensions of the flow domain were R 5 10.4 cm
and L 5 24.3 cm, where L is the coordinate of the soil surface[L] (Figure 1). Dimensions and instrument locations wereidentical as used in the subsequent laboratory experiment. The
3.0-cm-long ceramic ring of the extraction device (l ) was inserted in the center of the soil column (r 5 0), with its center
located 6.3 cm below the soil surface (d). Three hypotheticaltensiometers were installed (see also in Figure 2) at the following positions: r 5 4.0 cm and z 5 18 cm (T 1 ); r 5 6.0 cm
and z 5 18 cm (T 2 ); and r 5 6.0 cm and z 5 0.3 cm (T 3 ).
Forward hypothetical simulations were done using a sandy
loam and silt soil, for which Carsel and Parrish [1988] presented
the hydraulic parameters. Parameter values are presented inTable 1. The three soils represent a wide range in soil hydraulicproperties, so that simulations for likely soil types could be
evaluated from these data. As the initial condition, hydraulicequilibrium was assumed with a soil water potential of 0 cm at
the bottom of the soil column, corresponding with a soil water
potential of 224.3 cm at the soil surface. The following vacuum
steps were applied to the extraction device: h ex 5 230 cm (0 ,t , 60 h), h ex 5 260 cm (60 , t , 120 h), h ex 5 2120 cm(120 , t , 240 h), h ex 5 2240 cm (240 , t , 360 h), andh ex 5 2480 cm (360 , t , 600 h).

INOUE ET AL.: IN SITU ESTIMATION

Soil matric potential at the three hypothetical locations (h 1 ,
h 2 , and h 3 ; or q j , j 5 1, 2, and 3 in (9)) and cumulative
extraction volume (Q, q 4 in (9)) were simulated for these soil
types for the given initial and boundary conditions. In the next
step the parameters a , n, u r , and K s were optimized fromdiscrete values of the simulated Q and soil matric potential
values h 1 through h 3 by inverse solution using identical initialand boundary conditions. The saturated water content, u s , wasassumed to be known. Four sets of initial parameter values foreach optimization were taken from Carsel and Parrish [1988]
using parameters for the sandy loam, loam, silt, and clay (truevalues for sand); sand, loam, silt loam, and clay (silt); and sand,sandy loam, silt, and clay (loam), which combined include thewhole range of possible parameter values. The correspondenceof each of the four sets of final parameter values with their true
values used in the forward simulations yields information on
the uniqueness of the parameters for the given experimental
conditions and soils.
2.5.

Laboratory Experiment

The flow regime was identical to that used in the numericalexperiments. The extraction device consists of a 6.0-cmoutside-diameter (OD) PVC pipe with a total vertical length of
11.5 cm, with a 3-cm-long porous ceramic ring (Soil Moisture
Equipment Corporation) glued in the tube, 6.3 cm from thetop. The inside (r i ) and outside (r 0 ) radiuses of the ceramic
ring were 2.70 and 3.02 cm, respectively. Tensiometers wereconstructed by cementing a 1-cm-long 0.635-cm-OD ceramiccup (Soil Moisture Equipment Corporation) to 0.63-cm-OD
acrylic tubing. After filling with water, a pressure transducer(Omega Technologies Corporation) was connected to the water-filled tubing using a short piece of tygon tubing. Transducers were connected to a 21X data logger using an AM416
multiplexer (Campbell Scientific Corporation), providing soil
matric potential measurements with an accuracy of about 2 cm.
The K s value for the ceramic ring material, K cer [LT21], wasindependently measured by imposing a pressure difference
across the ring, while submerged in water. The Darcy equationin radial coordinates was used to solve for K cer analytically,leading tor0ri
K cer 52 p l~h out 2 h in!
Q ln

(14)

where Q is the steady state volumetric flow rate [L 3 T 21 ]extracted through the ring of length l [L], and h out and h indenote the soil matric potential values at the outside and insideof the ring [L] at radial distances of r 0 and r i from the centerof the ring [L] (Figure 1), respectively. Replicated measurements yielded K cer 5 8.334 3 1024 cm h21.
A schematic of the laboratory experiment, which includes
the burette assembly for vacuum extraction and the boundary
conditions, is presented in Figure 2. The soil used was a Columbia fine sandy loam. The soil was air dried, sieved througha 2-mm screen, and packed uniformly to a dry bulk density of
1.46 g/cm3 in the soil column with porous plastic on the bottom. The soil was saturated from the bottom by applying aslight positive pressure to the water. After saturation, hydraulicequilibrium was established by applying a water potential of
about 0 cm at the bottom of the column using a Mariotte
device for 6 days. At this time, drainage rate was zero andtensiometer measurements indicated a zero total soil water

1039

Table 1. Soil Hydraulic Parameters for Selected Textural
Classes
Soil Type

us

ur

a,
cm21

n

K s,
cm h21

Sand
Sandy loam
Loam
Silt loam
Silt
Clay

0.43
0.41
0.43
0.45
0.46
0.38

0.045
0.065
0.078
0.067
0.034
0.068

0.145
0.075
0.036
0.02
0.016
0.008

2.68
1.89
1.56
1.41
1.37
1.09

29.7
4.42
1.04
0.45
0.25
0.02

Carsel and Parish [1988].

potential gradient, with a soil matric potential of 23.3 cm at
z 5 0.3 cm (T 3 ) and 221 cm at z 5 18 cm (average of T 1 andT 2 ). Therefore, assuming hydraulic equilibrium, h i is 23.0 cm
at the bottom and linearly decreases to 227.3 cm at the top of
the column. Measurement accuracy of the tensiometer readings was determined from separate calibration of each transducer and varied between 61–2 cm. A 0.7-cm-thick PVC plate
to prevent evaporation of the soil during the extraction experiment covered the column.
The water-filled chamber of the extraction device was connected to a burette by tygon tubing. While applying vacuum tothe air phase in the burette, the extracted soil water flows into
the burette with the pressure of the inflow end equal to theapplied vacuum. The volume of extracted water was measured
both manually and using a pressure transducer in the bottomof the burette (Figure 2). A film of oil to prevent water loss by
evaporation covered the water in the burette. The measurement error of the burette readings was about 0.2 mL. Whilevacuum was applied to the extraction device, the water pressure in the chamber was measured simultaneously using atensiometer connected to a pressure transducer. Average vacuum steps (as measured in the chamber of the extraction device) were h ex 5 235 cm (0 , t , 25 h), h ex 5 265 cm (25 ,t , 73 h), h ex 5 2125 cm (73 , t , 217 h), h ex 5 2240 cm(217 , t , 339 h), and h ex 5 2480 cm (339 , t , 605 h).
From weight measurements of the column at the end of theextraction experiment at 605 h (t end), after oven drying, and at
the initial time (t 0 ), values for column average volumetricwater content and dry soil bulk density were calculated. The
dry bulk density was 1.46 g cm23, whereas the column-averagevolumetric water content values at t 0 and t end were 0.37 and0.22 cm3 cm23, respectively. The volumes of water (cm3) corresponding with these two water content values (2939 and 1758mL) were included in the objective function, thereby insuringthat the initial and final conditions of the experiment andsimulation model were identical.
The soil water retention function was independently determined from a multistep outflow experiment [Eching and Hopmans, 1993] for the same soil with a slightly lower bulk densityof 1.42 g cm23. The soil was air dried and sieved, and a 216-cm3
sample was packed in a 6.4-cm-diameter core. Subsequently,
the soil was saturated with water in a Tempe pressure cell.
Various air pressure increments were applied to drain thesample, during which the soil matric potential in the center ofthe core and total drainage volume was measured.
2.6.

Field Experiment

A detailed overview of the field experiment is presented in
Figure 3. The field soil is a Yolo silt loam with an approximate

1040

INOUE ET AL.: IN SITU ESTIMATION

1995), soil samples representing the 8-cm soil depth were collected and volumetric water content was determined from ovendrying. These volumetric water content data were matchedwith soil water potential values of the T 1 tensiometers, resulting in the following u 2 h (cm) points: (0.39, 294), (0.36,2133), and (0.35, 2149) (Figure 10a). Also these three independently measured soil water retention points were includedin the OF (9). In addition to the tensiometers used to measure
response to soil water extraction, additional tensiometers were
installed at depths of 8 and 20 cm, at distances of more than 20
cm away from the ceramic extraction device. These matricpotential measurements were not influenced by water extraction, and together with the 40- and 55-cm-deep tensiometerscharacterized the draining of the soil profile by gravity forces.
Therefore these matric potential measurements in combination with water content measurements using a neutron probeat depths of 25, 40, 55, and 70 cm were used to obtain additional independent estimates of in situ u (h) and K( u ) datausing the instantaneous profile method [Green et al., 1986].

Figure 3. Experimental setup of in situ multistep extraction
experiment.

clay content of 22% [Eching and Hopmans, 1993]. Soil was
excavated to a depth of 60 cm and leveled. The ceramic extraction device (OD 6 cm) was installed in the center of theplot (r 5 0) into a hole of 5.1. cm OD, created beforehandwith a soil sampler, with the center of the ceramic ring at z 5
210 cm depth. The ceramic extraction device was wetted
before installation to insure the best possible contact between
the soil and ceramic. Similar tensiometers as applied in thelaboratory experiment were installed at the following positions:r 5 4.0 cm and z 5 210 cm (T 1 ); r 5 6.0 cm and z 5 215
cm (T 2 ); r 5 20.0 cm and z 5 240 cm (T 3 ); and r 5 20 cm
and z 5 255 cm (T 4 ). Vacuum was applied to the buretteusing a vacuum tank, which was evacuated at a predeterminedvacuum for each of the solution extraction step increments.
The plot was covered with a shelter to minimize temperaturefluctuations which would influence the vacuum in the tank.
The 1.2 m 3 1.2 m square plot was ponded with a constanthead of 0.5 cm water for about 3 days until the steady state
infiltration rate was 1 cm h21 (decreased from 1.3 to 0.85 cm
h21 during the 3-day period). After 3 days the total water
potential gradient (H 5 h 1 z), dH/dz, between T 3 and T 4
was approximately 1 cm cm21. Subsequently, the plot wascovered by a plastic sheet to prevent soil evaporation. The soil
was allowed to drain for a period of 46.5 h at which time thetotal head gradient was about 0.3 cm cm21. After the 46.5
hours of free drainage, the first vacuum extraction step was
applied (h ex 5 2195 cm for 46.5 , t , 71.2 h). Subsequentvacuum steps were h ex 5 2415 cm for 71.2 , t , 93.0 h, andh ex 5 2685 cm for 93.0 , t , 120 h.
From core samples at the 60-cm soil depth, measured saturated volumetric water content was 0.56 cm3 cm23 and K svaried between 2.12 and 2.94 cm h21. Dry bulk density of thesesamples varied between 1.27 and 1.35 g cm23. Time zero forthe computer simulations was at the conclusion of the infiltration test. During the extraction experiment (December 4 –9,

3.

Results and Discussion

3.1.

Numerical Experiments

Parameter optimization results for parameters a , n, u r , andK s ( u s ) of the silt, sand, and loam using four initial parameter
estimates are presented in Table 2. The daggers in the lastcolumn indicate optimized parameter values that are within
1% of the true values. Already, it must be pointed out that theparameter optimization method is tested for extreme conditions; that is, initial parameters were selected from the whole
range of soil textural classes to test convergence for any of thethree soil types. In reality, one would select as an initial estimate a parameter set close to the expected final parametervalues. In almost all optimizations for which parameter values
were different from their true values by more than 1%, simulated cumulative extraction volume and soil matric potential
values agreed well with the true values. In such cases theinverse solution converged to a local minimum with a value ofTable 2. Parameter Optimization Results for Synthetic
Cases
Final Estimates

a,
cm21

K s,
cm h21

OF

20.97
1.237
1.579
0.182

0.16e-04
0.77e-03
0.62e-03
0.20e-02

1.177
1.250
1.370
1.369

20.12
1.040
0.249
0.249

0.35e-01
0.57e-02
0.17e-04†
0.17e-04†

Loam
0.036
1.559
0.036
1.559
0.036
1.559
0.036
1.559

1.043
1.042
1.042
1.042

0.18e-04†
0.18e-04†
0.18e-04†
0.18e-04†

InitialEstimates*ur

Sandy loam
Loam
Silt
Clay

0.050
0.044
0.045
0.000

Sand
0.136
2.686
0.081
2.806
0.085
2.786
0.056
3.070

Sand
Loam
Silt loam
Clay

0.001
0.001
0.034
0.034

0.104
0.028
0.016
0.016

Sand
Sandy loam
Silt
Clay

0.078
0.078
0.078
0.078

n

Silt

*See Table 1.
†Optimized parameters are within 1% of true parameter values(Table 1).

INOUE ET AL.: IN SITU ESTIMATION

1041

Figure 4. Soil water retention and unsaturated hydraulic conductivity curves with corresponding cumulative
extraction volume (Q) and soil matric potential values (h) for (a) sand, (b) loam, and (c) silt as a function oftime during multistep extraction (T 1 , T 2 , and T 3 represent tensiometers located at (r, z) 5 (4, 18 cm), (6,18), and (6, 0.3), respectively; h, retention curve; K, hydraulic conductivity function).
the OF close to that of the global minimum. Consequently, theoptimized parameter set is nonunique and could be far fromthe real set. This was particular true for the optimization results of the sandy soil. The true hydraulic functions with thecorresponding simulated cumulative extraction volume (Q)and water potential head values (h 1 , h 2 , and h 3 ) are shown in
Figure 4 for all three investigated soils.
Although the last applied vacuum step for the sandy soil is2480 cm, corresponding soil matric potential values do notdecrease beyond 235 cm (Figure 4a), even at only 1-cm distance from the extraction device (T 1 ). Presumably, as the soil
is near u r at this water potential, the soil’s hydraulic conductivity becomes the limiting factor and much longer extraction
experiments would be needed to yield information for water
potential values smaller than 235 cm. Consequently, the optimized parameter values describe the hydraulic functions in

the limited range of 235 , h , 0 cm only, with their shapelargely controlled by the parameter n in the region. Inspectionof Tables 1 and 2 shows that all optimized n and u r values areclose to their real values of 2.68 and 0.045, respectively. Under
these hypothetical experimental conditions, the parameters K sand a appear insensitive to the OF, which can be concluded
from the good fit (OF , 1022 in Table 2) for a wide range ofvalues of these two parameters. Incidentally, the small changes
in water potential values in both space and time will require
highly accurate soil water potential measurements. Both Q andh i (i 5 1, 2, 3) responded weakly to the incremental changes
in applied vacuum. The stepwise changes in vacuum weremuch more clearly evident in the extraction volume and soil
matric potential responses for the other two soil types, where
the drop of K with matric potential (h) is much smaller thanfor the sand (Figure 4).

1042

INOUE ET AL.: IN SITU ESTIMATION

Optimization results were slightly more successful for the silt
soil (Figure 4c), for which two optimizations (using silt loam
and clay initial parameters) converged to the true solution.
Again, it is clear that n is the most sensitive parameter and K s
is the least sensitive parameter. In either case the choice ofinitial parameters is crucial to obtain convergence to the true
parameter values. The best results were obtained for the loam
soil (Figure 4b), for which convergence was achieved for any of
the four sets of initial parameter values. As Figure 4b demonstrates, soil water retention of the loam is steadily decreasingover the range of 0 to 2400 cm, whereas the rate of decreasein unsaturated K with water content is sufficient to cause aclear response of Q and h by step wise changes in applied
vacuum. Both of these characteristics make the extraction
method well suited for loamy-type soils under the given experimental conditions.
While reducing the number of optimized parameters to
three by fixing u r to its true value, none of the sandy soil
optimizations converged, whereas all the silt and loam optimizations converged to the true parameter values. These results
show that reducing the number of free parameters increases
uniqueness of the optimization problem, but at the same time
can lead to ill-posed conditions. The results in Table 2 also
show that the magnitude of the optimized K s is proportional toa . Intuitively, this is clear as large a values correspond withcoarser-textured soils of larger K s .
Parameter sensitivity for soil matric potential and extraction
volume measurements and the dependence of parameter sensitivity on tensiometer position were determined for all parameters using (11) and (13). Sensitivity of the measurement typeis optimal for maximum usu values, and differences in sensitivity
between parameters can be determined directly from comparison of s values. When comparing sensitivities to different
parameters, one should recognize that the effect of a 1%change in one parameter of the soil hydraulic functions could
be much different from that of the same 1% change of anotherparameter. Šimunek and van Genuchten [1996] showed that if
the sensitivity of matric potential is estimated as a function oftime during infiltration, parameter sensitivity is maximum
when the time rate of change of soil matric potential is high.
For suction infiltrometer measurements this occurred when
the wetting front passed the depth of the tensiometer location.
It is therefore anticipated that the sensitivity of soil matricpotential measurements will increase using step increments in
vacuum (multistep extraction) rather than using a single vacuum step, since a multitude of vacuum increments create aseries of periods with increased soil matric potential gradientsnear the extraction device. Figure 5 shows the sensitivity (s)(Equation (11)) of soil matric potential (h 2 ) at T 2 (Figure 2)
and cumulative extraction volume (Q) as a function of time forall optimized parameters for the loamy soil. When correcting
for differences in absolute magnitude between Q and h, it
becomes clear that the sensitivity of Q to any of the listed
parameters is at least 5 times as large than the sensitivity of h.
Figure 5 demonstrates that both cumulative extraction volumes and matric potentials are most sensitive to parameters aand n and least sensitive to u r and K cer. Since K s , contrary toparameters a and n, can change several orders of magnitude,
the sensitivity to its value is significant as well. Moreover,Figure 5 shows that the sensitivity to all parameters increases
with time.
Since the sensitivity is related to the time-rate of change ofthe measured flow variable, it is expected that the largest

time-averaged sensitivity (S) in (13a) be at the interface of theextraction device and the soil. However, this is not completely
so since the applied vacuum solely controls soil matric potential at that interface. In all, we find that the maximum sensitivity is near the extractor-soil interface and that sensitivity
decreases as the tensiometer is positioned at larger distancesfrom the extractor (loam and sand). However, in some cases(silt), sensitivity was high in the whole upper part of the flow
domain, which was attributed to the limited size of the flow
domain. In this part of the flow domain, soil matric potentialchanges are partly caused by the zero flux boundary condition,rather than by the applied vacuum in the extractor. The spatial
distribution of the time-averaged sensitivity of n is shown in
Figure 6 for each of the three soils.
Table 3 presents the time average and maximum sensitivity
for all four parameters and three soil types. As was alreadyanticipated from analysis of the parameter estimations, theparameter n is the most sensitive to the soil matric potential
data, and K s and u r are the least sensitive. Surprising, however,is the overall low sensitivity of all parameters for the loam
relative to that of the silt and sand, when considering theexcellent parameter optimization results of the loamy soil (Table 2). Sensitivity, as calculated by (13), reflects the behavior of
the objective function in the vicinity of the true parametervalues, that is, near the global minimum. Therefore higher
sensitivity means that the minimum is better defined andshould be estimated with higher precision once the global
minimum is identified. It does not, however, give any information about the other possible local minima in the objectivefunction (e.g., for sand and silt) elsewhere in the entire parameter space.
3.2.

Laboratory Experiment

The parameters of the soil hydraulic functions were optimized using the multistep extraction laboratory experiment.
Three sets of initial parameter values were arbitrarily selectedfrom Table 1. For each set of optimizations, Table 4 includes
the optimized parameter values as well as their uncertainty.
The latter is determined from the main diagonal of the parameter covariance matrix, and is expressed by the normalized
standard deviation (NSD), equal to the ratio of standard deviation and optimized parameter value. In this manner theuncertainties of the parameters can be evaluated by directcomparison of the NSD values. The objective function for thisexperiment included the water volume measurements at t 0 andt end, in addition to cumulative extraction volume and the matric potential values at the three tensiometer locations. Since
no saturated water content value ( u s ) was available, its value
was also optimized, bringing the total number of optimized
parameters to five. However, since the volume of water present
at t 0 was known and included in the OF, optimization of u s did
not compromise the numerical inversion. For the first set of
optimizations (“fixed K cer”), the measured value of the ceramic conductivity, K cer, was taken as a fixed parameter in the
optimization. The results of the optimizations are presented in
Figures 7a and 7b, and the corresponding optimized parametervalues are listed in Table 4. The contributions of cumulative
extraction volume, soil matric potential at three positions combined, and soil water storage to the total value of the OF arelisted separately in Table 5.
We obtained very good agreement between measured andcalculated cumulative extraction volumes for the fixed K cer
case (Figure 7a). However, it is apparent from Figure 7b that

INOUE ET AL.: IN SITU ESTIMATION

1043

Figure 5. Sensitivity (s i, j ) of (a) cumulative extraction volume (Q) and (b) soil matric potential values (h 2 )to the parameters u r , u s , a , n, K s , and K cer as a function of time for loamy soil.

the optimized soil water potential values do not match thecorresponding measured values well. We hypothesized that theindependently measured K cer was not representative for theexperimental conditions because of imperfect hydraulic contact between the ceramic and the surrounding soil. Since it is
impossible to estimate the contact area, K cer was also included
as a free parameter to be optimized. The results of these
optimizations are listed as “fitted K cer” in Table 4. Instead ofthe measured value of 0.000834 cm h21, its optimized value
was only 0.00025 cm h21. As a result of including K cer as anadditional free parameter, the value of the OF (equation (9))
decreased from 0.0753 to 0.0503 (Table 5), and parameter
uncertainty as expressed by NSD (Table 4) was reduced. However, the overall decrease in deviations in the OF was at theexpense of an increase in deviations between measured andoptimized cumulative extraction volume (see first column ofTable 5) and between measured and optimized soil matricpotential values for the first three vacuum steps. Consequently,
for the third set of optimizations, it was assumed that K cer was
equal to the measured value of 0.000834 cm h21 for the first
three vacuum increments (hydraulic contact is assumed to be
perfect) but decreased thereafter to a value to be optimized as
hydraulic contact decreases because of soil desaturation. The
results of these optimizations (“variable K cer”) are shown in
Figures 7c and 7d, and in the last column of Table 4. Convergence was obtained with K cer 5 0.00013 cm h21, corresponding
with a further decrease of parameter uncertainty (Table 4) andof the OF from 0.0503 to 0.0184 (Table 5). Although the1044

INOUE ET AL.: IN SITU ESTIMATION

Figure 6. Contour plots of time-averaged sensitivity (S(r, z, b j )) of the soil matric potential to theparameter n of multistep extraction experiment for all three soils.

cumulative extraction volume fit is not as good as for the firstcase with a fixed ceramic conductivity (compare Figures 7a and7c), an excellent fit was obtained for the soil matric potential
data for all but the last vacuum increment. Comparison of thethree cases in Table 4 shows that the increase in the optimized
K s value compensates for the corresponding reduction in K cer.
Thus the reduction in hydraulic contact area is artificially compensated for in the simulations by a decrease in the optimized
soil’s saturated hydraulic conductivity.
Close inspection of Table 5 shows that the type of weighting
employed in our optimizations results in OF values that are
dominated by soil matric potential measurements (72–98% of
total OF value). However, it is not very clear how differencesin weighting affects the OF value. For example, it is expectedthat deviations in soil matric potential values will increase if
measurement or model errors are larger. Although we havesome knowledge about measurement errors, the magnitude of
the model error and its contribution to the OF is unknown.
Also, soil heterogeneity will increase deviations between measured and simulated soil water potential values (point measurements) since the model assumes a homogeneous soil,
whereas cumulative extraction volume (whole domain measurement) is not affected by soil heterogeneity.
Table 4 also lists the soil hydraulic parameters obtained
independently by multistep outflow method with one tensiometer in the center of the core [Eching and Hopmans, 1993].
These were obtained by fitting core-average retention points tothe van Genuchten function (6). Estimated soil hydraulic functions for all three optimization options (fixed K cer, fitted K cer,and variable K cer) together with those determined independently are shown in Figure 8. The best agreement of the presented optimizations with the independently determined soil
hydraulic functions was achieved using the variable K cer optimization option. The difference in K s between the variable
K cer case and the independently estimated hydraulic functions
is about 1 order in magnitude (0.327 versus 4.2 cm h21),
whereas the other parameter values, especially the shape parameters a and n, are relatively close.
3.3.

Field Experiment

The selected diameter of the radial-symmetric computational domain for the field experiment was 70 cm. The domain
Table 3. Time-Averaged Sensitivity and Maximum
Sensitivity of Soil Matric Potential to a 1% Change of theParameters u r , a , n, and K s for the sand, silt, and loam
Sensitivity S

Soil
Type

ur

a

n

Ks

Sand
Silt
Loam

0.04 (0.093)
0.42 (0.088)
0.067 (0.493)

0.425 (1.96)
1.180 (2.74)
0.429 (1.51)

1.08 (3.04)
1.06 (2.15)
1.083 (2.28)

0.075 (0.344)
0.422 (1.029)
0.108 (0.660)

Maximum sensitivity given in parentheses.

INOUE ET AL.: IN SITU ESTIMATION

1045

Table 4. Parameter Estimation Results for Laboratory ExperimentParameter

Value

NSD, %

Value

NSD, %

Value

NSD, %

Columbia Sandy Loam,
Independently
Measureda , cm21nusur
K s , cm h21
K cer, cm h21

0.00596
2.219
0.37
0.100
0.02620.000834

4.2
6.2
0.1
16.7
15.5
zzz

0.00738
2.706
0.370
0.176
0.129
0.00025

3.5
5.8
0.07
4.3
24.4
7.9

0.0077
3.143
0.37
0.197
0.393
0.00013

1.4
3.7
0.04
1.3
19.3
7.7

0.0093
3.26
0.43
0.14
4.2
zzz

Fixed K cer

Fitted K cer

Variable K cer

NSD, normalized standard deviation, 100 s /b j (%), estimated from parameter variance.

size was large to justify the assumption of a zero flux boundary
condition along the lateral boundaries. The size of the domain
for which zero lateral flux can be assumed will depend on soil
type. It can be relatively small for sandy soils for which the fast
decrease in unsaturated hydraulic conductivity with increased
suction along the extraction device limits the size of the domain of influence but must be larger for loamy and clayey soils.
The measured decrease in soil water potential at the lower
boundary of the simulated soil domain (T 4 ) from the time thatinternal drainage start (t 0 ) was fitted with the following powerfunction:h lb 5 a~t 1 b! c 1 d(15)

The fitted power function through the measured soil water
potential data of the 55-cm-depth tensiometer is used as lower
boundary condition (h lb ) for the flow simulation. Fitted values
for the coefficients were: a 5 262.1, b 5 0.772, c 5 0.147,
and d 5 45.9, where time (t) and soil water potential (h lb )Figure 7. Comparison of measured (symbols) and optimized (lines) cumulative extraction (Q) and soil
matric potential values (h) using (a, b) a fixed K cer value and (c, d) optimized K cer (“variable K cer”).
Laboratory experiment.

1046

INOUE ET AL.: IN SITU ESTIMATION

Table 5. Contribution of Measurement Type to OF forLaboratory ExperimentContribution to OF
Case

Q

h i,
i 5 1, 2, 3

Water
Storage

Total
OF

Fixed K cer
Fitted K cer
Variable K cer

0.00130.01240.00520.07410.03790.01320.71e-07
0.986e-05
0.129e-04

0.07530.05030.0184are given in hours and cm, respectively. The value for h lb at t 0
was 213.9 cm. The measured hydraulic head gradient of 1.0
cm/cm was used to estimate h i ( z) at time t 0 . The data points
in Figure 9 correspond with the measured soil matric potential
values and cumulative extraction volume after infiltration
stopped and internal drainage started (t 0 5 0), with the firstvacuum step applied after t 5 46.5 h after which soil water
was extracted. Soil water potential values showed little response to vacuum increments, except for the tensiometer (T 1 ),which was positioned 1 cm from the extraction device. The
jumps in T 1 values correspond with the applied vacuum increments.
The soil water extractor was the same as used in the laboratory experiments with a measured K cer value of 0.0008334 cm
h21. However parameter optimizations using this value were
not very successful with relatively large differences between
optimized and measured extraction volume and soil water potential values. Hence the same approach as in the laboratory
was used, allowing K cer to be optimized. The final results areshown in Figures 9a and 9b, with the optimized parametervalues listed in Table 6 (fitted K cer). Using a K cer of 0.000282,
the fit for the field experiment was extremely good (Figures 9a
and 9b), reducing the OF value with a factor of 6 (Table 6).
From the information included into the objective function, thatis, the cumulative extracted volume and tensiometer readings,
it is not possible to estimate simultaneously both u r and u s ,since these two parameters are fully correlated. It is possible
either to fix one of these two parameters and to optimize theother or to optimize the water content interval D u 5 u s 2 u r .
Both approaches should result in similar results. We fixed u s at
the independently measured value (0.560) and estimated u r .
Note from the parameter estimation results given in Table 6that the water content interval Du was small and equal to only
0.120 cm3 cm23. The water content interval between the initialwater content and the water content corresponding to themeasured soil matric potential value at T 1 must have a close
correspondence with the cumulative extracted volume. The
small value of Du was the result of the relatively small amountof water extracted during the experiment (only about 0.5 L).
The results in Table 6 also show the large uncertainty of thesaturated hydraulic conductivity K s for which the value of NSD
is always larger than 60%.
The reduction of the optimized against measured hydraulicconductivities of the ceramic cylinder, K cer, is similar for bothlaboratory and field experiments. The ratio between optimized
K cer and K cer measured independently is 0.30 and 0.34 for thelaboratory and field experiments, respectively. K cer measurements were carried out while the ceramic was fully submergedin water. When similar measurements were performed in conductive saturated sand, the measured K cer was significantly
lower (0.000560 versus 0.0008334 cm h21).

In addition to cumulative extraction volume and soil matricpotential values at three locations, we also included the three
independently measured u (h) points during the extraction
experiment in the objective function. Inclusion of such information into the objective function breaks the mutual correlation of u r and u s so that now both parameters can be estimatedsimultaneously. Figures 9c and 9d show the comparison between measured and calculated cumulative extraction volumes
and tensiometer readings, respectively. The fit of experimental
with fitted cumulative extraction volumes and matric potential
values for T 2 and T 3 is similar to those for the previous optimizations (Figures 9a and 9b). However, deviations between
measured and fitted matric potential data for T 1 increased in
the final extraction step.
Figure 10 presents the optimized soil water retention andhydraulic conductivity functions for the fitted K cer case withand without the independently determined u (h) data, as wellas the soil hydraulic data estimated independently by the in-

Figure 8. (a) Estimated soil water retention and (b) hydraulic conductivity functions for three optimization options (Table4) for Columbia sandy loam, compared with independently
determined curves.

INOUE ET AL.: IN SITU ESTIMATION

1047

Figure 9. Comparison of measured (symbols) and optimized (lines) cumulative extraction and soil matricpotential values (h) using an optimized K cer value. Objective function includes (a), (b) cumulative extraction
volumes and soil matric potential measurements and (c), (d) three measured u (h) data points. Field experiment.
stantaneous profile method (diamond symbols). Figure 10aalso includes the three u (h) points (solid circles). Both optimized soil water retention and unsaturated hydraulic conductivity functions closely approximate the independently estimated retention and unsaturated hydraulic conductivity data in
the range between 2150 and 250 cm. However, the optimized
parameters and hydraulic functions are only valid for the soil
matric potential range from near saturation to about 2250 cm,
that is, the range over which the extraction experiment was

carried out. Continued soil water extraction after 5 days would
have generated soil matric potential data smaller than measured in the described field experiment.

4.

Summary and Conclusions

We have introduced a new method for estimating soil hydraulic parameters from a transient flow experiment. The experiment involves extraction of soil solution using successively

Table 6. Parameter Estimation Results for Field Experiment
Fixed K cer

Fitted K cer

Fitted K cer 1 u (h) data

Parameter

Value

NSD, %

Value

NSD, %

Value

NSD, %

a , cm21nusur
K s , cm h21
K cer, cm h21
OF(b)0.0331
1.746
0.560
0.389
8.932
0.0008334
0.0686

18.3
2.6
zzz
3.43
62.4
zzz
zzz

0.0231
1.688
0.560
0.441
3.28
0.000282
0.0109

23.4
2.2
zzz
3.34
74.4
3.5
zzz

0.0220
2.313
0.552
0.301
9.725
0.00026
0.0787

30.0
6.0
11.9
1.3
119.8
3.9
zzz

NSD, normalized standard deviation, 100 s /b j (%), estimated from parameter variance.

1048

INOUE ET AL.: IN SITU ESTIMATION

Figure 10. (a) Estimated soil water retention and (b) hydraulic conductivity functions for the two optimization options of Table 6 for Yolo clay loam, compared with independently determined retention andunsaturated hydraulic conductivity data.

increasing vacuum steps from an initially near-saturated soil.
The extracted volume and measured soil matric potential values at several locations near the soil water extraction device
are measured. The collected data are analyzed by a numericalcode, which combines the Levenberg-Marquardt optimization
algorithm with the unsaturated water flow code HYDRUS-2D.
The soil hydraulic parameters are obtained by minimization of
the objective function, which includes the deviations between
simulated and experimental data.

We evaluated the feasibility of the vacuum extraction technique using numerically generated data. We concluded thatthe method is well suited for loamy-textured soils but not forsandy soils. This is because the matric potential response to theapplied vacuum in the sandy soil is minimal. The success of theinversion procedure for the silt soil is dependent on initial
parameter values. Parameter sensitivity analysis showed thatthe method is most sensitive to the shape parameters a and n
and least sensitive to the residual water content u r and theINOUE ET AL.: IN SITU ESTIMATION

saturated hydraulic conductivity K s . The highest sensitivity ofthe measured soil matric potential head is close to the extraction device and decreases with increasing distance from theextractor. Therefore tensiometers should be located close to
the ceramic ring where the extraction vacuum is applied.
The method was further tested under well-defined experimental conditions in a laboratory for a Columbia sandy loam.
The objective function included the cumulative extraction volume, soil matric potential readings for three tensiometers, andthe initial and final total water volumes in the soil sample.
Parameter optimization was successful if the saturated hydraulic conductivity of the extractor (K cer) was optimized simultaneously with the soil hydraulic parameters rather than assuming its independently measured value. We hypothesize that
K cer is changing during the extraction experiment because of
reduction of hydraulic contact between the ceramic ring andthe surrounding soil as the soil desaturates.
Finally, the multistep extraction method was tested in situ
for a Yolo silt loam. Optimized soil water retention and unsaturated hydraulic conductivity data corresponded well withindependent estimates obtained from the instantaneous profile
method in the same experimental plot. However, care should
be taken in extrapolating the optimized hydraulic functions
beyond the water content range for which the experimental
data were obtained.
The problem of hydraulic contact between the ceramic
membrane and the surrounding soil is of critical importancefor further applications of the extraction method. From ouranalysis it appears that although we took great care to assure
hydraulic contact when installing the extraction device, hydraulic contact was reduced during soil water extraction, thereby
affecting the optimization results. To overcome this problem,
improved devices such as those presented by Shani and Or[1995], which guarantee hydraulic contact with the surrounding
soil throughout the duration of the extraction experiment,
need to be developed. Moreover, similar to the experiments by
Gribb [1996], the experimental procedure would be greatly
simplified if the tensiometers and extraction device were combined into a single probe.
Acknowledgments. This research was supported by a grant fromthe Water Resources Center. We thank Jim MacIntyre for his tirelesswork assisting in the experimental part of the study. We also acknowledge Rien van Genuchten of the USDA-ARS Salinity Laboratory forhis support. The authors also thank Brent Clothier and another anonymous reviewer for their constructive criticism.

References
Bard, Y., Nonlinear Parameter Estimation, 341 pp., Academic, San
Diego, Calif., 1974.
Brooks, R. H., and A. T. Corey, Properties of porous media affecting
fluid flow, J. Irrig. Drain. Div. Am. Soc. Civ. Eng., 92, 61– 88, 1966.
Carrera, J., and S. P. Neuman, Estimation of aquifer parameters undertransient and steady state conditions, 1, Maximum likelihood
method incorporating prior information, Water Resour. Res., 22,
199 –210, 1986a.
Carrera, J., and S. P. Neuman, Estimation of aquifer parameters undertransient and steady state conditions, 2, Uniqueness, stability, andsolution algorithms, Water Resour. Res., 22, 211–227, 1986b.
Carsel, R. F., and R. S. Parrish, Developing joint probability distributions of soil water retention characteristics, Water Resour. Res.,
24(5), 755–769, 1988.
Celia, M. A., E. T. Bouloutas, and R. L. Zarba, A general massconservative numerical solution for the unsaturated flow equation,
Water Resour. Res., 26(7), 1483–1496, 1990.

1049

Clausnitzer, V., and J. W. Hopmans, Non-linear parameter estimation:
LM_OPT. General-purpose optimization code based on the Levenberg-Marquardt algorithm, Land, Air Water Resour. Pap. 100032,
Univ. of Calif., Davis, 1995.
Dirksen, C., Unsaturated hydraulic conductivity, in Soil Analysis: Physical Methods, edited by K. A. Smith and C. E. Mullins, pp. 209 –269,
Marcel Dekker, New York, 1991.
Doering, E. J., Soil water diffusivity by the one-step method, Soil Sci.,
99, 322–326, 1965.
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(5), 1167–1175, 1993.
Eching, S. O., J. W. Hopmans, and O. Wendroth, Unsaturated hydraulic conductivity from transient multi-step outflow and soil waterpressure data, Soil Sci. Soc. Am. J., 58, 687– 695, 1994.
Gardner, W. R., Calculation of capillary conductivity from pressure
plate outflow data, Soil Sci. Soc. Am. Proc., 20, 317–320, 1956.
Green, R. E., L. R. Ahuja, and S. K. Chong, Hydraulic conductivity,
diffusivity, and sorptivity of unsaturated soils: field methods, in
Methods of Soil Analysis, part I, monogr. 9, edited by A. Klute, 2nded., pp. 771–798, Soil Sci. Soc. of Am., Madison, Wis., 1986.
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.
Klute, A., Water retention: Laboratory methods, in Methods of Soil
Analysis, part I, monogr. 9, edited by A. Klute, 2nd ed., pp. 635– 662,
Soil Sci. Soc. of Am., Madison, Wis., 1986.
Klute, A., and C. Dirksen, Conductivities and diffusivities of unsaturated soils, in Methods of Soil Analysis, part I, monogr. 9, edited by
A. Klute, 2nd ed., pp. 687–734, Soil Sci. Soc. of Am., Madison, Wis.,
1986.
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 for 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.
Marquardt, D. W., An algorithm for least-squares estimation of nonlinear parameters, SIAM J. Appl. Math., 11, 431– 441, 1963.
Mualem, Y., A new model for predicting the hydraulic conductivity ofunsaturated porous media, Water Resour. Res., 12, 513–522, 1976.
Parker, J. C., J. B. Kool, and M. T. van Genuchten, Determining soil
hydraulic properties from one-step outflow experiments by parameter estimation, II, Experimental studies, Soil Sci. Soc. Am. J., 49,
1354 –1359, 1985.
Shani, U., and D. Or., In situ method for estimating subsurface unsaturated hydraulic conductivity, Water Resour. Res., 31, 1863–1870,
1995.
Šimunek, J., and M. T. van Genuchten, Estimating unsaturated soil
hydraulic properties from tension disc infiltrometer data by numerical inversion, Water Resour. Res., 32(9), 2683–2696, 1996.
Šimunek, J., and M. T. van Genuchten, Parameter estimation of soil
hydraulic properties from multiple tension disc infiltrometer data,
Soil Sci., 6, 383–398, 1997.
Šimunek, J., M. Sejna, and M. T. van Genuchten, The HYDRUS-2D
software package for simulating water flow and solute transport in
two-dimensional variably saturated media, version 1.0, IGWMCTPS-53, Int. Ground Water Modeling Cent., Colo. School of Mines,
Golden, 1996.
Toorman, A. F., P. J. Wierenga, and R. G. Hills, Parameter estimationof soil hydraulic properties from one-step outflow data, Water Resour. Res., 28, 3021–3028, 1992.
Valiantzas, J. D., and D. G. Kerkides, A simple iterative method forthe simultaneous determination of soil hydraulic properties fromone-step outflow data, Water Resour. Res., 26, 143–152, 1990.
van Dam, J. C., J. N. M. Stricker, and P. Droogers, From one-step tomulti-step determination of soil hydraulic functions by outflow experiments, Rep. 7, Dep. of Water Resour., Agric. Univ., Wageningen, Netherlands, 1990.
van Dam, J. C., J. N. M. Stricker, and P. Droogers, Inverse method for
determining soil hydraulic functions from multi-step outflow experiments, Soil Sci. Soc. Am. J., 58, 647– 652, 1994.
van Genuchten, M. T., A closed-form equation for predicting the1050

INOUE ET AL.: IN SITU ESTIMATION

hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44,892– 898, 1980.
Whisler, F. D., and K. K. Watson, One-dimensional gravity drainage of
uniform columns of porous materials, J. Hydrol., 6, 277–296, 1968.
Yeh., W. W.-G., Review of parameter identification procedures ingroundwater hydrology: The inverse problem, Water Resour. Res.,
22(2), 95–108, 1986.
V. Clausnitzer and J. W. Hopmans, Hydrology Program, Dept.
LAWR, University of California, Davis, CA 95616. (e-mail:jwhopmans@ucdavis.edu)

M. Inoue, Arid Land Research Center, Tottori University, Tottori,
Japan.
J. Šimunek, USDA Salinity Laboratory, USDA-ARS, 450 Big
Springs Road, Riverside, CA 92507.

(Received July 14, 1997; revised January 15, 1998;
accepted January 22, 1998.)