WATER RESOURCES RESEARCH, VOL. 41, W04015, doi:10.1029/2004WR003379, 2005

Estimating soil hydraulic parameters from transient flow
experiments in a centrifuge using parameter optimization
technique
Jirka Šimůnek
Department of Environmental Sciences, University of California, Riverside, California, USA

John R. Nimmo
U.S. Geological Survey, Menlo Park, California, USA
Received 28 May 2004; revised 20 January 2005; accepted 7 February 2005; published 27 April 2005.

[1] A modified version of the Hydrus software package that can directly or inversely

simulate water flow in a transient centrifugal field is presented. The inverse solver forparameter estimation of the soil hydraulic parameters is then applied to multirotationtransient flow experiments in a centrifuge. Using time-variable water contents measured at
a sequence of several rotation speeds, soil hydraulic properties were successfully
estimated by numerical inversion of transient experiments. The inverse method was thenevaluated by comparing estimated soil hydraulic properties with those determined
independently using an equilibrium analysis. The optimized soil hydraulic properties
compared well with those determined using equilibrium analysis and steady stateexperiment. Multirotation experiments in a centrifuge not only offer significant timesavings by accelerating time but also provide significantly more information for theparameter estimation procedure compared to multistep outflow experiments in a
gravitational field.
Citation: Šimůnek, J., and J. R. Nimmo (2005), Estimating soil hydraulic parameters from transient flow experiments in a centrifuge
using parameter optimization technique, Water Resour. Res., 41, W04015, doi:10.1029/2004WR003379.

1. Introduction[2] The main information needed for evaluating, analyzing, and predicting flow and solute transport in soils
concerns the soil hydraulic properties. Traditional methodsof measuring soil hydraulic properties are tedious and time
consuming [Dane and Topp, 2002]; therefore faster, easiertechniques are constantly being sought [e.g., Hopmans etal., 2002]. One way of accelerating the measurement
process is by carrying out needed experiments in a centrifuge. The steady state conditions required for many types ofmeasurements can be achieved in much shorter time [e.g.,
Nimmo et al., 1987, 2002; Nimmo and Mello, 1991; Concaand Wright, 1992].
[3] Over the last century, high-speed centrifuges havebecome standard in scientific and technical fields of soil
physics, petroleum industry and environmental engineering
for measuring water contents, and saturated and unsaturated
hydraulic conductivities. Briggs and McLane [1907] used acentrifuge to measure what they called the ‘‘moisture equivalent,’’ a single-number index of a soil’s capacity to retain
water. Later, the water retention curve was the measurement
objective, for example, in the work of Gardner [1937],
Russell and Richards [1938], and Hassler and Brunner[1945]. In these methods the water in an unsaturated sample
equilibrates with the centrifugal field, much as in the
hanging-column method [Dane and Hopmans, 2002] it
Copyright 2005 by the American Geophysical Union.
0043-1397/05/2004WR003379$09.00

equilibrates with a gravitational field. Development of suchtechniques has continued to the present [Forbes, 1994;
Paningbatan, 1980; Sigda and Wilson, 2003]. The use ofcentrifuge-scaled physical modeling of dynamic unsaturated
flow has also been explored recently [e.g., Culligan et al.,
1997; Nakajima et al., 2003]. For measuring unsaturated
hydraulic conductivity by the steady state centrifuge method,
two different types of apparatus have been used, as explainedby Nimmo et al. [2002]. In the internal flow control (IFC)
apparatus the means of controlling flow is internal tothe centrifuge bucket [Nimmo et al., 1987], while in theunsaturated flow apparatus (UFA), the flow is controlled
outside the centrifuge [Conca and Wright, 1992, 1998]. Eachof these implementations has its own advantages anddisadvantages, for example larger samples and greater
speeds with the IFC apparatus, faster operation and lessdemand on the operator with the UFA apparatus [Nimmo etal., 2002]. Some experiments have investigated transient
instead of steady flow [e.g., Alemi et al., 1976; Nimmo,
1990], and even unstable flow [Culligan et al., 1997;
Griffioen et al., 1997]. So far these transient methods
have not led to widely used techniques for soil property
measurement, though they provide significant insights into
the nature of unsaturated flow.
[4] Recent advances allow measurement of state variables
inside the rotating sample. Nimmo [1990] measured electricalconductivity using a four-electrode technique and converted
the measured values into water contents. Savvidou andCulligan [1998] monitored pore pressures and solute concentrations in a spinning sample with miniature pressure

W04015

1 of 9

W04015ŠIMŮNEK AND NIMMO: TRANSIENT ANALYSES OF FLOW EXPERIMENTS IN A CENTRIFUGE

transducers and resistivity probes. Investigators at the Idaho
National Engineering and Environmental Laboratory areexperimenting with measuring pressure heads using
miniature tensiometers rotating together with soil samples[Mattson et al., 2003].
[5] Direct solvers of the Richards equation modified toaccount for a centrifugal field have been developed in thelast two decades. Bear et al. [1984] presented a theory of
fluid flow through a deformable porous medium in thepresence of a centrifugal force and then solved the resulting
flow and deformation equations using the finite difference
method. They discussed implementation of several types ofboundary conditions for different applications and presented
calculated water content and pressure head profiles for
hypothetical simulations. In their simulations they showedthat the porosity changes were less than 0.01% after
12 minutes of centrifuging, and thus can be neglected in
majority of applications. Nimmo [1990] also developed adirect solver of the Richards equation for a centrifugal
field and successfully compared numerical simulations
using independently determined soil hydraulic properties
to measured results of transient flow experiments.
[6] An alternative to steady state experiments is analysisof transient experiments using parameter estimation techniques [Hopmans et al., 2002]. Steady state calculations areoften based on assumption of constant water contents andzero, or linear, pressure head gradients, which is not always
the case. On the other hand, parameter estimation techniques, such as those based on the numerical solution of the
Richards equation, do not require such assumptions, and are
constrained only by assumptions inherent in the derivationof the transient Richards equation. The one-step outflow
method [Gardner, 1956; Kool et al., 1985] and themultistep outflow method [van Dam et al., 1994; Echingand Hopmans, 1993; Durner et al., 1999] are frequently
used in combination with a parameter estimation technique
to determine soil hydraulic parameters in a gravitational field[Hopmans et al., 2002]. During the experiment, pressures orsuctions are applied at the top or bottom of the soil sample
while the resulting outflow from the sample, and, optionally,
also the pressure head within the sample, are measured.
Cumulative outflow data in combination with measured
pressure heads, if available, are then used to define theobjective function against which the numerical model iscalibrated to estimate soil hydraulic parameters.
[7] Similar procedures can be followed in the centrifuge
using step changes in rotational speed rather than appliedpressure head. Different rotation speeds then, when thelower boundary condition is specified by the zero pressurelevel at an imposed water table established within a porous
plate below the soil boundary, will result in different
suctions applied at the bottom of the sample. Similarly tothe multistep outflow experiment in a gravitational field,
sequentially increasing rotational speeds will result in
stepwise draining of the sample. This information can beused to define the objective function for model calibration.
Alemi et al. [1976] presented a theory for analyzing thevolumetric outflow of water from a soil core when the speedof centrifugation is suddenly increased. Their approach,however, was not based on numerical inversion of the
Richards equation, and they did not present or analyze
experimental data.

W04015[8] When a multistep outflow experiment in a gravitational field is carried out for a long enough time at eachpressure step to reach steady conditions, the resulting pairsof water contents and pressure heads define data points of
the retention curve. Similar experiments in a centrifugal
field can provide much more information. Since a centrifugal field results in pressure heads that vary greatly between
the top and bottom of the sample, each steady state
condition can provide multiple data points of the retention
curve depending on the number of positions at which thewater content is measured.
[9] In this manuscript we address the question of whether
one can obtain additional benefits by combining theadvantages of centrifuge-accelerated experiments withparameter estimation analysis. The objectives of this paper
are threefold. The first is to present direct and inverse
solvers of the Richards equation modified for a transient
centrifugal field. Second, to apply the inverse solver forparameter estimation of the soil hydraulic parameters using
multirotation transient flow experiments in a centrifuge.
Third, to evaluate the method by comparing estimated soil
hydraulic properties with those determined independently
using an equilibrium analysis.

2. Theory
2.1. Mathematical Model[10] Darcy-Buckingham’s law for one-dimensional flow
in a centrifugal field, neglecting gravity and hysteresis, canbe expressed as [Nimmo et al., 1987; Conca and Wright,
1992]
q ¼ K*ðyÞdy rw2 rdr

ð1Þ

where q is the Darcy flux [LT1], y is the matric potential[ML1 T2, Pa], r is the bulk density of water [ML3], w is
the angular speed [radT1], K* is the hydraulic conductivity[L3TM1], and r is the radial coordinate [L]. This
formulation leads to units for the hydraulic conductivity
that are not often used. Simple transformation of theequation leads to conductivity being expressed in morecommonly used units of [LT1] and Darcy-Buckingham’s
law being expressed in terms of the pressure head:
q ¼ K*rg

dy
rw2 r

rgdr
rg




¼ K

dh w2 r

dr
g


ð2Þ

where K(=K*rg) is the hydraulic conductivity in [LT1], g
is the gravitational acceleration [LT2] (i.e., 9.81 m s2) andh is the pressure head [L].
[11] While the Richards equation in a unit gravitational
field has the following formulation
 

@q
@
@h
¼
K
þ cos a
@t @x
@x

ð3Þ

where q is the water content [L3 L3], and a is the angle
between gravity and flow direction (0 for vertical flow), the2 of 9

W04015ŠIMŮNEK AND NIMMO: TRANSIENT ANALYSES OF FLOW EXPERIMENTS IN A CENTRIFUGE

modified version of the Richards equation for a centrifugal
field is:
 

@q
@
@h w2 r
¼
K

@t @r
@r
g

ð4Þ

W04015

while the lower boundary condition discretized Darcy’s law:dh w2 r
q1 ¼ K

dr
g



jþ1;k
 K1þ1=2

w2 r1þ1=2
h2jþ1;kþ1  hjþ1;kþ1
1

Dr
g

!

ð8Þ

2.2. Numerical Model[12] The numerical model for flow in a centrifugal field
was developed by modifying the Hydrus-1d code [Šimůnek
et al., 1998b]. The Richards equation (3) in Hydrus issolved numerically using a Galerkin type linear finite
element method combined with a mass lumping scheme,
leading to a numerical scheme equivalent to standard finite
differences. Integration in time is achieved using an implicit(backward) finite difference scheme for both saturated andunsaturated conditions. The water content term is evaluatedusing the mass conservative method proposed by Celia etal. [1990] that leads to the set of nonlinear algebraic
equations:q ijþ1;kþ1  q ji
Dt¼

jþ1;kþ1
jþ1;kþ1 hjþ1;kþ1 hjþ1;kþ1
jþ1;k hiþ1
jþ1;k hiii1Kiþ1=2
 Ki1=2
Dxi
Dxi1

1
Dx
þ

jþ1;k
jþ1;k
Kiþ1=2
 Ki1=2

Dx

After taking into account initial and boundary conditions,these are solved using the Picard iteration scheme. In (5), iis the node number in the spatial discretization, j is the index
for the temporal discretization, k is the iteration number in
the Picard iteration scheme, Dt is the time step [T], and Dx[=(xi+1  xi1)/2], Dxi [=xi+1  xi], and Dxi1 [=xi  xi1]
are spatial discretization steps [L].
[13] The modified Richards equation for flow in a centrifugal field (4) is discretized in a similar manner, leadingtojþ1;kþ1
jþ1;kþ1q ijþ1;kþ1  q ji
1 hjþ1;kþ1 hjþ1;kþ1
jþ1;k hiþ1
jþ1;k hiii1¼
Kiþ1=2
 Ki1=2
Dr
Dr
Dt
Dr
i
i1

2



w
g

Kijþ1;k þ ri

jþ1;k
jþ1;k !
Kiþ1=2
 Ki1=2

Dr

!

ð6Þ

where Dr [=(ri+1  ri1)/2], Dri [=ri+1  ri], and Dri1 [=ri ri1] are spatial discretization steps [L].
[14] Numerical solution requires adaptive time steps
and spatially variable discretization. Discretization at theoutflow boundary needs to be significantly finer to accommodate large gradients at early times after a change in therotation speed and to result in small mass balance errors.
[15] Special attention needs to be paid to flux boundary
conditions, since they have an important effect on stabilityof the numerical solution, as well as on the overall massbalance error. The upper boundary condition was specified
by discretization of the Richards equation as follows:q Njþ1;kþ1  q jN
2¼
Dt
DrN 1
þ KNjþ1;k
1=2

jþ1;k
qjþ1
N  KN 1=2

w2 rN 1=2
g!

2.3. Soil Hydraulic Properties[16] The unsaturated soil hydraulic properties in thispaper are assumed to be represented by the following
expressions [van Genuchten, 1980]:qðhÞ  qr
1¼
m
qs  qrð1 þ jahjn Þ

ð9Þ

h

m i2
K ðqÞ ¼ Ks Sel 1  1  Se1=m

ð10Þ

Se ðhÞ ¼

!

ð5Þ

cos a

where qN and q1 are the specified upper and lower boundary
fluxes [LT1] and indices N and 1 represent boundary
nodes.

hNjþ1;kþ1  hNjþ1;kþ1
1
DrN 1
ð7Þ

where Se is the effective water content (dimensionless), Ks is
the saturated hydraulic conductivity [LT1], qr and qs denoteresidual and saturated water contents [L3 L3], respectively,
l is a pore connectivity parameter (dimensionless), and a[L1], n (dimensionless), and m(=1  1/n) (dimensionless)are empirical parameters. The predictive K(q) model isbased on the capillary model of Mualem [1976] inconjunction with equation (9). The pore connectivity
parameter l in the hydraulic conductivity function was
estimated by Mualem [1976] to be 0.5 as an average formany soils. The hydraulic characteristics defined byequations (9) and (10) contain 5 unknown parameters: qr,
qs, a, n, and Ks. For a multirotation experiment carried out asa drying process, the hydraulic parameters in Equations (9)and (10) represent drying branches of the unsaturated
hydraulic properties.
2.4. Parameter Optimization[17] The soil hydraulic parameters of (9) and (10) are
obtained by calibrating the numerical model by minimizing
the objective function F defined as the sum of squared
deviations between measured and calculated water contents:
Fðb; qÞ ¼

nj
m X
X

q*j ðti Þ  qj ðti ; bÞ

2

ð11Þ

j¼1 i¼1

where m represents different locations of water content
measurements, nj is the number of measurements at aparticular location, q*(tj i) are specific water content measurements at time ti and location j, and qj(ti, b) are thecorresponding model predictions for the vector of optimized
parameters b (e.g., qr, qs, a, n, and Ks).
[18] Minimization of the objective function F is accomplished by using the Levenberg-Marquardt nonlinear minimization method [Marquardt, 1963], which has become acommonly used approach in nonlinear least squares fitting
among soil scientists and hydrologists [Šimůnek andHopmans, 2002]. The method combines the Newton andsteepest descend methods, and provides confidence intervals for the optimized parameters. Since the Levenberg-

3 of 9

W04015ŠIMŮNEK AND NIMMO: TRANSIENT ANALYSES OF FLOW EXPERIMENTS IN A CENTRIFUGE

Figure 1. Schematic of the soil sample, location ofelectrodes (Arabic numberals), and assumed locations forwater content measurements (Roman numerals). Water
contents for location I were calculated by analyzing signal
from electrodes 1, 2, 3, and 4. Similarly, location II
corresponds to electrodes 3, 4, 5, and 6; location III
corresponds to electrodes 5, 6, 7, and 8; location
IV corresponds to electrodes 1, 2, 5, and 6; location V
corresponds to electrodes 3, 4, 7, and 8; and location VI
corresponds to electrodes 1, 2, 7, and 8. See color versionof this figure in the HTML.
Marquardt method is a gradient method that seeks only the
local minimum, we restarted each optimization with several
different initial estimates of optimized parameters, andreport here those that resulted in the lowest value of theobjective function.
[19] Although in the presented analyses we assume
that measured water contents are point measurements, thenumerical code can also consider a finite measurement
volume by averaging calculated water contents over
specified volume. While in this paper we define the objective function only in terms of the water content, the Hydrus
model allows other types of information, such as pressure
heads and actual or cumulative outflow, to be considered aswell.

near the bottom of the ceramic plate (away from the axis ofrotation). For different rotation speeds the water table at thisposition established different constant pressure heads at theoutflow boundary of the soil sample.
[21] Eight electrodes in two vertical rows were installed
through the retainer into the soil sample. This configuration
allowed electrical conductivity measurements using
six different four-electrode sets (Figure 1). Three main
measurements were obtained by using four neighboring
electrodes (small square configuration), and thus electrodes 1,
2, 3, and 4 composed the lower electrode set (measurements I); 3, 4, 5, and 6 the middle set (measurements II);and 5, 6, 7, and 8 the upper set (measurements III). Threeother combinations of electrodes (larger rectangle configuration) could also be analyzed, resulting in measurement
sets IV (electrodes 1, 2, 5, and 6), V (electrodes 3, 4, 7, and8), and VI (electrodes 1, 2, 7, and 8). Electrical conductivity
observations are thus at mean depths of 2.67, 1.905, 1.141,
2.279, 1.520, 1.905 cm for measurement sets I through VI.
Nimmo [1990] conducted supplementary measurements
showing that about 90% of the measurement influence isconfined to a region within 9 mm of the center of eachsquare. The published results of these supplementarymeasurements can be used in evaluating the uncertainty
that results from assuming that data measured with thesquare electrode array represent point measurements at thecenter of each square. The calibration procedure to obtainwater contents from measured electrical conductivity
measurements is described by Nimmo [1990].
[22] Experiments were carried out on Oakley sand packed
to a bulk density of 1.83 Mg/m3 and porosity 0.333. The
soil hydraulic properties, i.e., the hydraulic conductivities
and retention data were measured using the steady state
centrifuge method by Nimmo and Akstin [1988]. The soil
sample had a height of 38 mm and a diameter of 25.4 mm.
The soil was saturated before experimental runs using
selenate solution (0.01 N CaSO4 and 0.01 CaSeO4) withelectrical conductivity equal to 1.43 dS/m at the laboratory
temperature of 22C.
[23] The sequence of rotation speeds followed in three
separate runs from resaturation is given in Table 1. Table 1
also indicates the centrifugal force for the applied rotation
speed at the center of the soil sample. Note that in our

Table 1. Sequences of Centrifuge Speedsa
Run 1

3. Materials, Methods, and Data
[20] Details about the centrifuge apparatus, as well as
about the soil sample used in this study are given by Nimmo[1990]. Here we only summarize the main features ofcompleted experiments and their corresponding analyses.
The centrifuge was specially modified by the manufacturer,
Beckman Instruments, to incorporate a Meridian Laboratory
MM12 high-speed, low-noise rotating electrical contractor
on top of the rotor. (Brand names are given for identification
only and do not indicate any endorsement by UCR or
USGS.) The soil sample was placed in a retainer on a
porous ceramic plate 9.2 mm thick (Figure 1). The effective
hydraulic conductivity of a plate was measured to be 3 107 m s1 [Nimmo, 1990]. A short lip of the impervious
plate supporting the ceramic plate maintained a water table

W04015

1

Run 2
1

t, s

w, s

g

t, s

w, s

10,24964,12972,35380,352
92,299144,744
169,222

46.1
57.6
70.2
74.3
104.7
136.1
230.4

38.1
59.5
88.3
99.2
197
333
953

5
40
65
90
4,389
9,028
13,73725,21975,87995,969115,576

50.3
33.8
24.1
16.8
14.7
22.0
29.3
38.7
46.1
57.6
104.7

Run 3
g

t, s

w, s1

g

45.5
20.5
10.4
5.1
3.9
8.7
15.4
26.9
38.1
59.5
197

10
15
40
60
6,393
10,345
17,82533,067
83,677
93,988104,234
122,717379,884

6.3
50.3
45.0
35.6
29.3
33.5
36.7
38.7
41.9
46.1
52.4
57.6
70.1

0.71
45.5
36.4
22.8
15.4
20.2
24.1
26.9
31.5
38.1
49.2
59.5
88.3

a
Here t, final time for given rotations; w, rotation speed (rad s1); g,
multiple g acceleration.

4 of 9

W04015ŠIMŮNEK AND NIMMO: TRANSIENT ANALYSES OF FLOW EXPERIMENTS IN A CENTRIFUGE

Table 2. Soil Hydraulic Parameters Obtained From Equilibrium
Analysesa
Run

qr

qs

a, m1

n

R2

1
2
3
2a
3a

0.047 (0.0023)
0.00
0.015 (0.021)
0.062 (0.0049)
0.051 (0.011)

0.280b0.280b0.280b0.180 (0.0035)
0.200 (0.017)

3.28 (0.49)
19.6 (4.9)
3.45 (0.38)
1.22 (0.082)
1.52 (0.26)

2.03 (0.11)
1.30 (0.029)
1.76 (0.17)
3.45 (0.47)
2.65 (0.47)

0.94
0.82
0.93
0.95
0.93

a 2

R is the regression coefficient. Standard errors are given in parentheses.
Fixed.

bcalculations we could take into account not only the maintime segments when a particular rotation speed was applied
over a period of several hours, but also short time segments
lasting only several seconds or minutes.

4. Results[24] The collected experimental data can be analyzed inseveral ways. In the following sections we will first apply asimple equilibrium analysis and then the more complex
inverse analyses using parameter estimation in combination
with a numerical solver of the unsaturated flow equation.
4.1. Equilibrium Analyses[25] First, we analyze measured water contents at the endof each main time segment, assuming that hydrostatic
equilibrium was achieved before the next change of speed.
We do realize that especially at greater rotation speeds, it isunlikely that the equilibrium conditions were fully reached.
In the low-water-content conditions of these runs, however,q is very weakly sensitive to h, so that deviations fromequilibrium are minor, and thus errors from this assumption
are minor as well. Numerical evaluation of equation (1)shows that equilibrium conditions are reached relatively
quickly (within minutes) at low rotations, while at large
rotation speeds the equilibration process can last many days
[Nimmo, 1990] (see section 4.2).
[26] Under the assumption of equilibrium conditions,pressure heads for particular positions, r, can be calculated
from [Nimmo, 1990]:
hðrÞ ¼yðrÞ w2  2¼
r  r02
rg
2g

ð12Þ

W04015where r0 is the radius of location with zero potential [L](=20.4 cm). These calculated pressure heads were thenpaired with corresponding measured water contents at the
six locations I through VI (Figure 1) at the end of aparticular time segment. The resulting retention data werefitted with the analytical model of van Genuchten [1980]
using the RETC code [van Genuchten et al., 1991]. Twosets of optimizations were done on each data set. First, we
optimized three retention parameters qr, a, and n, and fixed
the saturated water content qs to a value reported by Nimmo[1990]. Then we fitted all four retention parameters
including qs. The resulting soil hydraulic parameters forboth optimizations are given in Table 2, while the graphs ofboth the measured retention data and their fits are shown in
Figure 2. Note that since for the first run the first rotation
speed was relatively high and thus the first retention data
point was for the pressure head of 80 cm, we could not
successfully fit the retention model without fixing thesaturated water content. For the two other runs, with the firstrotation speeds of 14.7 and 29.3 s1, the largest pressure
heads were equal to 9 and 35 cm, respectively, and thus
the retention curve was much better defined close tosaturation, resulting in better fits. For run 2, however,measured water contents close to saturation for different
locations display relatively large differences.
4.2. Transient Analyses[27] Second, we analyzed transient data using parameter
estimation in combination with a numerical solver ofthe unsaturated flow equation (1). Optimizations for eachcentrifugal run were again run in two sets. First, with all fivesoil hydraulic parameters optimized, and second with theeffective saturated water content fixed at the observed valueof 0.28 [Nimmo, 1990]. Figure 3 compares measured andfitted water contents for all three runs and for all 6
measurement locations. While measured water contents
seem to be reaching the equilibrium state for slower speeds,
as expected, for larger speeds water contents were stilldecreasing at the end of the run, and thus equilibrium wasclearly not reached. In these cases, however, the measured
water contents may not be significantly different from thevalues they would have at a true equilibrium state. Noticethat there are significant differences in water contents within
the soil sample at each speed, with the difference between
the minimum and maximum water contents often larger

Figure 2. Retention curves resulting from the equilibrium analyses. The saturated water contents were
either fixed (denoted Fitted) or fitted (Fitted a). See color version of this figure in the HTML.
5 of 9

W04015ŠIMŮNEK AND NIMMO: TRANSIENT ANALYSES OF FLOW EXPERIMENTS IN A CENTRIFUGE

W04015Figure 3. Measured and optimized water contents for three centrifuge runs 1, 2, and 3. See color
version of this figure in the HTML.
than 0.05. That means that at each rotational speed theoutflow experiment in a centrifugal field, especially atintermediate and larger speeds, provides more information
than the similar experiment in the gravitational field (inwhich a much smaller range of water content is represented).
[28] The numerical model was relatively successful indescribing the pattern in measured data, both the earlystages of water content changes after the change in therotational speed, as well as reaching of the equilibrium stage(Figure 3). Agreement is better for the measurements at
higher speeds, consistent with the observation by Nimmo[1990] that various factors cause greater experimental
uncertainty in the lower-speed data. The optimized soil
hydraulic parameters and their standard errors are given in
Table 3 together with regression coefficients of the fit.

5. Discussion[29] Figure 4 shows calculated pressure head and watercontent profiles for all three runs at times when rotational
speed was changed. A parabolic decrease of the pressure
head from the bottom of the sample toward the top would
represent an equilibrium status. The pressure head is
expected to decrease parabolically within the porous ceramic,
which is assumed to be fully saturated at all times; thus
propagation of pressure should be relatively quick. Visually
it seems that equilibrium pressure head profiles were notreached only for the last three rotational speeds for run 1(the three largest speeds) (Figure 4a), and for the last two
speeds for run 2 (Figure 4b), for which the duration was
relatively short (Figure 3, middle). A large increase inrotational speed may lead to sudden desaturation of thebottom of the soil sample. Associated decrease of thehydraulic conductivity in the soil close to the interface

may establish a significant resistance to draining of thesample and reaching of equilibrium conditions.
[30] The comparison of transient optimized retentionfunctions with water contents obtained under the assumption of equilibrium state at the end of each rotation speed for
all three runs is shown in Figure 5. Notice that there is avery good correspondence between optimized retention
curves for intermediate pressure heads between 0.5 andabout 3 m. There is a significant scatter of water content
values close to saturation, probably caused by slightly
different initial conditions due to air entrapment during
resaturation. By applying substantial acceleration from thevery beginning, little information is collected for the experimental range close to saturation. This experimental range isprone to large uncertainties even for experiments in thegravitational field [Hopmans et al., 2002] and accelerationof the experiment in this range is usually not recommended[Wendroth et al., 1993; Šimůnek et al., 1998a]. Unless thesaturation status is measured from the very beginning of
the experiment or the mass balance information, such as the
cumulative outflow, is collected (which was not the case forour experimental data), the fitted parameters that characterize the soil hydraulic properties close to saturation (that is,
the saturated water content and the saturated hydraulicconductivity) are subject to large uncertainties. Since there
are only very few water contents measurements with values
above 0.2, the optimized soil hydraulic parameters aremainly determined by the soil hydraulic properties in theintermediate and dry saturation range. The saturated watercontent and the saturated hydraulic conductivity values
should be viewed mainly as fitting parameters without clear
physical meaning (see also Schaap and Leij [2000] and their
discussion on the physical meaning of the fitted saturated
hydraulic conductivity). Simultaneous fitting of these two

Table 3. Soil Hydraulic Parameters Obtained From Transient Analysesa
Run

qr

qs

a, m1n

Ks, m s1

R2

1
2
3
1a
2a
3a

0.043 (0.00052)
0.055 (0.00081)
0.061 (0.00034)
0.040 (0.00065)
0.000 (0.0028)
0.038 (0.001)

0.220 (0.0025)
0.184 (0.00044)
0.191 (0.00078)
0.280b0.280b0.280b

1.76 (0.051)
1.19 (0.0081)
1.39 (0.010)
2.07 (0.094)
11.5 (0.30)
3.14 (0.031)

2.38 (0.038)
3.69 (0.063)
3.38 (0.038)
2.13 (0.032)
1.39 (0.0082)
1.98 (0.017)

7.30e-07 (6.10e-08)
1.00e-07 (6.0e-09)
2.29e-07 (1.65e-08)
3.87e-06 (1.8e-07)
2.31e-05 (2.5e-06)
1.01e-05 (5.2e-07)

0.92
0.94
0.96
0.91
0.81
0.94

R is the regression coefficient. Standard errors are given in parentheses. Read 7.30e-07 as 7.30  107.
Fixed.

a 2b

6 of 9

W04015ŠIMŮNEK AND NIMMO: TRANSIENT ANALYSES OF FLOW EXPERIMENTS IN A CENTRIFUGE

W04015Figure 4. (top) Pressure heads and (bottom) water contents at times when rotation speed was changed:(a) run 1, (b) run 2, and (c) run 3. See color version of this figure in the HTML.
parameters together with the shape parameters, however, is
recommended [Hopmans et al., 2002] since it provides
numerical code with additional flexibility to describe soil
hydraulic properties in the measuring range for which mostinformation is collected, i.e., intermediate and dry range.
[31] There are slight differences also at the dry end, where
the optimization routine returned values of qr in the range of
0.043 and 0.061. In run 1, a much greater rotation speed
was used for the last step (see Table 1), and thus muchsmaller pressure heads and water contents were registered inthis experiment compared to other two runs. Consequently,
while the optimization program for run 1 could use actual
information about water contents and pressure heads in a
relatively dry range, in the other two runs that used slower
rotational speeds, the optimization program only extrapolated soil hydraulic properties in this range frommeasurement at higher-pressure heads.
[32] Figure 6 shows retention curves resulting from boththe equilibrium and transient analyses. All retention curves
are in a relatively narrow range. There are again larger
differences in obtained retention curves only close tosaturation and at the dry end. These differences are similar
to those reported in Figure 5.
[33] Finally, Figure 7 shows comparison of hydraulicconductivity functions obtained using numerical inversionof transient centrifugal experiments and independently
measured hydraulic conductivities using steady state centrifuge experiments [Nimmo et al., 1987]. Figure 7 also shows
the analytical function fitted to measured hydraulicconductivities by Nimmo et al. [1987]. Notice that all
hydraulic conductivity functions pass in a relatively narrow

range around measured values, and that they deviatesignificantly between themselves only for pressure heads
for which there are no measurements. Extrapolation of
optimized soil hydraulic functions beyond the range ofmeasurements is usually expected to be associated with ahigh level of uncertainty [e.g., Hopmans et al., 2002].

Figure 5. Optimized retention curves for three runs
compared to water contents obtained assuming equilibrium
at each time of change of rotational speed. See color versionof this figure in the HTML.

7 of 9

W04015ŠIMŮNEK AND NIMMO: TRANSIENT ANALYSES OF FLOW EXPERIMENTS IN A CENTRIFUGE

W04015

Hydraulic conductivity functions fitted to transient centrifugal data show similar shapes and values.

6. Summary[34] High-speed centrifuges have been often used during
the last two decades to measure water retention, andsaturated and unsaturated hydraulic conductivities [e.g.,
Nimmo et al., 1987; Conca and Wright, 1992] of soil androck samples. However, a steady state analysis is generally
used to compute K(q) from the measurements. Recent
advances in centrifuge experiments allow measurements ofmultiple transient variables, such as pressure heads, watercontents, concentrations, and outflow rates during theexperiment, i.e., variables that can be used to advantage
with parameter estimation techniques.
[35] In this paper we presented a modified version of
the Hydrus software package that can directly or inversely
simulate water flow in a centrifugal field, and the application of this new model to previously published data of
Nimmo [1990]. Using time-variable water contents measured at a sequence of several rotational speeds, we
successfully estimated soil hydraulic properties by numerical inversion of transient experiments. The optimized soil
hydraulic properties compared well with those determinedusing equilibrium analysis and steady state experiments.
[36] Similar parameter estimation procedures are often
used in soil physics for analyzing one-step and multistep
outflow experiments [Hopmans et al., 2002]. Multirotation
experiments in a centrifuge not only offer significant timesavings by accelerating time, but also provide significantly
more information for the parameter estimation procedure.
Multistep outflow experiments in a gravitational field, if
carried out for a sufficiently long time to reach steadyconditions, provide a single pair of the water content andpressure head at each pressure step that defines one data
point of the retention curve. Multirotation experiments, on
the other hand, can provide multiple pairs of pressure heads
and water contents for a single rotational speed. This is

Figure 7. Hydraulic conductivity functions obtained using
the parameter estimation analyses of transient centrifugal
experiments and hydraulic conductivities measured using
steady state analyses [Nimmo et al., 1987]. The thick solidline was fitted to measured data obtained from steady stateexperiments. See color version of this figure in the HTML.
because the pressure head changes significantly between thebottom and the top of the sample in a centrifugal field, while
it is almost uniform for the multistep experiment in a
gravitational field. The advantage of collecting multiple
retention data points at one speed is partly overcome byincreasing nonlinearity of water content distribution withina measurement window over which the water content isaveraged compared to the usual gravity system. In transientconditions a wide range of pressure heads also occurs for
multistep outflow experiments and these pressure heads canbe measured with tensiometers and used in the inverse
analyses. A detailed sensitivity analyses on both measured
and synthetic data, however, would then be required toreveal which experiment provides more information and ismore reliable to determine soil hydraulic properties. It
seems plausible that if multirotation experiments are carriedout for a sufficiently long time at each rotational speed toreach steady conditions, the collected data can provide
multiple points of the retention curve that would warrant the
optimization of free form analytical functions such as thosesuggested by Bitterlich et al. [2004].
[37] Acknowledgment. The senior author would like to thank Earl
Mattson and Carl Palmer from Idaho National Laboratory for providing
funding to develop the numerical solver for the centrifugal conditions.

References

Figure 6. Retention curves resulting from the equilibrium(SS) and transient (Transient) analyses. Number in thelegend refers to run number. See color version of this figurein the HTML.

Alemi, M. H., D. R. Nielsen, and J. W. Biggar (1976), Determining thehydraulic conductivity of soil cores by centrifugation, Soil Sci. Soc.
Am. J., 40, 212 – 218.
Bear, J., M. Y. Corapcioglu, and J. Balakrishna (1984), Modeling ofcentrifugal filtration in unsaturated deformable porous media, Adv. Water
Resour., 7, 150 – 167.
Bitterlich, S., W. Durner, S. C. Iden, and P. Knabner (2004), Inverse
estimation of the unsaturated soil hydraulic properties from column
outflow experiments using free-form parameterizations, Vadose Zone J.,
3, 971 – 981.
Briggs, L. J., and J. W. McLane (1907), The moisture equivalent of soils,
U.S. Bur. Soils Bull., 45, 23 pp.
Celia, M. A., E. T. Bououtas, and R. L. Zarba (1990), A general massconservative numerical solution for the unsaturated flow equation, Water
Resour. Res., 26, 1483 – 1496.
Conca, J. L., and J. Wright (1992), A new technology for direct measurements of unsaturated transport, paper presented at the Nuclear and

8 of 9

W04015ŠIMŮNEK AND NIMMO: TRANSIENT ANALYSES OF FLOW EXPERIMENTS IN A CENTRIFUGE

Hazardous Waste Management Spectrum ’92 Meeting, An. Nucl. Soc.,
Boise, Idaho.
Conca, J. L., and J. A. Wright (1998), The UFA method for rapid, direct
measurement of unsaturated transport properties in soil, sediment, androck, Aust. J. Soil Res., 36, 291 – 315.
Culligan, P. J., D. A. Barry, and J. Y. Parlange (1997), Scaling unstableinfiltration in the vadose zone, Can. Geotech. J., 34(3), 466 – 470.
Dane, J. H., and J. W. Hopmans (2002), Hanging water column, in Methodsof Soil Analysis, part 4, Physical Methods, edited by J. H. Dane and G. C.
Topp, pp. 680 – 683, Soil Sci. Soc. of Am., Madison, Wis.
Dane, J. H., and G. C. Topp (Eds.) (2002), Methods of Soil Analysis, part 1,
Physical Methods, 3rd ed., Soil Sci. Soc. of Am., Madison, Wis.
Durner, W., E. B. Schultze, and T. Zurühl (1999), State-of-the-art in inverse
modeling of inflow/outflow experiments, in Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media,
edited by M. T. van Genuchten et al., pp. 661 – 681, Univ. of Calif.,
Riverside.
Eching, S. O., and J. W. Hopmans (1993), Optimization of hydraulic functions from transient outflow and soil water pressure data, Soil Sci. Soc.
Am. J., 57, 1167 – 1175.
Forbes, P. L. (1994), Simple and accurate methods for converting centrifuge
data into drainage and imbibition capillary pressure curves, Log Anal.,
35, 31 – 53.
Gardner, R. (1937), A method of measuring the capillary tension of soil
moisture over a wide moisture range, Soil Sci., 43, 277 – 283.
Gardner, W. R. (1956), Calculation of capillary conductivity from pressure
plate outflow data, Soil Sci. Soc. Am. Proc., 20, 317 – 320.
Griffioen, J. W., P. J. Culligan, D. A. Barry, and K. Banno (1997),
Centrifuge scaling of unstable infiltration, Recent Res. Dev. Soil Sci.,
1, 29 – 41.
Hassler, G. L., and E. Brunner (1945), Measurement of capillary pressuresin small core samples, Trans. Am. Inst. Min. Metall. Pet. Eng., 160, 114 –123.
Hopmans, J. W., J. Šimůnek, N. Romano, and W. Durner (2002), Inverse
modeling of transient water flow, in Methods of Soil Analysis, part 1,
Physical Methods, edited by J. H. Dane and G. C. Topp, pp. 963 – 1008,
3rd ed., Soil Sci. Soc. of Am., Madison, Wis.
Kool, J. B., J. C. Parker, and M. T. van Genuchten (1985), 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.
Marquardt, D. W. (1963), An algorithm for least-squares estimation of
nonlinear parameters, SIAM J. Appl. Math., 11, 431 – 441.
Mattson, E. D., K. E. Baker, C. D. Palmer, R. W. Smith, and J. Šimůnek(2003), One-dimensional solute transport in variably saturated soil using
a geocentrifuge apparatus, Eos Trans. AGU, 84(46), Fall Meet. Suppl.,
Abstract H22A-0906.
Mualem, Y. (1976), A new model for predicting the hydraulic conductivityof unsaturated porous media, Water Resour. Res., 12, 513 – 522.
Nakajima, H., E. D. Mattson, and A. T. Stadler (2003), Unsaturated
hydraulic properties determined from geocentrifuge tests, Eos Trans.
AGU, 84(46), Fall Meet. Suppl., Abstract H22A-0908.
Nimmo, R. J. (1990), Experimental testing of transient unsaturated flow
theory at low water content in a centrifugal field, Water Resour. Res.,
26(9), 1951 – 1960.
Nimmo, J. R., and K. C. Akstin (1988), Hydraulic conductivity of a sandy
soil at low water content after compaction by various methods, Soil Sci.
Soc. Am. J., 52, 303 – 310.

W04015

Nimmo, R. J., and K. A. Mello (1991), Centrifugal techniques formeasuring saturated hydraulic conductivity, Water Resour. Res., 27(6),1263 – 1269.
Nimmo, J. R., J. Rubin, and D. P. Hammermeister (1987), Unsaturated flow
in a centrifugal field: Measurement of hydraulic conductivity and testingof Darcy’s law, Water Resour. Res., 23(1), 124 – 137.
Nimmo, J. R., K. S. Perkins, and A. M. Lewis (2002), Steady-state
centrifuge, in Methods of Soil Analysis, part 1, Physical Methods,
3rd ed., edited by J. H. Dane and G. C. Topp, pp. 903 – 916, Soil Sci.
Soc. of Am., Madison, Wis.
Paningbatan, E. P. (1980), Determination of soil moisture characteristic andhydraulic conductivity using a centrifuge, Ph.D. thesis, Univ. of Calif.,
Davis.
Russell, M. B., and L. A. Richards (1938), The determination of soil
moisture energy relations by centrifugation, Soil Sci. Soc. Am. Proc.,
3, 65 – 69.
Savvidou, C., and P. J. Culligan (1998), The application of centrifuge
modeling to geo-environmental problems, Proc. Inst. Civ. Eng. Geotech.
Eng., 131, 152 – 162.
Schaap, M. G., and F. J. Leij (2000), Improved prediction of unsaturated
hydraulic conductivity with the Mualem – van Genuchten model, Soil Sci.
Soc. Am. J., 64, 843 – 851.
Sigda, J. M., and J. L. Wilson (2003), Are faults preferential flow paths
through semiarid and arid vadose zones?, Water Resour. Res., 39(8),1225, doi:10.1029/2002WR001406.
Šimůnek, J., and J. W. Hopmans (2002), Parameter optimization and nonlinear fitting, in Methods of Soil Analysis, part 1, Physical Methods,
3rd ed., edited by J. H. Dane and G. C. Topp, pp. 139 – 157, Soil Sci.
Soc. of Am., Madison, Wis.
Šimůnek, J., O. Wendroth, and M. T. van Genuchten (1998a), A parameter
estimation analysis of the evaporation method for determining soil
hydraulic properties, Soil Sci. Soc. Am. J., 62, 894 – 905.
Šimůnek, J., M. Šejna, and M. T. van Genuchten (1998b), The HYDRUS1D software package for simulating the one-dimensional movement ofwater, heat, and multiple solutes in variably-saturated media, version 2.0,
Rep. IGWMC-TPS-70, 202 pp., Int. Ground Water Model. Cent., Colo.
Sch. of Mines, Golden, Colo.
van Dam, J. C., J. N. M. Stricker, and P. Droogers (1994), Inverse method
to determine soil hydraulic functions from multistep outflow experiment,
Soil Sci. Soc. Am. J., 58, 647 – 652.
van Genuchten, M. T. (1980), A closed-form equation for predicting thehydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44,892 – 898.
van Genuchten, M. T., F. J. Leij, and S. R. Yates (1991), The RETC codefor quantifying the hydraulic functions of unsaturated soils, Rep. EPA/
600/2-91-065, U.S. Enviorn. Prot. Agency, Washington, D. C.
Wendroth, O., W. Ehlers, J. W. Hopmans, H. Kage, J. Halbertsma, andJ. H. M. Wösten (1993), Reevaluation of the evaporation method for
determining hydraulic functions in unsaturated soils, Soil Sci. Soc. Am.
J., 57, 1436 – 1443.




J. R. Nimmo, U.S. Geological Survey, 345 Middlefield Road, MS-421,
Menlo Park, CA 94025, USA. (jrnimmo@usgs.gov)
J. Šimůnek, Department of Environmental Sciences, University ofCalifornia Riverside, Riverside, CA 92521, USA. (jiri.simunek@ucr.edu)

9 of 9