European Journal of Soil Science, March 2001, 52, 13±24

Non-equilibrium water ¯ow characterized by means of
upward in®ltration experiments
J . SÏ I M UÊ N E K a , O . W E N D R O T H b , N . W Y P L E R b & M . T . V A N G E N U C H T E N aaGeorge E. Brown Jr Salinity Laboratory, USDA-ARS, Riverside, CA 92507, USA, and bInstitute for Soil Landscape Research, ZALF,
15374 MuÈncheberg, Germany

Summary
Upward in®ltration experiments under tension were used to demonstrate the presence of non-equilibrium¯ow in soils, the phenomenon that has important implications for the accelerated movement of fertilizers,
pesticides, non-aqueous liquids, and other pollutants. Data obtained from these experiments wereanalysed using the single-porosity Richards equation, as well as a variably saturated, dual-porosity modeland a dual-permeability model for characterizing non-equilibrium water ¯ow. The laboratory
experiments were carried out on 0.10-m-long soil cores having an internal diameter of 0.10 m.
Constant pressure heads of ±0.10 and ±0.01 m were used as the lower boundary condition. Eachin®ltration was followed by a single-rate evaporation experiment to re-establish initial conditions, and toobtain the drying soil hydraulic properties. Pressure heads inside the cores were measured using ®vetensiometers, while evaporative water loss from the top was determined by weighing the soil samples.
The data were analysed to estimate parameters using a technique that combined a numerical solution of
the governing ¯ow equation (as implemented in a modi®ed version of the Hydrus-1D software) with aMarquardt±Levenberg optimization. The objective function for the parameter estimation was de®ned interms of pressure head readings, the cumulative in®ltration rate, and the ®nal total water volume in thecore during upward in®ltration. The ®nal total water volume was used, as well as the pressure headreadings during the evaporation part. Analysis of ¯ow responses obtained during the in®ltration
experiment demonstrated signi®cant non-equilibrium ¯ow. This behaviour could be well characterizedusing a model of physical non-equilibrium that divides the medium into inter- and intra-aggregate poreswith ®rst-order transfer of water between the two systems. The analysis also demonstrated the importanceof hysteresis.

Introduction
Local soil heterogeneity often results in the preferential
movement of water and chemicals through macropores.
Whereas in unsaturated fractured rock water moves preferentially through fractures and ®ssures and bypasses much of therock matrix, in aggregated and macroporous soils water moves
both through pores between aggregates, decayed root channels,
and earthworm burrows and through the soil matrix itself. An
important characteristic of preferential ¯ow is that during
wetting the moisture front can propagate quickly to signi®cant
depths while bypassing a large part of the matrix pore space.
Non-equilibrium preferential movement has received muchattention by soil scientists since it has important implicationsin terms of the accelerated movement of fertilizers, pesticides,
non-aqueous liquids, and other pollutants applied to the soil's
surface. Water and solutes move to far greater depths, andCorrespondence: J. SÏimuÊnek. E-mail: jsimunek@ussl.ars.usda.gov
Received 17 April 2000; revised version accepted 8 September 2000
# 2001 Blackwell Science Ltd

much faster, than would be predicted using the Richards
equation based on area-averaged moisture contents andpressure heads (Beven, 1991).
Preferential ¯ow in structured media is usually describedusing dual-porosity or dual-permeability models (Pruess &Wang, 1987; Gerke & van Genuchten, 1993; Jarvis, 1999). In
both the porous medium is assumed to consist of two
overlapping interacting regions, one associated with theinter-aggregate, macropore, or fracture system, and the otherwith less permeable micropores inside soil aggregates (intraaggregates), or the rock matrix. While in dual-porosity models
the water in the matrix is stagnant, dual-permeability models
allow for water to ¯ow in the matrix as well. Models of thistype have long been applied to solute transport studies.
Especially popular have been dual-porosity models in which
distinct mobile and immobile liquid ¯ow regions are assumed(e.g. van Genuchten & Wierenga, 1976), although dualpermeability models in which water can move in both theinter- and intra-aggregate pores are also becoming more
13

14 J. SÏimuÊnek et al.
popular (Pruess & Wang, 1987; Gerke & van Genuchten,1993).
The main disadvantage of dual-porosity or dual-permeability models is that, contrary to models based on a single
pore region, they require many more input parameters tocharacterize both pore systems. Little guidance is available as
to how to obtain these parameters, either by direct measurement, a priori estimation, or some calibration technique(Beven, 1991; Clothier et al., 1995; Jaynes et al., 1995).
Hence, we need to design experiments or devices that provide(either by inverse modelling or directly) parameters for these
relatively complex models. Solute transport parameters ofdual-porosity models are often obtained from column experiments on the assumption that water ¯ow is steady (NkediKizza et al., 1984), with codes for parameter estimation suchas CXTFIT (Toride et al., 1995) and/or STANMOD (SÏimuÊnek
et al., 1999b) that ®t analytical solutions of the transport
equation to experimental breakthrough curves. More sophisticated experiments and models are needed also to characterizeparameters when ¯ow is unsteady and the soil is unsaturated.
These parameters cannot usually be obtained by direct
measurement or independent calibration.
Methods for optimizing parameters are well established for
estimating hydraulic properties of single-porosity systems(Kool et al., 1987). They make it possible to estimate
simultaneously the retention and hydraulic conductivity
functions from data on transient ¯ow (Kool et al., 1987).
While they may be applied in many scenarios, numericalinversion of the Richards equation has thus far been limited
mostly to one-dimensional experiments in the laboratory (Koolet al., 1985; Russo et al., 1991; Hudson et al., 1996), usually bymeans of one-step (Kool et al., 1985) or multistep (van Dam
et al., 1994) out¯ow experiments, but increasingly also using
evaporation experiments (Santini et al., 1995; SÏimuÊnek et al.,1998b).
Experiments on out¯ow and evaporation both represent
desorption and thus lead to parameters associated with drying
branches of the soil's hydraulic functions. Often nonequilibrium ¯ow is associated with wetting and not drying,
which means that experiments designed to study nonequilibrium ¯ow should involve in®ltration. Several applications of inverse analyses to in®ltration have been reported,
including a hypothetical one-dimensional ponded in®ltration
experiment (Russo et al., 1991), axi-symmetrical in®ltrationfrom a tension disc permeameter (SÏimuÊnek et al., 1999c), andin®ltration from a modi®ed cone permeameter (Gribb et al.,
1998; SÏimuÊnek et al., 1999a). Of these, only SÏimuÊnek et al.
(1999c), in their experiment on in®ltration from a tension disc,observed signi®cant non-equilibrium ¯ow.
Hudson et al. (1996) suggested an upward in®ltration
experiment for laboratory conditions. They imposed a constant¯ux of water at the bottom of the soil sample, and measured
pressure heads inside the sample using tensiometers. To
maximize information for the inverse analysis, we suggest here
# 2001 Blackwell Science Ltd, European Journal of Soil Science, 52, 13±24

to initiate in®ltration by a certain tension at the bottom of thesample, rather than by imposing a boundary ¯ux. Since for ¯ux
in®ltration the precise amount of water in the sample is always
known, the only information characterizing the soil's hydraulicproperties is the shape of the wetting front as measured withtensiometers. For tension in®ltration, the soil also controls thetotal amount of water being taken up, thus providing additional
information and control for the numerical inversion. We
believe that the way water ¯ow responds during transientin®ltration can be used to calibrate models for aggregatedsoils. However, the question remains as to whether such data
are enough to identify fairly complex models of ¯ow required
to describe the non-equilibrium ¯ow.
The objective of this study was to identify both thedrying and wetting branches of the soil's hydraulicproperties from consecutive experiments in which thepressure head controlled upward in®ltration and evaporation
followed. We also examined the validity or reproducibilityof hydraulic properties inversely estimated in this mannerby doing the in®ltration experiments at the same anddifferent boundary conditions of pressure head, and by
several evaporation experiments. In addition, we used anexperiment on upward in®ltration under tension to demonstrate the presence of non-equilibrium ¯ow in the soil. The
data were analysed using both equilibrium and nonequilibrium ¯ow models. Non-equilibrium behaviour wascharacterized using a model of physical non-equilibriumwith two regions that divided the pore system into interand intra-aggregate pores with ®rst-order transfer of water
between the two pore systems. Both dual-porosity and dualpermeability non-equilibrium models were used to analyse
the data.

Experimental
Soil
Cylindrical undisturbed soil samples (0.10 m high and 0.10 m
internal diameter) were taken from the 0.05±0.15 m depth
layer of the Ap horizon of a sandy loam ®eld site in Lietzen,
Brandenburg, Germany. The soil was under conservation
tillage, and had a fairly large dry bulk density, b, of 1.65 g
cm3, which corresponds to an air-®lled porosity of about 0.377(= 1 ± b/2.65). The soil consisted of 62% sand, 32% silt, and
6% clay. The organic carbon content was 0.75%. The structureof this soil was very compact, typical for sandy soils that are
intensively tilled and compacted afterwards by traf®c.
Although aggregates could not be identi®ed, the soil matrix
is heterogeneous with respect to bulk density, as was shown bycomputer tomography. Apparently, ploughing often breaks
compacted zones, thus causing development of a fragmented
matrix with locally dense regions that alternate with looser
ones (H. Rogasik, personal communication).

Non-equilibrium water ¯ow
Table 1 Summary of experimental conditions

Experiment

Duration/days

Average®nal 
/m3m±3

Initial conditionIn®ltration I, h0 = ±0.10 m
Redistribution I
Evaporation I
Redistribution II
In®ltration II, h0 = ±0.01 m
Redistribution III
Evaporation II
Redistribution IV
In®ltration III, h0 = ±0.01 m

NA
2.54382.41535.15739.97130.96931.90774.40281.86380.9743

0.068
0.200
0.199
0.107
0.103
0.231
0.229
0.104
0.106
0.256

In®ltrated
volume/mm

15

quality that change much more smoothly than pressure heads,
which re¯ect the movement of moisture fronts. Measuring thepressure heads at longer intervals could quickly lead tomissing information about advancing moisture fronts in thesample, but this is less probable when measuring sample
weights.

13.12

Evaporation experiments
12.65

14.83

In®ltration experiments
Prior to the experiments, ®ve pressure transducer tensiometers
were installed at 10, 30, 50, 70, and 90 mm depth in the soil
(henceforth denoted as locations 1 to 5). The tensiometers,
consisting of ceramic cells with an inner diameter of 6 mm,
were inserted horizontally into the soil across a length of
50 mm. Good contact between the soil and tensiometers was
ensured by using a paste from a fraction of the soil drilled outof the sample before installation. The reaction time of ourtensiometers to a 10% change of pressure head was about 3 s at
a pressure h = ±0.1 m and up to 1 minute for h = ±6.0 m. The
samples during upward in®ltration were positioned on a porous
membrane that was used as a suction plate. At the intended
pressure head, the surface of the membrane was in hydraulicequilibrium with a Mariotte reservoir bottle that was placed ona balance. Before and after the in®ltration, the soil samples
were weighed to determine total amount of water in thesamples at the start and at the end. During the ®rst 30 h of the®rst upward in®ltration, weight losses from the reservoir due towater uptake of the soil sample were recorded every 5 minutes
using a video camera. Afterwards, when the weight losses
became smaller, they were recorded manually. Tensiometerreadings were recorded automatically every 2 s. After eachin®ltration experiment, the sample was sealed and water wasallowed to redistribute before starting the evaporation.
Additional water was allowed to in®ltrate into the soil sample
during the redistribution process for one experiment.
Three upward in®ltration experiments (I, II, and III) were
made, with lower boundary conditions of h0 = ±0.10, ±0.01,
and ±0.01 m, respectively. Each in®ltration was followed by anevaporation stage, and subsequently by a period of equilibration of the soil sample to re-establish similar initial conditions.
Weights of the reservoir bottle during in®ltration experiments
II and III were recorded with a data logger every 30 s. The
larger intervals for weighting the reservoir bottle duringin®ltration, as compared with measuring pressure heads, wasbased on the assumption that weights represent a cumulative

Evaporation experiments similar to those by Wendroth et al.
(1993) and SÏimuÊnek et al. (1998b) were done using the same
soil column as for the upward in®ltration experiments. Onemodi®cation from the earlier experiments was that evaporation
was allowed to proceed slowly, as dictated by the laboratory
conditions. The sample was weighed before and after theexperiments. During evaporation, tensiometers were recorded
automatically every 5 minutes, and the sample weighed every
4 h, except for two intervals during the night.
Prior to each evaporation experiment, after each in®ltration,
when pressure head changes ceased, hydraulic equilibrium wasassumed, and deviations in pressure head from equilibrium
were corrected in the calibration functions of the pressuretransducers of the ®ve tensiometers by the method ofWendroth et al. (1993). Moreover, the calibration of thepressure transducer was divided for each speci®c transducer
into two segments, one between 0 and 1.50 m pressure head,
and the other between 1.50 and 6.50 m pressure.
Table 1 summarizes the experimental conditions, including
duration of particular parts of each experiment, average watercontents in the sample at the beginning and end of eachexperiment, and the in®ltrated volumes. The ®rst in®ltration
experiment with a supply pressure of ±0.10 m lasted 2.5 days.
The second and third experiments with applied heads of±0.01 m lasted about 1 day. Evaporation experiments I and II
lasted about 5.1 and 4.4 days, respectively. Again, in®ltrationand evaporation experiments were separated by equilibration
periods of between 2 and 10 days. During these both thebottom and the top of the sample were sealed to prevent water
loss, thereby allowing water to redistribute inside the sample.

Methods
Below we give the governing equations for three models
describing unsaturated water ¯ow in soils. We ®rst discuss theequilibrium ¯ow of water using the standard single-porosity
model, after which we give the governing equations for themodels of dual-porosity and dual-permeability non-equilibrium ¯ow. Models of dual-porosity and dual-permeability
require many more input parameters than ones of singleporosity. For example, dual-porosity models may require
retention parameters for both the mobile and immobile
regions, hydraulic conductivity parameters for the mobile
region, and parameters of mass transfer between regions. Dualpermeability models require even more parameters since they
# 2001 Blackwell Science Ltd, European Journal of Soil Science, 52, 13±24

16 J. SÏimuÊnek et al.
must also account for ¯ow in the matrix. They need twofunctions for water retention, one for the matrix and one for the
pores between the aggregates. Also two or three hydraulicconductivity functions are required: one for the inter-aggregate
pores, another for the matrix, and sometimes yet another for
the interface between the aggregates and the matrix (Gerke &van Genuchten, 1993). We discuss below several options that
could minimize the number of parameters needed to describe
particular non-equilibrium models.
Equilibrium ¯ow model
Water ¯ow under equilibrium conditions may be describedusing a mixed formulation of the Richards equation:



@
@
@h

K h
1 ;
@t @x
@x

1

where x is the vertical coordinate (L), taken here to be positive
upward, t is time (T),  is the water content (L3L±3), h is
the pressure head (L), and K is the hydraulic conductivity(LT±1).
Soil hydraulic properties
We describe the soil hydraulic properties using the van
Genuchten±Mualem model (van Genuchten, 1980):
Se h  h ÿ r1s ÿ r1 j h jn m

2

and
m 2
K   Ks Sle 1 ÿ 1 ÿ S1=me   ;

3

where Se is effective ¯uid saturation (dimensionless), Ks is thesaturated hydraulic conductivity (LT±1), r and s denote theresidual and saturated water contents (L3L±3), respectively; l isa pore-connectivity parameter (dimensionless), and (l±1), andn and m (= 1 ± 1/n) are empirical shape parameters. The above
hydraulic functions contain six unknown, independent parameters: r, s, , n, l and Ks.




@m
@
@h

K h
 1 ÿ ÿ w;
@t
@x
@x

where m is the inter-aggregate water content (dimensionless)and ÿ w is the transfer rate for water from the inter- to the intraaggregate pores (T±1). This rate is described using the ®rstorder rate equation
ÿw # 2001 Blackwell Science Ltd, European Journal of Soil Science, 52, 13±24@im
im
 !Sme ÿ Se ;@t5

where im is the intra-aggregate water content, ! is a ®rst-order
rate coef®cient (T±1), and Sem and Seim are effective ¯uid
saturations in the inter- and intra-aggregate regions, respectively.
To minimize the number of parameters in the model, and to
distinguish our dual-porosity model from the more comprehensive dual-permeability model described below, we assumed
that mass transfer between the inter- and intra-aggregate
regions is proportional to the difference in their effective watercontents, rather than to their pressure heads (Gerke & van
Genuchten, 1993). By using Equation (5) we decreased thenumber of soil hydraulic parameters, since we need not know
the retention and hydraulic conductivity functions for the intraaggregate region explicitly. Thus, the dual-porosity nonequilibrium model contains nine parameters: the same six
parameters (r, s, , n, l and Ks) as for the equilibrium model(but now pertaining to the inter-aggregate region), two
additional factors characterizing the intra-aggregate region(i.e. its residual, rim, and saturated, sim, water contents), andthe ®rst-order mass transfer coef®cient !. By further assumingthat the residual water content of the inter-aggregate region isequal to zero (hence that residual water is present only in theimmobile region), we can decrease the number of model
parameters further to eight.
Dual-permeability non-equilibrium ¯ow model
In the dual-permeability model Richards' equations can beapplied to each of the two pore regions (Gerke & van
Genuchten, 1993, 1996). The ¯ow equations for the inter(subscript m) and intra-aggregate (subscript im) pore systems
are, respectively,

Dual-porosity non-equilibrium ¯ow model
In dual-porosity models the water in the pores within theaggregates does not move vertically; only that in the pores
between the aggregates ¯ows vertically. Thus, intraaggregate pores represent immobile pockets that can
exchange, retain and store water, but are not subject toconvective vertical ¯ow. The dual-porosity formulation usedhere is based on a mixed formulation of the Richards
equation to describe water ¯ow in the pores between theaggregates:

4

@m
@
@hm
ÿw

 Km  ÿ
Km
@t
@z
wm
@z

6

@im
@
@him
ÿw

 Kim  
;
Kim
@t
@z
1 ÿ wm@z

7

andwhere wm is equal to the ratio of the inter-aggregate and total
water contents, m/, and the rate of exchange of waterNon-equilibrium water ¯ow
between the inter- and intra-aggregate regions, ÿ w (T±1), is
described asÿw 

w hm ÿ him ;8

in which w is a ®rst-order mass transfer coef®cient (L±1T±1).
Unlike for the dual-porosity model, the mass transfer, ÿ , inEquation (8) is assumed to be proportional to the difference inpressure heads between the two pore regions. This approachrequires estimating retention curves for both pore regions. For
porous media with well-de®ned geometries, the ®rst-ordermass transfer coef®cient, w, can be de®ned as follows (Gerke& van Genuchten, 1996):w a2

Ka w ;

9

where is a shape factor that depends on the geometry. The
value ofranges from 3 for rectangular slabs to 15 forspherical aggregates. In Equation (9), a is the characteristic
length of the aggregate (L) (e.g. the radius of a spherical orsolid cylindrical aggregate, or half the width of a rectangularaggregate), and w (= 0.4) is a dimensionless scaling factor.
The effective hydraulic conductivity Ka (LT±1) of the interface
between a fracture and the matrix can be simply evaluated inarithmetic terms of both hm and him:
Ka(h) = 0.5{Ka(hm) + Ka(him)}.

(10)

The use of Equation (9) implies that the medium contains
geometrically well-de®ned cylindrical, rectangular or othertypes of macropores or fractures; see, for example, Edwards
et al. (1979) and van Genuchten & Dalton (1986). While
geometrically based models are conceptually attractive, they
may be too complicated for routine applications because
aggregated soils and rocks usually contain mixtures ofaggregates of various sizes and shapes. Hence, rather thanusing Equation (9) directly, we lump , a, and w into thede®nition of the effective hydraulic conductivity Ka* at theinterface to givew  Ka h:11

The above dual-permeability model has, in its fullcomplexity, many parameters. The porous regions between
and within the aggregates are each characterized by thesame six soil hydraulic parameters as the equilibrium
system. Additionally, estimates of the parameters Ksa*, a,and na are needed to characterize the hydraulic conductivity
at the interface when using the Mualem±van Genuchten
model. A previous assumption (Gerke & van Genuchten,
1993, 1996) is that the relative conductivity functions of
the interface and the intra-aggregate region are the same,requiring thus only the coef®cient Ksa* to scale the relative
conductivity function. If the pore connectivity parameters, l,

17

are equal to 0.5 (Mualem, 1976) and, as before, theresidual water content of the mobile region is zero, thenumber of parameters decreases to 10: sm, rim, sim, m,mm
and Ksa*. Alternatively, the two
im, nim, nim, Ks , Ks
parameters s and wm could be used instead of sm andsim.
Initial and boundary conditions
The lower and upper boundary conditions for an upwardin®ltration induced by tension are as follows:h  h0x  0; t > 0



@h
1 0
ÿK h
@x

x  L; t > 0;

and
12

where h0 is the pressure head applied at the bottom of thesample, and L is the height of the sample. The boundary
conditions for the evaporation experiment are


@h
1 0
ÿK h
@x

x  0; t > 0

and
ÿK h



@h
 1  qevap t
@xx  L; t > 0;13

where qevap is the evaporation rate, and for the equilibration period between the in®ltration and evaporation
experiments


@h
1 0
ÿK h
@x

x  0; t > 0



@h
1 0
ÿK h
@x

x  L; t > 0:

and
14

The initial condition for all three experimental phases isgiven in terms of pressure heads, as measured with thetensiometers:h  hi x

x > 0; t  0:15Equations (1), (4), or (6) and (7), subject to initial andboundary conditions speci®ed in Equations (12) to (15), weresolved using a modi®ed version of the Hydrus-1D numericalcode (SÏimuÊnek et al., 1998a). Modi®cations involved imple# 2001 Blackwell Science Ltd, European Journal of Soil Science, 52, 13±24

18 J. SÏimuÊnek et al.

F ; qm  

m
X
j1

"
vj

nj
X
i1

#
wij fqj ti  ÿ qj ti ;

2

g ;

16

where m represents different sets of measurements; nj is thenumber of measurements in a particular set, qj*(ti) is the
speci®c measurement at time ti for the jth measurement set,is the vector of optimized parameters (e.g. r, s, , n, Ks, and lfor equilibrium problems, and rim, sm, sim, m, m, nm, nim,
Ksm, Ksim, and Ksa* for non-equilibrium problems), qj(ti, )represents the corresponding model predictions for parameter
vector , and vj and wij are weights associated with a particularmeasurement set j or a measurement i within set j,
respectively. We assume that the weighting coef®cients wij
in Equation (16) are equal to one, that is, the variances of theerrors inside a particular measurement set are all the same. The
weighting coef®cients vj are given by
vj 

1
:
nj 2j17

This approach views the objective function as the average
weighted squared deviation that is normalized by measurement
variances j2.

Data analysis
Analysis of upward in®ltration experiments

Figure 1 Measured pressure heads and cumulative in®ltration rates
for the three in®ltration experiments (I, II, and III) subject toapplied pressure heads of ±0.10 (a), ±0.01 (b), and ±0.01 (c) m,
respectively.

menting the models of dual-porosity and dual-permeability¯ow.
Parameter estimation
We minimized the objective (or merit) function that measures
the agreement between measured and modelled data using theLevenberg±Marquardt non-linear method (Marquardt, 1963).
The objective function was de®ned by means of weighted
least-squares estimator as follows:# 2001 Blackwell Science Ltd, European Journal of Soil Science, 52, 13±24

Figure 1 shows measured pressure heads and cumulative
upward in®ltration rates for the three in®ltration experiments
subject to applied pressure heads of ±0.10 (top), ±0.01(middle), and ±0.01 (bottom) m. The upper two graphs alsopresent the ensuing redistribution phase of the experiment.
Notice that different time scales are used in each ®gure. The
cumulative in®ltration curves have a typical shape during theearly phase of the in®ltration until all ®ve tensiometers registera fairly steady pressure head corresponding to the equilibrium
pro®le. At that moment, according to theory, in®ltration shouldstop, so that the in®ltration rate becomes zero. However, we
continued monitoring the water intake by the soil sample andnoticed a continuing and almost constant in®ltration rate ineach of the three experiments. The rates stabilized at values ofabout 2, 8, and 5 mm day±1 for experiments I, II, and III,
respectively (Figure 2). Tensiometer readings during thissecond stage of the experiment did not change signi®cantly.
Although tensiometers integrate pressure heads along thesurface of the tensiometer cup, the readings will have a strong
tendency to be close to the largest pressure (closest to zero) in
contact with the surface of the tensiometer cup. The continuingin®ltration without corresponding changes in the pressure
heads suggests that water is redistributing within the soil
sample during this second phase of the in®ltration. A likely
explanation for this redistribution is that water is beingNon-equilibrium water ¯ow

19

Figure 3 Measured pressure heads during the redistribution process.

Figure 2 Measured in®ltration rates for the three in®ltration
experiments (I, II, and III) with applied pressure heads of ±0.10 (a),±0.01 (b), and ±0.01 (c) cm, respectively.

transferred from the larger (inter-aggregate) to the smaller(intra-aggregate) pores.
The sample was sealed at both sides at the ends of thein®ltration experiments, and water was allowed to redistributewithin the sample. Since in®ltration occurred from the bottom,
the ®nal pressure head pro®le should correspond to theequilibrium state, and no changes, in pressure head readings
should be registered afterwards. In our ®rst two experiments,
however, the pressure head decreased (suctions increased) inboth cases by about 0.22 m throughout the soil sample (Figure

1a,b). The smaller pressure heads and accompanying smaller
water contents again suggest slow transfer of water from theinter-aggregate to the intra-aggregate pores.
The redistribution of water between the two pore regions
was further studied in detail on a second sample that wassealed also at both sides at the end of the in®ltration
experiment. In this case we allowed water to redistributewithin the sample for about 10 days, while small amounts ofwater were added several times during redistribution to the topof the sample. Figure 3 shows the observed tensiometerreadings. Notice that the pressure heads steadily decreased
during the ®rst 5 days of redistribution until they reachedalmost ±0.50 m at the top tensiometer. At that time about
0.8 mm of water was added to the top of the sample, whichcaused the pressure heads to increase almost immediately toabout ±0.10 m. Similar increases occurred after the otheradditions of water. We again think that the decrease in pressure
is caused by redistribution of water between different poreregions within the soil sample, with the tensiometers
registering predominantly pressure heads of the larger andinitially wetter pores.
One can estimate the average mobile water content from theadvancing moisture front, assuming that there is no transfer ofwater between the two pore regions. Table 2 presents theaverage mobile water contents for all three in®ltration
experiments calculated for each depth where a tensiometer
was installed. Since for these calculations we assume thatwater ¯ows only in the inter-aggregate region, without
interaction with the intra-aggregate region, we calculated theaverage mobile water contents by simply dividing the
cumulative in®ltration at a time when the centre of the
moisture front reached a particular tensiometer by the height of
the tensiometer above the bottom of the sample. We assumed
that the centre of the water content front corresponded with thevalue of the pressure head equal to the arithmetic average ofthe initial and ®nal head. Note from Table 2 that, based onthese calculations, only about 6 volume percents participatedin the in®ltration. This value is between one third and one# 2001 Blackwell Science Ltd, European Journal of Soil Science, 52, 13±24

20 J. SÏimuÊnek et al.
Table 2 Mobile water contents, im, calculated from cumulativein®ltration and tensiometer readings
Experiment
Tensiometer /Depth

I

II

III

1 /10 mm
2 /30 mm
3 /50 mm
4 /70 mm
5 /90 mmAverage imOverall average im
Average ®nal 

0.066
0.052
0.058
0.076
0.066
0.0636

0.065
0.050
0.054
0.061
0.042
0.05440.0602
0.231

0.101
0.060
0.058
0.059
0.038
0.0632

0.200

0.256

fourth of the ®nal water content at the end of the in®ltration
experiments. Extending this value over the entire sample(height of 0.10 m) shows that only 6 mm of water was neededto `saturate' the sample fully. This is less than half of the waterthat actually in®ltrated for each in®ltration experiment (13.1,
12.7, and 14.8 mm; see Table 1). These numbers clearly
indicate that water had to redistribute within the sample fromthe larger (wetter) inter-aggregate pores (whose pressure heads
are measured mostly by the tensiometers) to the smaller intraaggregate pores within the soil matrix.
The ®nal water contents (Table 1) for all threein®ltration experiments were signi®cantly less than theestimated air-®lled porosity of 0.377. Signi®cant parts of
the intra-aggregate pores must therefore have been left
unsaturated at the end of the in®ltration experiments, withpossibly signi®cant amounts of entrapped air in the intraaggregate pores. These pores were being gradually ®lledwith water during redistribution (equilibration), a processthat must have desaturated the larger pores. It is also
generally accepted that the ®eld-saturated (or satiated)water content is much smaller than the porosity becauseof entrapped and dissolved air. Natural saturation in thelaboratory is often only about 0.8±0.9 of the porosity(Klute, 1986), although much smaller values have beenreported in the literature.
Numerical analysis of upward in®ltration assuming
equilibrium ¯ow
First, we optimized the soil hydraulic parameters for all threein®ltration experiments using the assumption of equilibrium
water ¯ow, i.e. Equation (1). Measured and optimized pressure
heads and cumulative in®ltrations for the ®rst two upwardin®ltration experiments with applied pressure heads of ±0.10
and ±0.01 m are presented in Figure 4. The third in®ltration
experiment produced very similar results. Optimized parameters for the three in®ltration experiments are given in
Table 3.
# 2001 Blackwell Science Ltd, European Journal of Soil Science, 52, 13±24

Figure 4 Measured (symbols) and optimized (lines) pressure heads
and cumulative in®ltration rates for the ®rst (a) and second (b)upward in®ltration experiments having applied pressure heads of±0.10 and ±0.01 cm, respectively. Optimized values were calculatedusing the assumption of equilibrium water ¯ow.

Reasonable correspondence was obtained between measured and calculated pressure heads (Figure 4). As was
expected in view of the discussions above, the equilibrium
model could not reproduce the cumulative in®ltration curves
after all tensiometers essentially showed equilibrium with thebottom boundary condition (pressure heads of either ±0.10 or±0.01 m). Whereas the measured in®ltration rates remained
constant for an extended time, the numerical model predicted,
as expected, in®ltration rates that quickly decreased to zero
once the tensiometer readings reached equilibrium. Since thisobserved behaviour could not be described using theequilibrium model, non-equilibrium models were used asshown below.
Numerical analysis of evaporation experiments
The two evaporation experiments carried out as part of thisstudy were analysed in terms of the equilibrium water ¯owNon-equilibrium water ¯ow

21

Table 3 Summary of optimization results
Experiment

Model



r

s

/m±1

n

Ks
/m day±1lEvaporation I
Evaporation II
In®ltration I
In®ltration II
In®ltration III

Equilibrium
Equilibrium
Equilibrium
Equilibrium
Equilibrium

0.002390.002731
0.09166
0.028240.02586

0.07840.08570.00450.09240.1276

0.202
0.228
0.211
0.234
0.267

1.07
1.20
2.82
2.97
2.39

2.09
2.26
1.24
1.41
1.88

0.03320.07100.02150.02160.0165

0.09360.0753±1.59±2.09±0.0003

Experiment

Model



r

s

/m±1

n

Ks
/m day±1lIn®ltration I
In®ltration II
In®ltration III

Dual-porosity
Dual-porosity
Dual-porosity

0.02903
0.03234
0.026180.000100.000180.00337

0.239
0.306
0.416

4.13
3.94
2.60

1.43
1.37
1.61

0.00803 ±1.03
0.0296 ±2.55
0.0253 ±0.00695

Experiment

Model



s

/m±1

nim

Ksim
/m day±10.05490.05910.03930.025880.283
0.377*
0.377*
0.377*

0.629
1.37
0.651
0.316

1.46
1.38
1.56
1.32

0.149e±70.1e±90.377e±70.149e±4In®ltration I
In®ltration I
In®ltration II
In®ltration III

Dual-permeability
Dual-permeability
Dual-permeability
Dual-permeability

im

m

/m±1
3.09
6.49
9.45
3.64

w

!
/day±1

0.300
0.457
0.378

0.160
0.09180.0616nm

Ksm/m day±1

w

Ksa*/m±1 day±1

1.63
1.72
1.77
1.64

0.04100.0815
0.441
0.113

0.240
0.267
0.208
0.334

0.123e±40.178e±40.694e±50.164e±4*Kept constant during optimization.

model, Equation (1), using the estimation approach suggestedby SÏimuÊnek et al. (1998b). Measured and optimized pressure
heads for both experiments are presented in Figure 5. An
excellent correspondence is seen between measured andoptimized values. The soil's hydraulic parameters estimatedfor the two evaporation experiments are also given in
Table 3. The two sets of parameters are very similar withthe exception of the parameter s, mainly a result of the greater
initial water content in the second evaporation experiment. The
soil hydraulic functions (both the retention and unsaturated
hydraulic conductivity functions) are also almost identical,
except for a region close to saturation (Figure 6). This region,
however, was beyond the measurement range since bothexperiments started, because of internal redistribution of water
within the sample after the end of the in®ltration experiment,with initial pressure heads of about ±0.40 m. Nevertheless,Figure 6 shows good consistency and repeatability of results
obtained with the two evaporation experiments.
Numerical analysis of upward in®ltration assuming nonequilibrium ¯ow
All three upward in®ltration experiments were next analysedusing both models of non-equilibrium water ¯ow. The
optimized hydraulic parameters for all inversion runs are
given in Table 3. Measured and ®tted pressure heads andcumulative in®ltration volumes using the dual-permeability
model are presented in Figure 7, again only for the ®rst two

in®ltration experiments. Results with the dual-porosity model
were very similar. Calculated pressure heads in Figure 7
correspond reasonably well with the measured values, not only
during in®ltration, but also during the equilibration or
redistribution phase. Notice in particular that the nonequilibrium model accurately described the decrease in thepressure head for all tensiometers during redistribution, when,
as we think, water moved from the inter-aggregate to the intraaggregate pores. The non-equilibrium model also correctly
reproduced the continued in®ltration when tensiometers
already registered constant pressure heads.
For our calculations of non-equilibrium ¯ow we ®xed thevalue of the total water content at saturation equal to theporosity. The fraction of inter-aggregate pores, w, wasestimated to be between 0.3 and 0.45 when using the dualporosity model, and between 0.24 and 0.33 when using thedual-permeability model. The values were slightly smaller forthe dual-permeability model since the intra-aggregate poresalso contribute somewhat to vertical water ¯ow.
Although the con®dence in the optimized parameters for the
non-equilibrium models is fairly small (because of their large
number), the parameters do re¯ect our conceptual picture ofthe ¯ow process. For example, the values of for the intraaggregate pore system are about one ®fth of those for thesystem of pores between aggregates. Small
values are
typical for ®ne-textured soils. Also, the saturated hydraulicconductivity of the inter-aggregate system was found to be 3 to7 orders of magnitude larger than Ks for the intra-aggregate# 2001 Blackwell Science Ltd, European Journal of Soil Science, 52, 13±24

22 J. SÏimuÊnek et al.

Figure 6 Soil hydraulic properties obtained by analysing datacollected with the evaporation experiments.

Figure 5 Measured and optimized pressure heads for the ®rst (a)and second (b) evaporation experiments.

system and 8 to 10 orders of magnitude larger than Ks for themass exchange term. Soil hydraulic functions for the interaggregate pore region optimized using the dual-permeability
model are presented in Figure 8. Notice that the retention andhydraulic conductivity functions have similar shapes, but withlarger differences in the estimated values of the saturated watercontent, sm, of the inter-aggregate pore region.

Discussion
The dual-porosity and dual-permeability models require large
numbers of quasi-empirical parameters that, in general, mustbe ®tted to experimental data. While both models are likely
improvements over existing equilibrium ¯ow models, theyremain simpli®cations of reality with parameters that may be
dif®cult to determine independently. Also, one can hardly
expect that complex aggregated pore systems be composed ofonly two classes of pores, one containing mobile water and one
immobile (or less mobile) water. These problems make it
dif®cult to identify parameters in dual-media models by meansother than ®tting to experimental data.
# 2001 Blackwell Science Ltd, European Journal of Soil Science, 52, 13±24

Figure 7 Measured (symbols) and optimized (lines) pressure heads
and cumulative in®ltration rates for the ®rst (a) and second (b)upward in®ltration experiments having applied pressure heads of±0.10 and ±0.01 cm, respectively. Optimized values were calculatedusing the assumption of non-equilibrium water ¯ow (dualpermeability model).

Non-equilibrium water ¯ow

23

prevented development of non-equilibrium water ¯ow, as wasapparent in our experiments.

Summary and conclusions

Figure 8 Soil hydraulic properties of the inter-aggregate pore region
obtained by analysing the in®ltration/redistribution data using thedual-permeability model.

Many studies have shown that hydraulic parameters of thesoil at equilibrium can be obtained by ®tting the traditional
Richards equation to data on transient water ¯ow. Five (r, s,, n, and Ks) or even six (plus l) soil hydraulic parameters
were successfully estimated from, for example, multistepout¯ow (van Dam et al., 1994), evaporation (SÏimuÊnek et al.,1998b), and upward in®ltration (Hudson et al., 1996) experiments. Nevertheless, questions about uniqueness of these
optimized parameters still persist. Applying models of dualporosity or dual-permeability non-equilibrium water ¯owcomplicates the inverse problem by invoking additionalparameters that make the problem less well-posed, if not illposed. Although in our analysis we limited ourselves to eightor nine ®tted parameters, the optimization runs always wererestarted with different initial guesses of the parameters toensure as far as possible reaching the global minimum. Table 3
presents parameters for runs having the smallest values of theobjective function. Still, we have no complete guarantee thatthe global minimum of the objective function was obtained ineach case.
Hudson et al. (1996) used upward in®ltration experiments
for estimating soil hydraulic parameters, apparently without
any problems with non-equilibrium water ¯ow. They carriedout their experiments on repacked samples of a loamy ®nesand, a soil with a fairly large value of n (> 5) typical of soil
having a narrow pore size distributions. Arti®cially packed
sieved soils generally do not have any identi®able structure.
The applied ¯ux at the bottom of the experiment by Hudson
et al. (1996) was about two orders of magnitude less than the
optimized saturated hydraulic conductivity (i.e. the extrapolated value was far beyond the used measurement range). The
much slower in®ltration used by Hudson et al. (1996), and thefact that their soil samples were arti®cially packed, probably

Soil's hydraulic properties are often measured in the laboratory
using various methods of steady-state or transient drying, suchas the pressure plate method, single- or multistep out¯ow, andevaporation methods. Before such experiments are initiated,samples are typically saturated with the common goal ofachieving full saturation (Klute, 1986). Parameters obtained
from these measurements are then often used in simulation
models to predict water ¯ow near the soil surface.
Unfortunately, very little attention is being paid to the factthat the degree of saturation obtained in the laboratory is
seldom achieved in the ®eld, and that preferential ¯ow sotypical for wetting in the ®eld (Flury et al., 1994) can rarely beobserved in experiments based on drying. We have suggested
above a simple laboratory method that may be able to capture
non-equilibrium ¯ow caused by the internal redistribution ofwater in an aggregated macroporous soil from the interaggregate pores to the intra-aggregate pores. The method isbased on modi®cation of the upward in®ltration experiment(Hudson et al., 1996). In®ltration is initiated by a constant orvariable tension (close to zero) at the bottom of the soil
sample. Water then ¯ows fairly quickly upward through thepores between the aggregates from where it subsequently
moves into the aggregates themselves. Fairly fast upward
movement of water was observed using tensiometers at several
depths.
The dual-porosity and dual-permeability models weresuccessfully used to describe the observed non-equilibrium
behaviour. However, our data were too few to optimize aunique set of soil hydraulic parameters needed for calibrationof the fairly complex non-equilibrium ¯ow models. Additional
measurements or modi®cations of the upward tension in®ltration experiment or both may be needed to obtain supplementary information, especially about the mass transfer process.
Initiating the in®ltration afresh after an equilibration period
could be one possible improvement of the method. Testing themethod for soils having different degrees of structuraldevelopment (from a massive repacked soil to a highly
aggregated or macroporous soil) is also desirable.

References
Beven, K. 1991. Modeling preferential ¯ow: an uncertain future? In:
Preferential Flow (eds T.J. Gish & A. Shirmohannadi), pp. 1±11.
American Society of Agricultural Engineers, St Joseph, MI.
Clothier, B.E., Heng, L., Magesan, G.N. & Vogeler, I. 1995. The
measured mobile-water content of an unsaturated soil as a functionof hydraulic regime. Australian Journal of Soil Research, 33,
397±414.
Edwards, W.M., van der Ploeg, R.R. & Ehlers, W. 1979. A numerical# 2001 Blackwell Science Ltd, European Journal of Soil Science, 52, 13±24

24 J. SÏimuÊnek et al.
study of the effects of noncapillary-sized pores upon in®ltration.
Soil Science Society of America Journal, 43, 851±856.
Flury, M., FluÈhler, H., Jury, W.A. & Leuenberger, J. 1994.
Susceptibility of soils to preferential ¯ow of water: a ®eld study.
Water Resources Research, 30, 1945±1954.
Gerke, H.H. & van Genuchten, M.T. 1993. A dual-porosity model for
simulating the preferential movement of water and solutes instructured porous media. Water Resources Research, 29, 305±319.
Gerke, H.H. & van Genuchten, M.T. 1996. Macroscopic representation of structural geometry for simulating water and solute
movement in dual-porosity media. Advances in Water Resources,
19, 343±357.
Gribb, M.M., SÏimuÊnek, J. & Leonard, M.F. 1998. Development of acone penetrometer method to determine soil hydraulic properties.
American Society of Civil Engineers Journal of Geotechnical andGeoenvironmental Engineering, 124, 820±829.
Hudson, D.B., Wierenga, P.J. & Hills, R.G. 1996. Unsaturated
hydraulic properties from upward ¯ow into soil cores. Soil Science
Society of America Journal, 60, 388±396.
Jarvis, N. 1999. Modeling the impact of preferential ¯ow on nonpoint
source pollution. In: Physical Nonequilibrium in Soils: Modelingand Application (eds H.M. Selim & L. Ma), pp. 195±221. Ann
Arbor Press, Chelsea, MI.
Jaynes, D.B., Logsdon, S.D. & Horton, R. 1995. Field method formeasuring mobile/immobile water content and solute transfer ratecoef®cient. Soil Science Society of America Journal, 59, 352±356.
Klute, A. 1986. Water retention: laboratory methods. In: Methods ofSoil Analysis: Part 1, 2nd edn (ed. A. Klute), pp. 635±662.
Agronomy Monograph 9, American Society of Agronomy and Soil
Science Society of America, Madison, WI.
Kool, J.B., Parker, J.C. & van Genuchten, M.T. 1985. Determining
soil hydraulic properties from one-step out¯ow experiments byparameter estimation: I. Theory and numerical studies. Soil Science
Society of America Journal, 49, 1348±1354.
Kool, J.B., Parker, J.C. & van Genuchten, M.T. 1987. Parameter
estimation for unsaturated ¯ow and transport models ± a review.
Journal of Hydrology, 91, 255±293.
Marquardt, D.W. 1963. An algorithm for least-squares estimation of
nonlinear parameters. SIAM Journal of Applied Mathematics, 11,431±441.
Mualem, Y. 1976. A new model for predicting the hydraulicconductivity of unsaturated porous media. Water Resources
Research, 12, 513±522.
Nkedi-Kizza, P., Biggar, J.W., Selim, H.M., van Genuchten, M.T.,
Wierenga, P.J., Davidson, J.M. & Nielsen, D.R. 1984. On theequivalence of two conceptual models for describing ion exchange
during transport through an aggregated oxisol. Water Resources
Research, 20, 1123±1130.
Pruess, K. & Wang, J.S.Y. 1987. Numerical modeling of isothermal
and non-isothermal ¯ow in unsaturated fractured rock ± a review.
In: Flow and Transport through Unsaturated Fractured Rock (eds
D.D. Evans & T.J. Nicholson), pp. 11±22. Geophysics
Monograph42, American Geophysical Union, Washington, DC.

# 2001 Blackwell Science Ltd, European Journal of Soil Science, 52, 13±24

Russo, D., Bresler, E., Shani, U. & Parker, J.C. 1991. Analyses ofin®ltration events in relation to determining soil hydraulic properties by inverse problem methodology. Water Resources Research,
27, 1361±1373.
Santini, A., Romano, N., Ciollaro, G. & Comegna, V. 1995.
Evaluation of a laboratory inverse method for determining
unsaturated hydraulic properties of a soil under different tillage
practices. Soil Science, 160, 340±351.
SÏimuÊnek, J., SÏejna, M. & van Genuchten, M.T. 1998a. The HYDRUS1D software package for simulating water ¯ow and solute transport
in two-dimensional variably saturated media, Version 2.0.
IGWMC-TPS-70, International Ground Water Modeling Center,
Colorado School of Mines, Golden, CO.
SÏimuÊnek, J., Wendroth, O. & van Genuchten, M.T. 1998b. A
parameter estimation analysis of the evaporation method for
determining soil hydraulic properties. Soil Science Society ofAmerica Journal, 62, 894±905.
SÏimuÊnek, J., KodesÏovaÂ, R., Gribb, M.M. & van Genuchten, M.T.
1999a. Estimating hysteresis in the soil water retention function
from cone permeameter experiments. Water Resources Research,
35, 1329±1345.
SÏimuÊnek, J., van Genuchten, M.T., SÏejna, M., Toride, N. & Leij, F.J.
1999b. The STANMOD computer software for evaluating solutetransport in porous media using analytical solutions of convectiondispersion equation, Versions 1.0 and 2.0. IGWMC-TPS-71,
International Ground Water Modeling Center, Colorado School ofMines, Golden, CO.
SÏimuÊnek, J., Wendroth, O. & van Genuchten, M.T. 1999c. Estimating
unsaturated soil hydraulic properties from laboratory tension discin®ltrometer experiments. Water Resources Research, 35,2965±2979.
Toride, N., Leij, F.J. & van Genuchten, M.T. 1995. The CXTFIT codefor estimating transport parameters from laboratory or ®eld tracer
experiments, Version 2.0. Research Report 137, US Salinity
Laboratory, US Department of Agriculture, Agricultural ResearchService, Riverside, CA.
Van Dam, J.C., Stricker, J.N.M. & Droogers, P. 1994. Inverse method
to determine soil hydraulic functions from multistep out¯owexperiment. Soil Science Society of America Proceedings, 58,
647±652.
Van Genuchten, M.T. 1980. A closedform equation for predicting thehydraulic conductivity of unsaturated soils. Soil Science Society ofAmerica Journal, 44, 892±898.
Van Genuchten, M.T. & Wierenga, P.J. 1976. Mass transfer studies insorbing porous media, I. Analytical solutions. Soil Science Societyof America Journal, 40, 473±481.
Van Genuchten, M.T. & Dalton, F.N. 1986. Models for simulating salt
movement in aggregated ®eld soils. Geoderma, 38, 165±183.
Wendroth, O., Ehlers, W., Hopmans, J.W., Kage, H., Halbertsma, J. &WoÈsten, J.H.M. 1993. Reevaluation of the evaporation method for
determining hydraulic functions in unsaturated soils. Soil Science
Society of America Journal, 57, 1436±1443.