M

Modeling Nonequilibrium Flow andTransport Processes Using HYDRUS

S

S

:V

Z

Jirka Šimůnek* and Mar nus Th. van Genuchten

Accurate process-based modeling of nonequilibrium water flow and solute transport remains a major challenge invadose zone hydrology. Our objec ve here was to describe a wide range of nonequilibrium flow and transport modeling approaches available within the latest version of the HYDRUS-1D so ware package. The formula ons range fromclassical models simula ng uniform flow and transport, to rela vely tradi onal mobile-immobile water physical andtwo-site chemical nonequilibrium models, to more complex dual-permeability models that consider both physical andchemical nonequilibrium. The models are divided into three groups: (i) physical nonequilibrium transport models, (ii)chemical nonequilibrium transport models, and (iii) physical and chemical nonequilibrium transport models. Physicalnonequilibrium models include the Mobile-Immobile Water Model, Dual-Porosity Model, Dual-Permeability Model,
and Dual-Permeability Model with Immobile Water. Chemical nonequilibrium models include the One Kine c Site
Model, the Two-Site Model, and the Two Kine c Sites Model. Finally, physical and chemical nonequilibrium transportmodels include the Dual-Porosity Model with One Kine c Site and the Dual-Permeability Model with Two-Site Sorp on.
Example calcula ons using the diﬀerent types of nonequilibrium models are presented. Implica ons for the formulaon of the inverse problem are also discussed. The many diﬀerent models that have been developed over the years fornonequilibrium flow and transport reflect the mul tude of o en simultaneous processes that can govern nonequilibrium and preferen al flow at the field scale.

T

increasing evidence that ﬂow and transport processes in soils often cannot be described using classicalmodels that assume uniform ﬂow and transport (e.g., NkediKizza et al., 1984; Hendrickx and Flury, 2001; Pot et al., 2005;
Köhne et al., 2006). Many laboratory and ﬁeld experiments havedemonstrated the presence of nonequilibrium ﬂow and transportconditions in soils. Nonequilibrium water ﬂow and solute transport in the unsaturated zone can be simulated at present by meansof a large number of models of various degrees of complexity anddimensionality. Modeling approaches range from relatively simple
analytical solutions for solute transport (e.g., van Genuchten,
1981; Toride et al., 1993) to complex numerical codes (e.g.,Šimůnek et al., 2005; Jacques and Šimůnek, 2005). While suchprograms as STANMOD (Šimůnek et al., 1999) that implement
analytical solutions undoubtedly will remain useful for simpliﬁed analyses of solute transport during steady-state ﬂow (e.g., forJ. Šimůnek, Dep. of Environmental Sciences, Univ. of California, Riverside,
CA 92521; M.Th. van Genuchten, USDA-ARS, U.S. Salinity Lab., 450 West
Big Springs Rd., Riverside, CA 92507. Received 17 Apr. 2007. *Corresponding author (Jiri.Simunek@ucr.edu).
Vadose Zone J. 7:782–797
doi:10.2136/vzj2007.0074© Soil Science Society of America
677 S. Segoe Rd. Madison, WI 53711 USA.
All rights reserved. No part of this periodical may be reproduced ortransmi ed in any form or by any means, electronic or mechanical,
including photocopying, recording, or any informa on storage andretrieval system, without permission in wri ng from the publisher.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

782

analyzing solute breakthrough curves measured in the laboratory, or for initial or approximate analysis of ﬁeld-scale transport
problems), numerical models are now increasingly being usedsince they can be applied more readily than analytical models
to realistic laboratory and ﬁeld problems. The use of numericalmodels has been further popularized during the last 20 yr or so
because of the appearance of appropriate software packages inboth the public and commercial domains and the developmentof increasingly sophisticated graphics-based interfaces that can
simplify their use tremendously.
Attempts to describe nonequilibrium transport have traditionally been developed along two lines: physical and chemical
nonequilibrium models (van Genuchten and Cleary, 1979).
While physical nonequilibrium models assume that nonequilibrium ﬂow or transport is caused by physical factors (e.g., vanGenuchten and Wierenga, 1976), chemical nonequilibrium
models assume that chemical factors are the cause of nonequilibrium transport (e.g., Selim et al., 1976; van Genuchten andWagenet, 1989). Only a few researchers (e.g., Brusseau et al.,
1989; Selim et al., 1999; Pot et al., 2005) have combined thephysical and chemical nonequilibrium approaches to accountfor both possible causes of nonequilibrium, thus improving thedescription of solute transport in soils.
Over the years, several publicly available numerical codes
have been developed that consider a number of options for
simulating nonequilibrium water ﬂow and solute transport (e.g.,
Pruess, 1991; Jarvis, 1994; van Dam et al., 1997). Unique tothe HYDRUS-1D software package (Šimůnek et al., 2005,2008) is the wide range of approaches that can be selected forF . 1. Conceptual physical nonequilibrium models for water flow and solute transport. In the plots, θ is the water content, θmo and θim in(b) and (c) are water contents of the mobile and immobile flow regions, respec vely; θM and θF in (d) are water contents of the matrix andmacropore (fracture) regions, respec vely, and θM,mo, θM,im, and θF in (e) are water contents of the mobile and immobile flow regions of thematrix domain, and in the macropore (fracture) domain, respec vely; c are concentra ons of corresponding regions, with subscripts having
the same meaning as for water contents, while S is the total solute content of the liquid phase.

simulating nonequilibrium processes. Our objective here is todescribe the large number of physical and chemical nonequilibrium approaches available in the latest version of HYDRUS-1D(Šimůnek et al., 2008). The models range from classical models
simulating uniform ﬂow and transport, to traditional dualporosity physical and two-site chemical nonequilibrium models,
to complex dual-permeability models that consider both physical and chemical nonequilibrium. Since our focus is mainly
on solute transport and since the nonequilibrium models forwater ﬂow have been reviewed relatively recently (Šimůnek
et al., 2003), we only very brieﬂy review here the governing
water ﬂow equations to deﬁne variables later used in the solutetransport equations. We will present several examples calculated with the diﬀerent nonequilibrium approaches to show
the eﬀect of various transport and reaction parameters, and todemonstrate the consequences of increased complexity in themodels. Implications for the formulation of the inverse problem
are also discussed.

of various water contents and concentrations used in the diﬀerent
models are given in the Appendix.
A hierarchical set of physical nonequilibrium flow andtransport models can be derived from the Uniform Flow Model.
The equilibrium ﬂow and transport model can be modiﬁed byassuming that the soil particles or aggregates have their own
microporosity and that water present in these micropores isimmobile (the Mobile–Immobile Water Model in Fig. 1b and2b). While the water content in the micropore domain is constant in time, dissolved solutes can move into and out of thisimmobile domain by molecular diﬀusion (e.g., van Genuchten
and Wierenga, 1976). This simple modiﬁcation leads to physicalnonequilibrium solute transport while still maintaining uniform
water ﬂow.
The mobile–immobile water model can be further expandedby assuming that both water and solute can move into and outof the immobile domain (Šimůnek et al., 2003), leading to theDual-Porosity Model in Fig. 1c and 2c. While the water content
inside of the soil particles or aggregates is assumed to be constant in the Mobile–Immobile Water Model, it can vary in theDual-Porosity Model since the immobile domain is now allowed
to dry out or rewet during drying and wetting processes. Water
ﬂow into and out of the immobile zone is usually described using
a ﬁrst-order rate process. Solute can move into the immobile
domain of the Dual-Porosity Model by both molecular diﬀusion
and advection with ﬂowing (exchanging) water. Since water can

Overview of Conceptual Models
A large number of alternative physical and chemical nonequilibrium models can be formulated. Figures 1 and 2 show
schematics of a range of possible physical equilibrium and nonequilibrium models for water ﬂow and solute transport. Figure
3 shows similar schematics of various chemical nonequilibrium models that have been incorporated into HYDRUS-1D.
Traditional ﬂow and transport models are
based on the classical description of uniform
ﬂow and transport in soils (the Uniform Flow
Model in Fig. 1a and 2a). In this model, theporous medium is viewed as a collection ofimpermeable soil particles (or of impermeablesoil aggregates or rock fragments), separated
by pores or fractures through which ﬂow
and transport takes place. Variably saturated
water ﬂow through such a porous system isusually described using the Richards equation and solute transport using the classicaladvection–dispersion equation. Deﬁnitions

F . 2. Conceptual physical nonequilibrium models for water flow and solute transport.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008783

F . 3. Conceptual chemical nonequilibrium models for reac ve solute transport. In the plots, θ is the water content, θmo and θim in (d) arewater contents of the mobile and immobile flow regions, respec vely; θM and θF in (e) are water contents of the matrix and macropore (fracture) regions, respec vely; c are concentra ons of the corresponding regions, se are sorbed concentra ons in equilibrium with the liquid
concentra ons of the corresponding regions, and sk are kine cally sorbed solute concentra ons of the corresponding regions.

move from the main pore system into the soil aggregates and viceversa, but not directly between the aggregates themselves, waterin the aggregates can be considered immobile from a larger scalepoint of view.
The limitation of water not being allowed to move directly
between aggregates is overcome in a Dual-Permeability Model(e.g., Gerke and van Genuchten, 1993a,b). Water and solutes
in such models also move directly between soil aggregates, as
shown in Fig. 1d and 2d. Dual-permeability models assume thatthe porous medium consists of two overlapping pore domains,with water ﬂowing relatively fast in one domain (often called themacropore, fracture, or interporosity domain) when close to fullsaturation, and slow in the other domain (often referred to as
the micropore, matrix, or intraporosity domain). Like the DualPorosity Model, the Dual-Permeability Model allows the transfer
of both water and solutes between the two pore regions.
Finally, the Dual-Permeability Model can be further reﬁned
by assuming that inside of the matrix domain an additional
immobile region exists into which solute can move by molecular
diﬀusion (the Dual-Permeability Model with MIM in Fig. 1eand 2e).
Chemical nonequilibrium models implemented into
HYDRUS-1D are schematically shown in Fig. 3. The simplest
chemical nonequilibrium model assumes that sorption is akinetic process (the One Kinetic Site Model in Fig. 3a), usuallydescribed by means of a ﬁrst-order rate equation. This model canbe expanded into a Two-Site Sorption model by assuming thatthe sorption sites can be divided into two fractions (e.g., Selim etal., 1976; van Genuchten and Wagenet, 1989). The simplest twosite sorption model arises when sorption on one fraction of thesorption sites is assumed to be instantaneous, while kinetic sorption occurs on the second fraction (Two-Site Model in Fig. 3b).
This model can be further expanded by assuming that sorption
on both fractions is kinetic and proceeds at diﬀerent rates (theTwo Kinetic Sites Model in Fig. 3c). The Two Kinetic Sites Model
reduces to the Two-Site Model when one rate is so high that it
can be considered instantaneous, to the One Kinetic Site Model
when both rates are the same, or to the chemical equilibrium
model when both rates are so high that they can be considered
instantaneous.
The diﬀerent models discussed thus far involve either physical or chemical nonequilibrium. Many transport situations will
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

involve both types of nonequilibrium. One obvious example
occurs during transport through an aggregated laboratory soilcolumn involving steady-state water ﬂow (either fully saturated
or unsaturated) when both a conservative tracer (no sorption)
and a reactive solute are used. The collected tracer breakthroughcurve may then display typical features reﬂecting nonequilibrium, such as a relatively rapid initial breakthrough followed byextensive tailing. Since the tracer is nonreactive, this nonequilibrium must be caused by physical factors. When the reactive
solute is additionally sorbed kinetically to the solid phase (anindication of a chemical nonequilibrium), the use of a model isrequired that simultaneously considers both physical and chemical nonequilibrium.
The resulting combined physical and chemical nonequilibrium approach may be simulated with HYDRUS-1D using
the Dual-Porosity Model with One Kinetic Site (Fig. 3d). Thismodel considers water ﬂow and solute transport in a dual-porosity system (or a medium with mobile–immobile water) whileassuming that sorption in the immobile zone is instantaneous.
Following the two-site kinetic sorption concept, however, thesorption sites in contact with the mobile zone are now divided
into two fractions, subject to either instantaneous or kinetic
sorption. Since the residence time of solutes in the immobile
domain is relatively large, equilibrium probably exists between
the solution and the sorption complex here, in which case there
is no need to consider kinetic sorption in the immobile domain.
The model, on the other hand, assumes the presence of kinetic
sorption sites in contact with the mobile zone since water canmove relatively fast in the macropore domain and thus prevent
chemical equilibrium.
Finally, chemical nonequilibrium can also be combined with
the Dual-Permeability Model. This last nonequilibrium option
implemented into HYDRUS-1D (the Dual-Permeability Modelwith Two-Site Sorption in Fig. 3e) assumes that equilibrium andkinetic sites exist in both the macropore (fracture) and micropore(matrix) domains.
A complete list of the diﬀerent models summarized here, anddescribed in more detail below, is given in Table 1, including thespeciﬁc equations used for the water ﬂow and solute transportmodels. The table additionally lists the various parameters, andtheir total number, that will appear in any particular solute transport model when used for steady-state ﬂow conditions.
784

1. Equilibrium and nonequilibrium models in HYDRUS-1D, their governing equa ons, and the number of solute transport parameters
for each model under steady-state water flow condi ons.

T

Model name

Water flowequa on(s)Uniform flow and transport
Mobile–immobile water model
Dual-porosity model
Dual-permeability model
Dual-permeability model with immobile water
One kine c site modelTwo-site modelTwo kine c sites model
Dual-porosity model with one kine c site
Dual-permeability model with two-site
sorp onSolute transportequa on(s)At steady-state water flow
No. ofparameters

Equilibrium model[6]
4
Physical nonequilibrium models[1][2], [11]7[2], [3][11]7[4], [5][12]9[4], [5][13], [14]12Chemical nonequilibrium models[1][6], [15]5[1][16], [17]
6
[1]
[18]
7
Physical and chemical nonequilibrium models[2], [3][2], [20], [21]9[4], [5]
[22]
13
[1]

Parameters
θ, q, λ, Kdθ mo, θ im, qmo, λ mo, Kd, fmo, ω mim
θ mo, θ im, qmo, λ mo, Kd, fmo, ω mim
θ m, θ f, qm, qf, λ m, λf, Kdm, Kdf, ω dp
θ m,m, θ im,m, θ f, qm, qf, λ m, λf, Kdm, Kdf, fm, ω dp, ω pdm
θ, q, λ, Kd, α k
θ, q, λ, Kd, fe, α k
θ, q, λ, ka1, kd1, ka2, kd2
θ mo, θ im, qmo, λ mo, Kd, fmo, fem, ω ph, α ch
θ m, θ f, qm, qf, λ m, λf, Kdm, Kdf, fm, f f, ω dp, α ch,m, α ch,fSpecific Models for Water Flowtion leads to two-region type ﬂow and transport models (vanGenuchten and Wierenga, 1976) that partition the liquid phase
into mobile (ﬂowing, interaggregate), θmo, and immobile (stagnant, intraaggregate), θim, regions [L3 L−3]:

We will brieﬂy review the governing equations for both equilibrium (uniform) water ﬂow and for nonequilibrium ﬂow in
dual-porosity and dual-permeability systems. Subsequently wewill describe various solute transport models that consider either
physical or chemical nonequilibrium and models that consider
simultaneously both physical and chemical nonequilibrium.

θ = θ mo + θ im

The dual-porosity formulation for water ﬂow can be based on amixed formulation of the Richards Eq. [1] to describe water ﬂow
in the macropores (the preferential ﬂow pathways) and a massbalance equation to describe moisture dynamics in the matrix asfollows (Šimůnek et al., 2003):

Uniform Flow Model

Numerical models for water ﬂow in soils (Fig. 1a and 2a) areusually based on the following equation:⎛ ∂h⎞⎤∂θ( h )∂ ⎡⎢ K ( h )⎜⎜ + 1⎟⎟⎟⎥ − S
=
[1]
⎝ ∂z
⎠⎦⎥
∂t
∂ z ⎣⎢or its extensions (e.g., for two- and three-dimensional systems). In
Eq. [1], often referred to as the Richards equation, z is the verticalcoordinate positive upward [L], t is time [T], h is the pressure
head [L], θ is the water content [L3 L−3], S is a sink term representing root water uptake or some other source or sink [T−1],
and K(h) is the unsaturated hydraulic conductivity function,
often given as the product of the relative hydraulic conductivity,
Kr (dimensionless), and the saturated hydraulic conductivity, Ks[L T−1]. Solutions of the Richards Eq. [1] require knowledge ofthe unsaturated soil hydraulic functions made up of the soil water
retention curve, θ(h), which describes the relationship between
the water content θ and the pressure head h, and the unsaturated
hydraulic conductivity function, K(h), which deﬁnes the hydraulic conductivity K as a function of h or θ. HYDRUS-1D considers
both relatively traditional models (Brooks and Corey, 1964; van
Genuchten, 1980) for the hydraulic functions, as well as more
recent alternative single- (e.g., Kosugi, 1996) and dual-porosity(Durner, 1994) models.

∂θ mo (hmo )
∂t
⎛ ∂h
⎞⎤
∂ ⎡
⎢ K (hmo )⎜ mo + 1⎟⎟⎟⎥ − S mo (hmo )−Γ w
=
⎜
⎝ ∂z
⎠⎦⎥
∂ z ⎣⎢[3]∂θ im (him )= −S im (him )+ Γ w∂twhere Smo and Sim are sink terms for the mobile and immobile regions, respectively [T−1], and Γw is the transfer rate forwater between the inter- and intraaggregate pore domains [T−1].
Šimůnek et al. (2003) and Köhne et al. (2004) discussed diﬀerent
formulations that can be used to evaluate the mass transfer rate Γw.
Dual-Permeability ModelDiﬀerent dual-permeability approaches (Fig. 1d, 1e, 2d, and2e) may be used to describe ﬂow and transport in structured
media. While some models invoke similar equations for ﬂow in
the fracture and matrix regions, others use diﬀerent formulations for the two regions. A typical example of the ﬁrst approach,
implemented in HYDRUS-1D, is the work of Gerke and van
Genuchten (1993a,b, 1996), who applied the Richards equation to each of the two pore regions. The ﬂow equations for themacropore (fracture) (subscript f ) and matrix (subscript m) poresystems in their approach are given by

Dual-Porosity Model

Dual-porosity models (Fig. 1c and 2c) assume that waterﬂow is restricted to the macropores (or interaggregate pores andfractures), and that water in the matrix (intraaggregate poresor the rock matrix) does not move at all. This conceptualizawww.vadosezonejournal.org · Vol. 7, No. 2, May 2008

[2]

785

∂θ f (hf ) ∂ ⎡
⎛ ∂h
⎞⎤
Γ
⎢ K f (hf )⎜ f + 1⎟⎟⎟⎥ − S f (hf )− w
=
⎜
⎝ ∂z
⎠⎥⎦
∂t
∂ z ⎢⎣w∂θ m (hm )
∂t
⎛ ∂h
⎞⎤
∂ ⎡
Γ
⎢ K m (hm )⎜⎜ m + 1⎟⎟⎟⎥ − S m (hm )+ w
=
⎝ ∂z
⎠⎥⎦
∂ z ⎢⎣
1− w

In our examples below we will be using Eq. [7] for linear sorption and applying the resulting deﬁnition of R to all of theliquid domains involved (e.g., the total liquid phase, the mobile
and immobile regions, or the matrix and macropore domains)with their appropriate parameters. Although considered in
HYDRUS-1D, the eﬀect of molecular diﬀusion will be neglectedin the various examples. The dispersion coeﬃcient D accounting
only for hydrodynamic dispersion [L2 T−1] is thus deﬁned as

[4]

q
[9]
θ
where λ is the dispersivity [L] and v the average pore velocity [LT−1]. The same deﬁnition of the dispersion coeﬃcient will be
used for all mobile phases. HYDRUS-1D additionally considers
molecular diﬀusion for transport in the gaseous phase, which will
not be discussed here either. Finally, in the examples to follow wewill use the time T needed to reach one pore volume of eﬄuent.
This time is given by
D = λv = λwhere w is the ratio of the volumes of the macropore or fracturedomain and the total soil system (dimensionless). This approachis relatively complicated in that the model requires characterization of water retention and hydraulic conductivity functions(potentially of diﬀerent form) for both pore regions, as well as a
hydraulic conductivity function of the fracture–matrix interface.
Note that the water contents, θf and θm in Eq. [4], have diﬀerent
meanings than in Eq. [3], where they represented water contentsof the total pore space (i.e., θ = θmo + θim), while here they referto water contents of the two separate (fracture or matrix) poredomains such thatθ = θ F + θ M = w θ f + (1− w ) θ m

θL
[10]
q
where L is some distance (e.g., column length) within the transport domain being considered. For the physical nonequilibrium
models discussed below, the pore volume for diﬀerent domains is
deﬁned using corresponding water contents and ﬂuxes.
T=[5]

Hence, lowercase subscripts in the dual-permeability model referto the local (pore-region) scale, while uppercase subscripts referto the global (total soil medium) scale.

Physical Nonequilibrium Transport Models

Specific Models for Solute Transport

Mobile–Immobile Water and Dual-Porosity Models

The concept of two-region, dual-porosity type solute transport (Fig. 1b, 1c, 2b, and 2c) was implemented already in earlier
versions (1.0 and 2.0) of HYDRUS-1D to permit consideration of physical nonequilibrium transport. While the physicalnonequilibrium transport model in the earlier versions was
combined only with uniform water ﬂow Eq. [1], Version 3.0 of
HYDRUS-1D was expanded to also consider the dual-porosity
water ﬂow model Eq. [3] with a transient immobile water content.
In both implementations, the governing solute transport equations are as follows:

Uniform Transport

Solute transport in numerical models is usually describedusing the relatively standard advection–dispersion equation (Fig.
1a and 2a) of the form∂θ c∂ s ∂ ⎜⎛∂ c ⎞⎟ ∂ qc
[6]
+ρ =
−φ
⎟−
⎜ θD
∂t
∂ t ∂ z ⎝ ∂ z ⎠⎟ ∂ z
or various extensions thereof (e.g., for two- and three-dimensional systems, or for multiple phases or components). In Eq.
[6], c is the solution concentration [M L−3], s is the sorbed concentration [M M−1], D is the dispersion coeﬃcient accounting
for both molecular diﬀusion and hydrodynamic dispersion [L2T−1], q is the volumetric ﬂuid ﬂux density [L T−1] evaluatedusing the Darcy–Buckingham law, and φ is a sink–source termthat accounts for various zero- and ﬁrst-order or other reactions[M L−3 T−1].
While HYDRUS-1D considers a general nonlinear sorption
equation that can be simpliﬁed into a Langmuir or Freundlich
isotherm (Šimůnek et al., 2005), for simplicity we assume here
only linear adsorption of the forms = Kdc∂θ moc mo∂s+ f moρ mo =
∂t
∂t
∂ ⎛⎜
∂ c mo ⎞⎟ ∂ q moc mo−φ mo −Γ s⎟−
⎜ θ mo Dmo∂z ⎝∂ z ⎠⎟∂zρK dθ

www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

∂θ imc im∂s+ (1− f mo )ρ im = Γ s −φ im
∂t
∂t

[11b]

Γ s = ω mim (c mo − c im )+ Γ w c *[11c]

in which solute exchange between the two liquid regions is modeled as the sum of an apparent ﬁrst-order diﬀusion process andadvective transport (where applicable). In Eq. [11], cmo and cim
are concentrations of the mobile and immobile regions [M L−3],
respectively; smo and sim are sorbed concentrations of the mobile
and immobile regions [M M−1], respectively; Dmo is the dispersion coeﬃcient in the mobile region [L2 T−1], qmo is thevolumetric ﬂuid ﬂux density in the mobile region [L T−1], φmoand φim are sink–source terms that account for various zero- and[7]

where Kd is the distribution coeﬃcient [L3 M−1]. Linear sorption leads to the following deﬁnition of the retardation factor R(dimensionless):
R = 1+

[11a]

[8]

786

ﬁrst-order or other reactions in both regions [M L−3 T−1]; fmo
is the fraction of sorption sites in contact with the mobile watercontent (dimensionless), ωmim is the mass transfer coeﬃcient[T−1], and Γs is the mass transfer term for solutes between themobile and immobile regions [M L−3 T−1]. Eauation [11a]describes solute transport in the mobile (macropore) zone, Eq.
[11b] is a mass balance for the immobile (micropore) domain,
while Eq. [11c] (Γs) describes the rate of mass transfer between
the mobile and immobile domains. The second (advective) termof Γs in Eq. [11] is equal to zero for the Mobile–Immobile Model
since that model does not consider water ﬂow between the tworegions. In the Dual-Porosity Model, c* is equal to cmo for Γw >0 and cim for Γw < 0.

∂θ m,mc m,m∂s+ ρ m f m m,m =
∂t
∂t
[14b]
∂ c m,m ⎞⎟ ∂ q mc m,m∂ ⎛⎜Γ
−φ m,m + s −Γ s *⎟⎟−⎜⎜ θ m,m Dm
1− w
∂z ⎝
∂z ⎠
∂z

Analogous to Eq. [4], the dual-permeability formulation forsolute transport is based on advection–dispersion type equations
for transport in both the fracture and matrix regions as follows(Gerke and van Genuchten, 1993a,b) (Fig. 1d and 2d):[12a]

∂θ mc m∂s
+ρ m =∂t∂t∂ ⎛⎜∂ c m ⎞⎟ ∂ q mc mΓ
θ m Dm−φ m + s
⎟⎟−
⎜
∂z ⎝
∂z ⎠
∂z
1− w

[12b]

Γ s = ω dp (1− w ) θ m (c f − c m )+ Γ w c *

[12c]

[14d]

Γ s * = ω dpm (c m,m − c im,m )[14e]

One Kine c Site Model

When sorption in the Uniform Transport Model is considered (Fig. 3a) to be kinetic, Eq. [6] needs to be supplemented
with an equation describing the kinetics of the sorption process.
This is usually done by assuming a ﬁrst-order process as follows:
∂θ c
∂s k
∂ ⎛⎜
∂ c ⎞⎟ ∂ qc
+ρ
=
θD
−φ
⎟−
⎜
∂t
∂t
∂ z ⎝ ∂ z ⎠⎟ ∂ z∂s k= α k ρ (s ek − s k )−φ k∂ts ek = K d cρDual-Permeability Model with Immobile Water

The Dual-Permeability Model with Immobile Water (Fig. 1eand 2e) assumes that the liquid phase of the matrix can be further
partitioned into mobile (ﬂowing), θm,m [L3 L−3], and immobile(stagnant), θim,m [L3 L−3], regions as follows:[15]

where s ek is the sorbed concentration that would be reached at
equilibrium with the liquid-phase concentration [M M−1], sk isthe sorbed concentration of the kinetic sorption sites [M M−1],
αk is a ﬁrst-order rate constant describing the kinetics of thesorption process [T−1], and φk represents a sink–source termthat accounts for various zero- and ﬁrst-order or other reactions
at the kinetic sorption sites [M L−3 T−1].

[13]

where θm is the volumetric water content of the matrix poresystem [L3 L−3]. The governing advection–dispersion equation
for transport in the matrix region (Eq. [12b]) is then replacedwith the modiﬁed equations (Eq. [11]) (e.g., Pot et al., 2005)
to yield∂θ f c f∂s∂ ⎜⎛∂c ⎞ ∂q cΓ+ ρf f =θ f Df f ⎟⎟⎟− f f −φ f − s
⎜
[14a]
∂t
∂t
∂z ⎝
∂z ⎠
∂z
w

www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

Γ s = ω dp (1− w ) θ m (c f − c m,m )+ Γ w c *Chemical Nonequilibrium Transport Models

The variables in Eq. [12] have similar meanings as in Eq. [11],except that they refer now to two overlapping domains, i.e., thematrix (subscript m) and fracture (subscript f ) domains. Equation[12a] describes solute transport in the fracture domain, Eq. [12b]transport in the matrix domain, and Eq. [12c] advective–dispersive mass transfer between the fracture and matrix domains.
Equation [12] assumes complete advective–dispersive transport
descriptions for both the fractures and the matrix. Van Genuchten
and Dalton (1986) and Gerke and van Genuchten (1996), among
others, discussed possible expressions for the ﬁrst-order solute
mass transfer coeﬃcient, ωdp [T−1].

θ m = θ m,m + θ im,m[14c]where cim,m and cm,m are solute concentrations in the immobile
and mobile zones of the matrix region [M L−3], respectively;
φm,m and φim,m represent various reactions in the mobile andimmobile parts of the matrix [M L−3 T−1], respectively; fm isagain the fraction of sorption sites in contact with the mobile
region of the matrix (dimensionless), ωdpm is the mass transfer
coeﬃcient between mobile and immobile zones of the matrix
region [T−1], and Γs* is the mass transfer term for solutes
between the mobile and immobile regions of the matrix domain[M L−3 T−1]. Equation [14a] now describes solute transport in
the fracture domain, Eq. [14b] transport in the mobile zone of
the matrix domain, Eq. [14c] is a mass balance for the immobile
zone of the matrix domain, Eq. [14d] describes mass transfer between the fracture and matrix domains, while Eq. [14e]describes mass transfer between the mobile and immobile zones
within the matrix domain.

Dual-Permeability Model∂θ f c f∂s∂ ⎛⎜∂ c f ⎞⎟ ∂ q f c f
Γ+ρ f =−φ f − s⎟⎟−⎜⎝ θ f Df
⎠
∂t
∂t
∂z
∂z
∂z
w

∂θ im,mc im,m∂s
+ ρ m (1− f m ) im,m = Γ s * −φ im,m∂t∂t

Two-Site Model

Similarly to the mobile–immobile water concept, the concept of two-site sorption (Selim et al., 1976; van Genuchten andWagenet, 1989) (Fig. 3b) was implemented already in Versions
1.0 and 2.0 of HYDRUS-1D to permit consideration of nonequilibrium adsorption–desorption reactions. The two-site sorption

787

concept assumes that the sorption sites can be divided into two
fractions:s = se + sk

blocking to attachment sites or depth-dependent straining thatare considered in HYDRUS-1D (e.g., Bradford et al., 2004).
It is easily shown that the formulation based on attachment–detachment coeﬃcients is mathematically identical to
the formulation using ﬁrst-order mass transfer coeﬃcients. For
example, by comparing Eq. [15] with [18] we have[16]

Sorption se [M M−1], on one fraction of the sites (Type 1 sites) is
assumed to be instantaneous, while sorption sk [M M−1], on theremaining (Type 2) sites is considered to be a ﬁrst-order kinetic
rate process. The system of equations describing the Two-Site
Model is given by:
∂θ c
∂s e
∂s k
∂ ⎛⎜
∂ c ⎞⎟ ∂ qc
+ρ
+ρ
=
−φ
⎟−
⎜ θD
∂t
∂t
∂t
∂ z ⎝ ∂ z ⎠⎟ ∂ z

[17a]

Kd =

s e = f e K dc

[17b]

∂s k
= α k ρ (s ek − s k )−φ k
∂t

[17c]

ρ

s ek =(1− f e )K d c∂s k= θ ka c − kd ρ s k = ρα k (K d c − s k )∂tρka = α k K dθ
kd = α k
ρ

θ ka
ρα k

Physical and Chemical Nonequilibrium Transport Models
Dual-Porosity Model with One Kine c Site

This model (Fig. 3d) is similar to the Dual-Porosity Model(Eq. [2]) in that the porous medium is divided into mobile andimmobile domains such that θ = θmo + θim. The current model,
however, additionally divides the sorption sites in contact withthe mobile zone, similarly to the Two-Site Model (Eq. [16]), into
two fractions involving instantaneous and kinetic sorption suchthat the total sorption concentration at equilibrium is given by[17d]where fe is the fraction of exchange sites assumed to be in
equilibrium with the liquid phase (dimensionless), and αk is aﬁrst-order rate constant [T−1]. Equation [17a] describes solutetransport in the total system, Eq. [17b] equilibrium sorption onto
the instantaneous sorption sites, Eq. [17c] is a mass balance of
the kinetic sorption sites (van Genuchten and Wagenet, 1989),
while Eq. [17d] represents the sorbed concentration of the kinetic
sites when equilibrium would be reached with the liquid-phase
concentration.

s = (1− f mo )s im + f mo s moek
= (1− f mo )s im + f mo (s mo+ s mo,e)= (1− f mo )K d c mo + f mo f em K d c mo

[20]

+ f mo (1− f em )K d c mo
= K d c mo

Two Kine c Sites Model

To facilitate simulations of the transport of colloids or
microorganisms (such as viruses and bacteria), Version 3.0 of
HYDRUS-1D also implemented a Two Kinetic Sites Model (Fig.
3c) using the attachment–detachment approach:ewhere s mo
is the sorbed concentration in equilibrium with theliquid-phase concentration of the mobile region of the Dualk
Porosity Model [M M−1], s mo,eis the sorbed concentration of
the kinetic sites in contact with the mobile region of the DualPorosity Model when at equilibrium [M M−1], fmo is the fractionof sorption sites in contact with mobile water (the remainder is incontact with immobile water), and fem is the fraction of sorption
sites in equilibrium with the mobile liquid phase (the remaining
kinetic sites are also in contact with the mobile liquid phase). The
complete Dual-Porosity Model with One Kinetic Site is describedusing the following equations:

∂θ c
∂s k
∂s k
∂ ⎜⎛
∂ c ⎞⎟ ∂ qc
+ρ 1 +ρ 2 =−φ⎟−⎜⎝ θ D∂t∂t∂t∂z∂ z ⎠⎟ ∂ z
∂ s1k
[18]
= ka1θ c − kd1ρ s1k −φ k1∂t∂s kρ 2 = ka2 θ c − kd2ρ s 2k −φ k2∂twhere s1k and s 2k are sorbed concentrations of the ﬁrst andsecond fractions of kinetic sorption sites [M M−1], respectively;
ka1 and ka2 are attachment coeﬃcients for the ﬁrst and second
fractions of kinetic sorption sites [T−1], respectively; kd1 andkd2 are detachment coeﬃcients for the ﬁrst and second fractionsof kinetic sorption sites [T−1], respectively; and φk1 and φk2
represent sink–source terms for the ﬁrst and second fractionsof kinetic sorption sites [M L−3 T−1], respectively. Note that
the Two Kinetic Sites Model can be used (and often is used) to
describe diﬀerent processes. While the ﬁrst kinetic process could
be chemical attachment, the second kinetic process could represent physical straining (e.g., Bradford et al., 2004; Gargiulo et al.,2007, 2008). Note that in Eq. [18] we do not give the nonlinear
blocking coeﬃcients accounting for, for example, Langmuirianρ

www.vadosezonejournal.org · Vol. 7, No. 2, May 2008[19]

∂θ moc mo∂s e+ f moρ mo =
∂t
∂t
⎛
∂ ⎜
∂c ⎞ ∂q c
θ mo Dmo mo ⎟⎟⎟− mo mo −φ mo −Γ s1 −Γ s2
⎜
∂z ⎝
∂z ⎠
∂z

[21a]

∂θ imc im∂s+ (1− f mo )ρ im = Γ s1 −φ im
∂t
∂t

[21b]

f moρ

k
∂ s mo= Γ s2 −φ mo,k∂t

Γ s1 = ω ph (c mo − c im )

788

[21c]
[21d]

k
k
Γ s2 = α ch ρ (s mo,e
− s mo
)

[21e]

e
s mo
= f em K d c mo

[21f ]ks mo,e= (1− f em )K d c momatrix and fracture domains, respectively; φf, φm, φf,k, and φm,krepresent reactions in the equilibrium phases of the fracture andmatrix domains and at the kinetic sites of the fracture and matrix
domains [M L−3 T−1], respectively; and αch,m and αch,f are again
ﬁrst-order rate constants for the matrix and fracture domains[T−1], respectively. Note that the distribution coeﬃcients can be
diﬀerent in the diﬀerent regions (i.e., Kdf ≠ Kdm).

[21g]

are ﬁrst-order rate constants [T−1] account-where ωph and αch
ing for physical and chemical rate processes, respectively; Γs1 is
the mass transfer term for solute exchange between the mobile
and immobile regions [M L−3 T−1]; Γs2 represents mass transfer
to the kinetic sorption sites in the mobile region [M L−3 T−1];
and φmo, φim, and φmo,k represent sink–source terms for theequilibrium phases in the mobile zone, the immobile zone, andthe kinetic sorption sites [M L−3 T−1], respectively. Equation[21a] describes transport in the mobile phase, Eq. [21b] is a massbalance for the immobile phase, and Eq. [21c] a mass balance
for the kinetic sorption sites in contact with the mobile zone.
Equations [21d] and [21e] describe mass transfer rates between
the mobile and immobile zones and to the kinetic sorption sites,
respectively, while Eq. [21f ] and [21g] represent sorption onto
the equilibrium and kinetic sorption sites in contact with themobile zone, respectively.

Numerical ImplementaƟon
All of the models presented here were implemented in the
HYDRUS-1D software package and as such are all solved withvery similar numerical techniques. The Galerkin-type linear ﬁnite
element method was used for spatial discretization of the governing partial diﬀerential equations, while ﬁnite diﬀerence methods
were used to approximate temporal derivatives. A fully implicit
ﬁnite diﬀerence scheme with Picard linearization was used tosolve the Richards equation, while a Crank–Nicholson ﬁnite
diﬀerence scheme was used for solution of the advection–dispersion equations. For the dual-permeability models, we always
ﬁrst solved the equations describing processes in the matrix, after
which the equations describing processes in the fractures weresolved. Complete details about the invoked numerical techniques
are provided in the HYDRUS-1D technical manual (Šimůnek
et al., 2008).

Dual-Permeability Model with Two-Site Sorp on

Finally, simultaneous physical and chemical nonequilibrium
processes are implemented in HYDRUS-1D by assuming applicability of the Dual-Permeability Model (Gerke and van Genuchten,
1993a; Šimůnek et al., 2003) and dividing the sorption sites ofboth the fracture and matrix domains into equilibrium andkinetic sites (Fig. 3e). This model leads to the following set of
equations (Pot et al., 2005):

ImplicaƟons for FormulaƟon of the
Inverse Problem
The physical and chemical nonequilibrium models both
involve relatively large numbers of parameters, many of whichcannot be (or cannot easily be) measured independently. Someparameters need to be obtained by calibrating a particular model
against laboratory or ﬁeld measurements. This is usually doneusing a parameter estimation procedure (e.g., Šimůnek andHopmans, 2002; Šimůnek et al., 2002) in which the sum ofsquared deviations between measurements and model predictions,organized in an objective function, is minimized. While the deﬁnition of the objective function for models that consider only
uniform ﬂow and transport may be relatively straightforward, theobjective functions for calibration of the nonequilibrium models
can become extremely complicated. For the equilibrium ﬂow andtransport models, one usually deﬁnes the objective function interms of measured pressure heads, water contents, solution concentrations, or actual or cumulative water or solute ﬂuxes. Notefrom Fig. 1a that the system is always described using uniquevalues of the water content, pressure head, or solution concentration. In addition to the resident concentration, the objectivefunction can also be deﬁned using the ﬂux concentration or the
total solute mass at a speciﬁed location. The total solute mass
must include not only the mass in the liquid phase but also themass sorbed to either instantaneous or kinetic sorption sites (thelatter if a kinetic sorption model is used).
By comparison, when physical nonequilibrium models areused, one can encounter diﬀerent water contents such as themobile water content, the immobile water content, the water content of the fracture domain, the water content of the matrix, orthe total water content. The nonequilibrium models typically will
also involve diﬀerent types of concentrations, including concentrations of the mobile and immobile zones, sorbed concentrations∂θ f c f∂s e+ ρf f =∂t∂t∂ ⎛⎜∂ c f ⎞⎟ ∂ q f c f
Γ−φ f − s −Γ f⎟⎟−⎜⎝ θ f Df
⎠
∂z
∂z
∂z
w
∂θ mc m∂s e
+ ρm m =∂t∂t∂ ⎛⎜∂ c m ⎞⎟ ∂ q mc mΓ
θ m Dm−φ m + s − Γ m⎟−∂ z ⎜⎝∂ z ⎠⎟
∂z
1− w

[22]
∂ s fk= Γ f −φ f,k∂t∂s kρ m = Γ m −φ m,k∂t
Γ s = ω dp (1− w ) θ m (c f − c m )+ Γ w c *ρΓ f = ρ f α ch,f ⎡⎢(1− f f )K df c f − s fk ⎤⎥⎣⎦⎡Γ m = ρ mα ch,m ⎢(1− f m )K dmc m − s mk ⎤⎥⎣⎦where s mk and s fk are sorbed concentrations of Type 2 (kinetic)sites in the matrix and fracture domains [M M−1], respectively;
fm and ff are fractions of the exchange sites assumed to be in
equilibrium with the solution phases (dimensionless) of thewww.vadosezonejournal.org · Vol. 7, No. 2, May 2008789

the case most resembling uniform transport in the mobile region,indeed reaches approximately 0.5 after about 4 d. Notice also that
this breakthrough curve shows the most tailing resulting from theslow release of solute from the mobile zone into immobile liquid.
The large value of the mass transfer coeﬃcient (10 d−1) leads tofast equilibration of the concentrations of the mobile and immobile
zones, and thus to a breakthrough curve resembling uniform transport in the entire pore system (T = 1.666 d and R = 4) characterized
by a relatively sigmoidal curve. If the mass transfer coeﬃcient had
been equal to zero (i.e., equilibrium ﬂow in the mobile region), theresulting breakthrough curves would have been the same as in Fig.
4. Figure 5 (right) demonstrates the eﬀect of an increasing fraction of immobile water and the corresponding fraction of sorption
sites. Following Nkedi-Kizza et al. (1983), the ratios of mobile tototal water (θmo/θ) and equilibrium to total sorption sites (f) wereassumed to be the same. A smaller fraction of mobile water leads toearlier solute arrival and more pronounced tailing. Notice that the
red lines in Fig. 5 (left and right) represent the same simulation.

associated with instantaneous and kinetic sorption sites, and total
concentrations. Since diﬀerent measurement techniques may leadto diﬀerent types of concentrations (e.g., resident, ﬂux, or total
concentrations), a ﬂexible inverse approach should also allow for
the diﬀerent modes of concentrations. Diﬀerent types of waterand solute ﬂuxes can furthermore be deﬁned in the nonequilibrium models, including total ﬂuxes, ﬂuxes in the mobile zone,
or ﬂuxes in the fracture and matrix domains. To provide users
with ﬂexibility during the model calibration process, these various water contents, concentrations, and water and solute ﬂuxes
must be considered in the objective function. Table 2 provides alist of variables in the diﬀerent equilibrium and nonequilibrium
models that can be incorporated in the objective function when
HYDRUS-1D is used for parameter estimation purposes.

ApplicaƟons
We now brieﬂy demonstrate the main features of selectedtransport models, including especially the sensitivity of calculated
breakthrough curves to several key transport parameters. As before,
we will focus only on solute transport, while assuming that waterﬂow is at full saturation and steady state and thus that the various
water contents and ﬂuxes are constant with depth and time. Steadystate ﬂow is initiated by ﬁxing all pressure heads in the system equal
to zero (both the initial and boundary conditions).

Dual-Permeability ModelIn the dual-permeability system, solute moves simultaneously through two overlapping porous regions, with the numberof model parameters further increasing. Solute transport during
steady-state water ﬂow now requires the following nine parameters: two water contents (θm and θf ), two ﬂuxes (qm and qf ),two dispersivities (λm and λf ), distribution coeﬃcients for eachdomain (Kdm and Kdf ), and a mass transfer coeﬃcient ωdp. The
exchange of solute between the two pore regions is proportional
to the mass transfer coeﬃcient ωdp. When this coeﬃcient is equal
to zero (black lines in Fig. 6), solute will move independently
through each of the two pore systems, in which case the transport
process reduces to that of the uniform transport model applied
separately to each of the two pore systems. Outﬂow concentrations from the matrix (thin line), fracture (intermediate line),
and the entire soil (thick line) are shown in Fig. 6. Note that
due to accelerated transport in the dual-permeability system, thetime scale in Fig. 6 is only 10 d, compared with 20 d in thepreceding examples. For the selected parameters, the retardation
factor for both regions is equal to 4 and the pore volumes areequal to 0.1666 and 1.666 d for the fracture and matrix domains,respectively. The solute front indeed arrives after about 0.6666
and 6.666 d (TR) in the fracture and matrix domains, respectively. The average outﬂow concentration can now be obtained
as a weighted average of the outﬂow concentrations from eachdomain:Uniform Transport

The uniform (equilibrium) solute transport model during
steady-state water ﬂow is fully deﬁned in terms of four parameters(Table 1). Of these, two parameters (the water content, θ, andthe ﬂuid ﬂux density, q) are related to water ﬂow, and only two
parameters (the dispersivity, λ, and the distribution coeﬃcient,Kd) are directly related to solute transport. Figure 4 shows thewell-known eﬀects of the dispersivity (left) and the distribution
coeﬃcient (right) on calculated breakthrough curves. Notice thatthe red lines in Fig. 4 (left and right) represent the same simulation. While higher values of the dispersivity lead to earlier arrivalof the solute and more pronounced tailing compared with lower
λ values, larger values of the distribution coeﬃcient delay thearrival of solute in the eﬄuent. Note that the selected distribution coeﬃcients of 0, 1, and 3 lead to retardation factors of 1, 4,
and 10, respectively. Also note that one pore volume (Eq. [10])for the selected parameters corresponds to 1 d.
Physical Nonequilibrium Transport Models
Mobile–Immobile Water and Dual-Porosity Models

Since the immobile water content is constant for conditionsof steady-state water ﬂow for both the Mobile–Immobile Waterand the Dual-Porosity models, the two models are identical withrespect to solute transport. These two models have three additional parameters compared with the Uniform Transport Model:
the immobile water content (θim), the fraction of sorption sites incontact with mobile water (fmo), and the mass transfer coeﬃcient(ωmim) (Table 1). Figure 5 demonstrates the eﬀect of the masstransfer coeﬃcient (left) and the fraction of mobile water (right)
on calculated breakthrough curves. For the breakthrough curves
in Fig. 5 (left), one pore volume (T = Lθmo/q) is equal to 1 d andthe retardation factor of the mobile phase (R = 1 + fmoρKd/θmo)is equal to 4. Notice that the eﬄuent concentration for the casewith the lowest value of the mass transfer coeﬃcient (0.1 d−1), i.e.,
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

wc f q f + (1− w ) c m q m[23]
wq f + (1− w ) q m
As shown in Fig. 6, the average outﬂow concentration quickly
increases initially, similarly as for the fracture domain, but thenstabilizes at about 0.5 since only about half of the outﬂow comesfrom the fracture domain. When solute arrives also in the outﬂow
from the matrix domain, the average outﬂow concentration gradually increases again until it reaches unit concentration. The solute
breakthrough curves start deviating from those calculated using the
Uniform Flow Model individually for the matrix and macropore
regions when the mass transfer coeﬃcient is increased (red lines inFig. 6) because of exchange of solute between the two regions. For
relatively large values of the mass transfer coeﬃcient, this causes
cT =

790

T

2. Variables that can be included in the objec ve func on.

Model
Uniform model

Dual-Porosity

Variable
Water content
Pressure head
Resident concentra on
Flux concentra onθhcDefini on

Total solute mass (uniform transport)
Total solute mass (One Kine c Site Model)
Total solute mass (Two-Site Model)
Total solute mass (Two Kine c Sites Model)c(θ + ρKd)cθ + ρskc(θ + fe ρKd) + ρsk

Water content
Pressure head
Mobile-zone resident concentra on
Immobile-zone resident concentra on
Mobile-zone flux concentra onθ = θ mo + θ im
hmo
cmo
cim

Average liquid concentra onc θ +c θ
c Tl = mo mo im imθTotal solute mass (no chemical nonequilibrium)c mo (θ mo + f moρ K d )+ c im ⎣⎡ θ im + (1− f mo )ρ K d ⎦⎤Total solute mass (kine c sorp on)

c mo (θ mo + f mo f em ρ K d )+c−Dθ ∂ cq ∂zc θ + ρ s1k + ρ s 2k

c mo −

Dmo θ mo ∂ c mo
q mo
∂z

k
+ c im ⎣⎡ θ im + (1− f mo )ρ K d ⎦⎤f moρ s mo

Dual-Permeability

Dual-Permeability with
MIM

Matrix water content
Fracture water content
Total water content
Matrix pressure head
Fracture pressure head
Matrix flux
Fracture flux
Total flux
Fracture concentra on
Matrix concentra on
Flux concentra onθ M = (1 − w)θ mθ F = wθ f
θ = wθ f + (1 − w)θ mhm
hf
qm(1 − w)
q fw
qm(1 − w) + qfw
cfcm

Average liquid concentra on

wc f θ f + (1− w )c m θ m
w θ f + (1− w )θ m

Total solute mass (no chemical nonequilibrium)
Total solute mass (kine c sorp on)

wc f(θ f + ρ fKdf ) + (1 − w)cm(θ m + ρ mKdm)

Matrix water content
Fracture water content
Total water content
Matrix pressure head
Fracture pressure head
Matrix flux
Fracture flux
Total flux
Fracture concentra on
Matrix immobile concentra on
Matrix mobile concentra on
Flux concentra onθ M = (1 − w)(θ m,m + θ im,m)θ F = wθ f
θ = wθ f + (1 − w)(θ m,m + θ im,m)hm
hfqm(1 − w)
q fw
qm(1 − w) + qfw
cf
cim,m
cm,m

Average liquid concentra on

wc f θ f + (1− w )(c m,m θ m,m + c im,m θ im,m )w θ f + (1− w )( θ m,m + θ im,m )Total solute mass (no chemical nonequilibrium)wc f(θ f + ρ fKdf) + (1 − w){cm,m(θ m,m + fm ρ mKdm) + cim,m[θ im,m + (1 − fm)ρ mKdm]}

www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

wc f q f + (1− w )c m q m
wq f + (1− w )q m

w ⎡⎢ c f (θ f + f f ρ f K df )+ ρ f s fk ⎤⎥ +⎣⎦⎡(1− w ) ⎢ c m (θ m + f m ρ m K dm )+ ρ m s mk ⎥⎤⎣⎦wc f q f + (1− w )c m,m q m
wq f + (1− w )q m

791

F . 4. Breakthrough curves calculated using the Uniform Transport Model for a 10-cm-long soil column and the following parameters: q = 5

cm d−1, θ = 0.5, and ρb = 1.5 g cm−3; on the le , Kd = 1 cm3 g−1 and λ = 0.1, 1, and 10 cm; and on the right, λ = 1 cm and Kd = 0, 1, and 3 cm3
g−1.

F . 5. Breakthrough curves calculated using the Mobile–Immobile Water Model for a 10-cm-long soil column and the following parameters: q= 3 cm d−1, θ = 0.5, λmo = 1 cm, Kd = 1 cm3 g−1, ρb = 1.5 g cm−3. On the le : fmo = 0.6, θmo = 0.3, θim = 0.2, and ωmim = 0.1, 0.5, 10 d−1. On theright: ωmim = 0.5 d−1 and (a) fmo = 0.4, θmo = 0.2, θim = 0.3; (b) fmo = 0.6, θmo = 0.3, θim = 0.2; and (c) fmo = 0.8, θmo = 0.4, θim = 0.1.

F . 6. Breakthrough curves calculated using the Dual-Permeability

Model for a 10-cm-long soil column and the following parameters:
qm = 3 cm d−1, qf = 30 cm d−1, θ = θm = θf = 0.5, w = 0.1, λm = λf =1 cm, Kdm = Kdf = 1 cm3 g−1, ρb = 1.5 g cm−3, and ωdp = 0, 0.1, and
0.5 d−1. Matrix, fracture, and total breakthrough curves are represented by thin, medium, and thick lines, respec vely.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008792

rapid equilibration of the concentrations of the two domains. The
breakthrough curves from the two regions will then start converging and will again resemble the results of the uniform transportmodel applied to the entire soil (blue lines in Fig. 6). For a fully
equilibrated system, R equals 4 and one pore volume equals 0.88
d, causing the solute to arrive at approximately 3.5 d.

requires only one additional parameter: the mass transfer coeﬃcient αk characterizing the sorption process. Three diﬀerent values
for αk of 0.1, 0.5, and 10 d−1 were considered. The larger valueof the mass transfer coeﬃcient (αk = 10 d−1) leads to relatively
fast equilibration between the liquid- and solid-phase concentrations, and therefore to a relatively sigmoidal breakthrough curve(blue line in Fig. 8) resembling uniform transport (red line in
Fig. 4). Lower values of the mass transfer coeﬃcient lead to lesssorption because of slower equilibration between the liquid andsolid phases, and consequently to earlier solute arrival and more
pronounced tailing (black line in Fig. 8).

Dual-Permeability with Mobile–Immobile Water Model

The simulations above (red lines in Fig. 6) were further
modiﬁed by assuming the presence of immobile water (varying
between water contents of 0 and 0.4) and a corresponding fraction of sorption sites in the matrix region. The resulting model
requires three additional parameters: the immobile water content
of the matrix domain (θim,m), the fraction of sorption sites incontact with the mobile zone of the matrix domain (fm), and themass transfer coeﬃcient characterizing solute exchange between
the mobile and immobile zones of the matrix domain (αm). Since
the same ﬂux is forced through an increasingly smaller part ofthe matrix domain, the average pore water velocity increases andsolute in the matrix region will arrive earlier at the end of thecolumn (Fig. 7). Increased concentrations in the mobile zone of
the matrix domain lead to smaller gradients between the fracture
and matrix domains, and correspondingly less solute exchange
between the two domains. This, in turn, causes less tailing in the
breakthrough curve of the fracture domain. When only 20% ofwater in the matrix domain is mobile (θm = 0.5, θim,m = 0.1, θm,m
= 0.4), the average pore water velocity is only half of that in thefracture domain (green lines in Fig. 7).

Two-Site Model

The same simulations as for the One Kinetic Site Model wererepeated with the Two-Site Model, assuming that 40% of the sorption sites were in equilibrium with the liquid-phase concentrations(Fig. 9, left). This model, compared with the One Kinetic Site
Model, requires only one additional parameter: the fraction of sorption sites (fe) in equilibrium with the liquid-phase concentration.
The eﬀect on the breakthrough curve obtained with a mass transfer
coeﬃcient, αk, of 10 d−1 was relatively small since this value by
itself represents rapid sorption. Other breakthrough curves deviated less from this simulation than for the One Kinetic Site Modelbecause of the relatively high fraction of sorption sites that were atequilibrium with the liquid phase. Similarly to the One Kinetic Site
Model, smaller values of the mass transfer coeﬃcient caused earlier
solute arrival and more prolonged tailing, but to a lesser extent.
Figure 9 (right) demonstrates the eﬀect of the fraction of sorption sites in equilibrium with the liquid phase on the computed
breakthrough curves. The red lines in Fig. 9 (left and right) again
represent the same simulation. A smaller fraction of sorption sites
in equilibrium with the liquid phase leads to less instantaneous
sorption, more kinetic sorption, and more pronounced tailing, andthus correspondingly to earlier solute arrival than when a larger
fraction of instantaneous sites is considered.

Chemical Nonequilibrium Transport Models
One Kine c Site Model

The simulation of Fig. 4 (red, uniform transport) was taken
as a basis for evaluating the eﬀect of kinetic sorption. Comparedwith the Uniform Flow Model, the One Kinetic Site ModelF . 7. Breakthrough curves calculated using the Dual-Permeability

Model with MIM for a 10-cm-long soil column and the following
parameters: qm = 3 cm d−1, qf = 30 cm d−1, θ = θm = θf = 0.5, w = 0.1,λm = λf = 1 cm, Kdm = Kdf = 1 cm3 g−1, ρb = 1.5 g cm−3, ωdp = 0.1 d−1,
ωdpm = 0.1 d−1, θim,m (thim) = 0.0, 0.1, 0.3, and 0.4, and fm = 1, 0.8,
0.4, and 0.2, respec vely. Matrix, fracture, and total breakthroughcurves are represented by thin, medium, and thick lines, respec vely.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

F . 8. Breakthrough curves calculated using the One Kine c Site

Model for a 10-cm-long soil column and the following parameters:
q = 5 cm d−1, θ = 0.5, λ = 1 cm, Kd = 1 cm3 g−1, ρb = 1.5 g cm−3, andαk = 0.1, 0.5, and 10 d−1.
793

F . 9. Breakthrough curves calculated using the Two-Site Kine c Model for a 10-cm-long soil column and the following parameters: q = 5 cm

d−1, θ = 0.5, λ = 1 cm, Kd = 1 cm3 g−1, and ρb = 1.5 g cm−3; on the le , fe = 0.4 and αk = 0.1, 0.5, and 10 d−1; and on the right, αk = 0.5 d−1
and fe = 0.1, 0.4, and 0.8.

Two Kine c Sites Modelthat the red lines in Fig. 11 (left and right) again represent thesame simulation. Results presented in Fig. 11 (left) display thesame characteristics as results shown in Fig. 9 for the Two Site
Kinetic Model. Lower values of the mass transfer coeﬃcient lead
to less sorption because of slower equilibration between the liquidand solid phases, and hence to earlier solute arrival and more
pronounced tailing (black line in Fig. 11, left). A smaller fractionof sorption sites in equilibrium with the liquid phase leads to lessinstantaneous sorption and more kinetic sorption, and hence toearlier solute arrival than when a larger fraction of instantaneoussites is considered (Fig. 11, right).

The colloid attachment–detachment model has the same
number of parameters as the One Kinetic Site Model, except thatthe parameters Kd and αk are replaced with attachment (ka) anddetachment (kd) coeﬃcients. Since we considered two diﬀerent
sorption sites, we have two pairs of attachment and detachment
coeﬃcients (ka1, kd1, ka2, kd2). Here we used Eq. [19] to convert
the One-Site Kinetic Model parameters (Kd = 1 cm3 g−1, αk =0.5 d−1; red line in Fig. 8) into attachment–detachment coefﬁcients (ka1 = 1.5 d−1, kd1 = 0.5 d−1) for the ﬁrst fraction ofsorption sites. Additionally, a ﬁnite (10-d) duration of the solute
pulse was considered, while varying the attachment–detachment
coeﬃcients (ka2 and kd2, respectively) for the second fraction ofsorption sites. In the remaining examples, a ﬁnite solute pulse isused so that both increasing and decreasing limbs of the breakthrough curves can be displayed. Increased attachment to anddetachment from the second fraction of sorption sites reduced theoutﬂow concentrations during the ﬁrst part of the breakthroughcurves, but increased them during the second part (Fig. 10).

Dual-Permeability Model with Two-Site Sorp on

This model, compared with the Dual-Permeability Model,
has four additional parameters characterizing kinetic sorption inboth regions, i.e., the fractions of sorption sites in equilibriumwith the liquid-phase concentrations of both regions (fm and ff )Physical and Chemical Nonequilibrium Flow and Transport Models
Dual-Porosity Model with One Kine c Site

Compared with the Dual-Porosity Model, this model has
two additional parameters characterizing kinetic sorption in themobile zone, i.e., the fraction of sorption sites (fem) in equilibrium with the liquid-phase concentration of the mobile zone,
and the mass transfer coeﬃcient (αch). The simulations in Fig. 5
obtained with the Dual-Porosity Model were modiﬁed by assuming the presence of kinetic sorption sites in the mobile zone. We
ﬁrst assumed that 40% of the sorption sites in the mobile phase(fem = 0.4) was at equilibrium with the mobile liquid phase andthat the remaining sites were kinetic, using mass transfer coefﬁcients, αch, of 0.1, 0.5, and 10 d−1 (Fig. 11, left). We next
assumed that αch was constant and equal to 0.1 d−1 while thefraction of sorption sites in the mobile phase at equilibrium withthe liquid concentration of the mobile water varied (fem = 0.1, 0.4,
0.7, and 1) (Fig. 11, right). We further assumed in both sets ofsimulations that the duration of the solute pulse was 10 d. Noticewww.vadosezonejournal.org · Vol. 7, No. 2, May 2008

F . 10. Breakthrough curves calculated using the Two Kine c Sites

Model for a 10-cm-long soil column and the following parameters:
solute pulse dura on = 10 d, q = 5 cm d−1, θ = 0.5, λ = 1 cm, Kd = 1
cm3 g−1, ρb = 1.5 g cm−3, ka1 = 1.5 d−1, kd1 = 0.5 d−1, and ka2 = 0.0,
0.3, and 3.0 d−1 with kd2 = 0.0, 0.1, and 1.0 d−1, respec vely.
794

F . 11. Breakthrough curves calculated using the Dual-Porosity Model with One Kine c Site for a 10-cm-long soil column and the following

parameters: solute pulse dura on = 10 d, q = 3 cm d−1, θ = 0.5, θmo = 0.3, θim = 0.2, λmo = 1 cm, Kd = 1 cm3 g−1, ρb = 1.5 g cm−3, fmo = 0.6, α =0.1 d−1; on the le , fem = 0.4 and αch = 0.1, 0.5, and 10 d−1 ; and on the right, αch = 0.1 d−1 and fem = 1.0, 0.7, 0.4, and 0.1.

and the mass transfer coeﬃcients for both regions (αch,m and
αch,f ). The simulation in Fig. 6 (red line) with a mass transfer
coeﬃcient, ωdp, of 0.1 d−1 was recalculated assuming the presence of kinetic and equilibrium sorption sites in both the matrix
and fracture domains, and using a ﬁnite duration (10 d) of theapplied solute pulse (Fig. 12). Outﬂow concentrations from thematrix (thin line), the fracture domain (intermediate line), andthe entire soil (thick line) are shown in Fig. 12. As before for theOne- and Two-Site models, kinetic sorption leads to slower sorption and thus to higher concentrations and earlier solute arrivalin the eﬄuent from both regions.

information also about the water content and concentration
distributions in the soil proﬁle and between diﬀerent phases(e.g., liquid and solid) and regions (e.g., mobile and immobile,
matrix and fracture), here we limited the examples to only calculated breakthrough curves. The calculations, among otherthings, show that a single breakthrough curve will probably notprovide enough information to parameterize models involving
both physical and chemical nonequilibrium. Nkedi-Kizza et al.
(1984) previously showed that the traditionally used physicaland chemical nonequilibrium transport models (i.e., the dualporosity model with mobile and immobile ﬂow regions and thetwo-site sorption model, respectively), are mathematically identical when applied to solute breakthrough curves. By analyzing one
breakthrough curve, it is hence often not possible to discriminate which process (physical or chemical) is responsible for theobserved nonequilibrium process since both types of nonequilibrium can produce the same or very similar results. In that case,one needs additional information, such as breakthrough curves
measured simultaneously for a tracer and a reactive chemical, tobetter analyze the underlying transport processes. Interesting
applications of some of these models were presented earlier by
Pot et al. (2005) for steady-state water ﬂow conditions and by
Köhne et al. (2006) for transient-ﬂow conditions. Both studies
were performed on undisturbed soil columns and involved thetransport of a nonreactive tracer as well as a reactive pesticide[isoproturon (N,N-dimethyl-N´[4-(1-methylethyl)phenyl]urea)].
While the tracer breakthrough curves were ﬁrst used in both studies to characterize the physical conditions of ﬂow process, thepesticide breakthrough curves subsequently provided required
information for characterizing the chemical conditions.
Similarly, Bradford et al. (2004) showed that experimental
breakthrough curves of colloids or pathogenic microorganisms
can be ﬁtted equally well using diﬀerent models with diﬀerent
assumptions (e.g., attachment vs. straining), and that additional
information (in their case, the spatial distribution of colloids) isneeded to fully discriminate between the diﬀerent processes ormodels. Thorough studies are still needed to evaluate how muchand what type of information is required to fully parameterize

Discussion
Example calculations for the diﬀerent types of nonequilibrium were presented. While HYDRUS-1D provides much

F . 12. Breakthrough curves calculated using the Dual-Permeability

Model with Two-Site Sorp on for a 10-cm-long soil column and thefollowing parameters: solute pulse dura on = 10 d, qm = 3 cm d−1, qf= 30 cm d−1, θ = θm = θf = 0.5, w = 0.1, λm = λf = 1 cm, Kdm = Kdf = 1
cm3 g−1, ρb = 1.5 g cm−3, ωdp = 0.1 d−1, αch,m = αch,f = 0.1 d−1, and ff= fm = 1, 0.7, and 0.4. Matrix, fracture, and total breakthrough curves
are represented by thin, medium, and thick lines, respec vely.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

795

cim

selected models. To facilitate such studies, the objective function
for the inverse problem in HYDRUS-1D can be formulated interms of a large number of variables involving not only boundary concentration ﬂuxes, but also water and solute distributions
within the soil proﬁle and in diﬀerent phases (Table 2). For thesereasons we believe that HYDRUS-1D is a very attractive tool for
analyzing both forward and inverse ﬂow and transport problems.

liquid-phase concentration in the immobile region of thedual-porosity model and of the mobile–immobile model[M L−3]cm
liquid-phase concentration in the matrix region of the dualpermeability model [M L−3]
cm,im liquid-phase concentration in the immobile part of thematrix region of the dual-permeability model [M L−3]
cm,mo liquid-phase concentration in the mobile part of the matrix
region of the dual-permeability model [M L−3]
cmo liquid-phase concentration in the mobile region of thedual-porosity and mobile–immobile models [M L−3]
S
solute content in the liquid phase [M L−3]essorbed concentration in equilibrium with the liquid-phase
concentration [M M−1]kssorbed concentration at kinetic sorption sites [M M−1]ks1sorbed concentration at the ﬁrst fraction of kinetic sorption sites [M M−1]ks2sorbed concentration at the second fraction of kinetic sorption sites [M M−1]es im sorbed concentration in equilibrium with the liquid-phase
concentration in the immobile region of the dual-porosity
model [M M−1]es mosorbed concentration in equilibrium with the liquid-phase
concentration in the mobile region of the dual-porosity
model [M M−1]ks mosorbed concentration at kinetic sorption sites in contact withthe mobile region of the dual-porosity model [M M−1]esmsorbed concentration in equilibrium with the liquid-phase
concentration in the matrix region of the dual-permeability model [M M−1]s mksorbed concentration at kinetic sorption sites in contact
with the liquid-phase concentration in the matrix regionof the dual-permeability model [M M−1]esfsorbed concentration in equilibrium with the liquid-phase
concentration in the macropore (fracture) region of thedual-permeability model [M M−1]ksfsorbed concentration at kinetic sorption sites in contact
with the liquid-phase concentration in the macropore (fracture) region of the dual-permeability model [M M−1]

Conclusions
We summarized a wide range of nonequilibrium water ﬂow
and solute transport models that are available in the latest versionof the HYDRUS-1D software package. The models range fromclassical models simulating uniform water ﬂow and solute transport, to traditional mobile–immobile water physical and two-sitechemical nonequilibrium models, to more complex dual-permeability models that consider both physical and chemical causes of
nonequilibrium. The presented models form a hierarchical systemfrom which diﬀerent formulations can be selected for diﬀerent
applications, depending on available information and data.
The various models were divided into three groups: (i) physicalnonequilibrium transport models, (ii) chemical nonequilibriumtransport models, and (iii) physical and chemical nonequilibrium transport models. Physical nonequilibrium models include
the Mobile–Immobile Water Model, Dual-Porosity Model,
Dual-Permeability Model, and a Dual-Permeability Model with
Immobile Water. Chemical nonequilibrium models include the
One Kinetic Site Model, a Two-Site Model, and a Two-Kinetic
Sites Model. Finally, physical and chemical nonequilibrium transport models include a Dual-Porosity Model with One Kinetic
Site and a Dual-Permeability Model with Two-Site Sorption. The
fact that a large number of models has been developed over theyears is in many ways reﬂective of the extremely complicatednature of ﬁeld-scale processes in which many diﬀerent physicaland chemical processes may combine to lead to nonequilibrium
ﬂow and transport.

Appendix
Variables for Water Contents and Concentra ons

water content [L3 L−3]
water content in the macropore (fracture) region of the dualpermeability model (local pore region scale) [L3 L−3]
θF
water content in the macropore (fracture) region of thedual-permeability model (global scale) [L3 L−3]
θim water content in the immobile region of the dual-porosityand mobile–immobile models [L3 L−3]
θmwater content in the matrix region of the dual-permeability
model (local pore region scale) [L3 L−3]
θM water content in the matrix region of the dual-permeability
model (global scale) [L3 L−3]
θm,im water content in the immobile part of the matrix regionof the dual-permeability model [L3 L−3]
θm,mo water content in the mobile part of the matrix region of
the dual-permeability model [L3 L−3]
θmo water content in the mobile region of the dual-porosityand mobile–immobile models [L3 L−3]cliquid-phase concentration [M L−3]cfliquid-phase concentration in the macropore (fracture)region of the dual-permeability model [M L−3]θθfwww.vadosezonejournal.org · Vol. 7, No. 2, May 2008

A

This work is based on work supported in part by the Terrestrial Sciences Program of the Army Research Oﬃce (Terrestrial Processes andLandscape Dynamics and Terrestrial System Modeling and Model Integration), by the National Science Foundation Biocomplexity Programs
no. 04-10055 and NSF DEB 04-21530, and by SAHRA (Sustainabilityof Semi-Arid Hydrology and Riparian Areas) under the STC Programof the National Science Foundation, Agreement no. EAR-9876800.

References
Bradford, S.A., M. Bettehar, J. Šimůnek, and M.Th. van Genuchten. 2004.
Straining and attachment of colloids in physically heterogeneous porous
media. Vadose Zone J. 3:384–394.
Brooks, R.H., and A.T. Corey. 1964. Hydraulic properties of porous media. Hydrol. Paper 3. Colorado State Univ., Fort Collins.
Brusseau, H.L., R.E. Jessup, and P.S.C. Rao. 1989. Modeling the transport ofsolutes inﬂuenced by multiprocess nonequilibrium. Water Resour. Res.
25:1971–1988.
Durner, W. 1994. Hydraulic conductivity estimation for soils with heterogeneous pore structure. Water Resour. Res. 30:211–233.
Gargiulo, G., S.A. Bradford, J. Šimůnek, P. Ustohal, H. Vereecken, and E.
Klumpp. 2007. Bacteria transport and deposition under unsaturated

796

conditions: The role of the matrix grain size and the bacteria surface protein. J. Contam. Hydrol. 92:255–273.
Gargiulo, G., S.A. Bradford, J. Šimůnek, P. Ustohal, H. Vereecken, and E.
Klumpp. 2008. Bacteria transport and deposition under unsaturated conditions: The role of water content and bacteria surface hydrophobicity. Vadose Zone J. 7:406–419 (this issue).
Gerke, H.H., and M.Th. van Genuchten. 1993a. A dual-porosity model for
simulating the preferential movement of water and solutes in structured
porous media. Water Resour. Res. 29:305–319.
Gerke, H.H., and M.Th. van Genuchten. 1993b. Evaluation of a ﬁrst-order water transfer term for variably saturated dual-porosity ﬂow models. Water
Resour. Res. 29:1225–1238.
Gerke, H.H., and M.Th. van Genuchten. 1996. Macroscopic representation of
structural geometry for simulating water and solute movement in dual-porosity media. Adv. Water Resour. 19:343–357.
Hendrickx, J.M.H., and M. Flury. 2001. Uniform and preferential ﬂow mechanisms in the vadose zone. p. 149–187. In Conceptual models of ﬂow andtransport in the fractured vadose zone. National Academy Press, Washington, DC.
Jacques, D., and J. Šimůnek. 2005. User manual of the Multicomponent Variably-Saturated Flow and Transport Model HP1. Description, veriﬁcation,and examples. Version 1.0. BLG-998 Rep. SCK⋅CEN, Mol, Belgium.
Jarvis, N.J. 1994. The MACRO model (version 3.1): Technical description andsample simulations. Rep. Diss. 19. Dep. Soil Sci., Swedish Univ. Agric.
Sci., Uppsala.
Köhne, J.M., S. Köhne, and J. Šimůnek. 2006. Multi-process herbicide transport in structured soil columns: Experiment and model analysis. J. Contam. Hydrol. 85:1–32.
Köhne, J.M., B. Mohanty, J. Šimůnek, and H.H. Gerke. 2004. Numericalevaluation of a second-order water transfer term for variably saturated
dual-permeability models. Water Resour. Res. 40:W07409, doi:10.1029/
2004WR003285.
Kosugi, K. 1996. Lognormal distribution model for unsaturated soil hydraulicproperties. Water Resour. Res. 32:2697–2703.
Nkedi-Kizza, P., J.W. Biggar, H.M. Selim, M.Th. van Genuchten, and P.J. Wierenga. 1984. Equivalence of two conceptual models for describing ion
exchange during transport through an aggregated Oxisol. Water Resour.
Res. 20:1123–1130.
Nkedi-Kizza, P., J.W. Biggar, M.Th. van Genuchten, P.J. Wierenga, H.M. Selim,
J.M. Davidson, and D.R. Nielsen. 1983. Modeling tritium and chloride 36
transport through an aggregated Oxisol. Water Resour. Res. 19:691–700.
Pot, V., J. Šimůnek, P. Benoit, Y. Coquet, A. Yra, and M.-J. Martínez-Cordón.
2005. Impact of rainfall intensity on the transport of two herbicides inundisturbed grassed ﬁlter strip soil cores. J. Contam. Hydrol. 81:63–88.
Pruess, K. 1991. Tough2: A general-purpose numerical simulator for multiphaseﬂuid and heat ﬂow. Rep. LBL-29400. Lawrence Berkeley Lab., Berkeley,CA.
Selim, H.M., J.M. Davidson, and R.S. Mansell. 1976. Evaluation of a two-site
adsorption–desorption model for describing solute transport in soil. p.
444–448. In Proc. Summer Computer Simulation Conf., Washington,
DC. 12–14 July 1976. Simulation Councils, La Jolla, CA..
Selim, H.M., L. Ma, and H. Zhu. 1999. Predicting solute transport in soils:
Second-order two-site models. Soil Sci. Soc. Am. J. 63:768–777.
Šimůnek, J., and J.W. Hopmans. 2002. Parameter optimization and nonlinear
ﬁtting. p. 139–157. In J.H. Dane and G.C. Topp (ed.) Methods of soil
analysis. Part 4. Physical methods. SSSA Book Ser. 5. SSSA, Madison,
WI.
Šimůnek, J., D. Jacques, J.W. Hopmans, M. Inoue, M. Flury, and M.Th. van
Genuchten. 2002. Solute transport during variably-saturated ﬂow: Inverse
methods. p. 1435–1449. In J.H. Dane and G.C. Topp (ed.) Methods ofsoil analysis. Part 4. Physical methods. SSSA Book Ser. 5. SSSA, Madison,
WI.
Šimůnek, J., N.J. Jarvis, M.Th. van Genuchten, and A. Gärdenäs. 2003. Review
and comparison of models for describing non-equilibrium and preferential
ﬂow and transport in the vadose zone. J. Hydrol. 272:14–35.
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2005. The HYDRUS-1D
software package for simulating the one-dimensional movement of water,
heat, and multiple solutes in variably-saturated media. Version 3.0. HYDRUS Software Ser. 1. Dep. of Environmental Sciences, Univ. of California, Riverside.

www.vadosezonejournal.org · Vol. 7, No. 2, May 2008Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2008. The HYDRUS-1D
software package for simulating the one-dimensional movement of water,
heat, and multiple solutes in variably-saturated media. Version 4.0. HYDRUS Software Ser. 3. Dep. of Environmental Sciences, Univ. of California, Riverside.
Šimůnek, J., M.Th. van Genuchten, M. Šejna, N. Toride, and F.J. Leij. 1999.
The STANMOD computer software for evaluating solute transport in porous media using analytical solutions of convection–dispersion equation.
Versions 1.0 and 2.0. IGWMC-TPS-71. Int. Ground Water Modeling
Ctr., Colorado School of Mines, Golden.
Toride, N., F.J. Leij, and M.Th. van Genuchten. 1993. A comprehensive set ofanalytical solutions for nonequilibrium solute transport with ﬁrst-orderdecay and zero-order production. Water Resour. Res. 29:2167–2182.
van Dam, J.C., J. Huygen, J.G. Wesseling, R.A. Feddes, P. Kabat, P.E.V. van
Valsum, P. Groenendijk, and C.A. van Diepen. 1997. Theory of SWAP,
version 2.0. Simulation of water ﬂow, solute transport and plant growthin the soil–water–atmosphere–plant environment. Dep. of Water Resour.
Rep. 71. Tech. Doc. 45. DLO Winand Staring Centre, Wageningen, theNetherlands.
van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.
van Genuchten, M.Th. 1981. Non-equilibrium transport parameters from miscible displacement experiments. Res. Rep. 119. U.S. Salinity Lab., Riverside, CA.
van Genuchten, M.Th., and R.W. Cleary. 1979. Movement of solutes in soils:
Computer-simulated and laboratory results. p. 349–386. In G.H. Bolt(ed.) Soil chemistry: B. Physico-chemical models. Dev. in Soil Sci. 5B.
Elsevier Scientiﬁc Publ., Amsterdam.
van Genuchten, M.Th., and F.N. Dalton. 1986. Models for simulating salt
movement in aggregated ﬁeld soils. Geoderma 38:165–183.
van Genuchten, M.Th., and R.J. Wagenet. 1989. Two-site/two-region models
for pesticide transport and degradation: Theoretical development and analytical solutions. Soil Sci. Soc. Am. J. 53:1303–1310.
van Genuchten, M.Th., and P.J. Wierenga. 1976. Mass transfer studies in sorbing
porous media: I. Analytical solutions. Soil Sci. Soc. Am. J. 40:473–481.

797