Published online August 24, 2006

Colloid-Facilitated Solute Transport in Variably Saturated Porous Media: Numerical
Model and Experimental Verification
Jirka Šimůnek,* Changming He, Liping Pang, and S. A. Bradford
et al., 1988; Noell et al., 1998), pesticides (Vinten et al.,
1983; Kan and Tomson, 1990; Lindqvist and Enfield,
1992), pharmaceuticals (Tolls, 2001; Thiele-Bruhn, 2003),
hormones (Hanselman et al., 2003), and other contaminants (Magee et al., 1991; Mansfeldt et al., 2004). Since
mobile colloids often move at rates similar or fasteras nonsorbing tracers, the potential of enhanced transport of colloid-associated contaminants can be very significant (e.g., McCarthy and Zachara, 1989). Failure toaccount for colloid-facilitated solute transport can severely underestimate the transport potential and risk assessment for these contaminants. Models that canaccurately describe the various mechanisms controlling
colloid and solute transport, and their mutual interactions and interactions with the solid phase, are essential
for improving predictions of colloid-facilitated transportof solutes in variably saturated porous media.
The transport behavior of dissolved contaminant
species has been studied for many years. By comparison,colloid transport and the mutual interactions among contaminants, colloids, and porous media are less well understood. While colloids are subject to similar subsurface
fate and transport processes as chemical compounds, they
are also subject to their own unique complexities (vanGenuchten and Šimůnek, 2004). Since many colloids andmicrobes are negatively charged, they are electrostatically repelled by negatively-charged solid surfaces. This
will lead to an anion exclusion process that can cause
slightly enhanced transport relative to fluid flow. The advective transport of colloids may similarly be enhancedby size exclusion, which limits their presence to the larger
pores (Bradford et al., 2003, 2006). In addition to being
subject to adsorption–desorption process at solid surfaces,
colloids are also affected by straining in the porous matrix(Bradford et al., 2003, 2006) and may accumulate at air–water interfaces (Wan and Wilson, 1994; Thompson andYates, 1999; Wan and Tokunaga, 2002). All of these
additional complexities require colloid transport models
to be more flexible than regular solute transport models.
Models that consider colloid-facilitated solute transport are based on mass balance equations for all colloid
and contaminant species. The various colloid-facilitatedsolute transport models that have appeared in the literature differ primarily in the manner which colloidtransport and contaminant interactions are handled. For
example, Mills et al. (1991) and Dunnivant et al. (1992)
assume that colloids are nonreactive with the solid phase,
Corapcioglu and Jiang (1993) and Jiang and Corapcioglu(1993) consider a first-order kinetic attachment of colloids, Saiers and Hornberger (1996) consider an irreversible nonlinear kinetic attachment of colloids, and van deWeerd and Leijnse (1997) describe colloid attachment
kinetics using the Langmuir equation. All colloidfacilitated transport models account for interactions be-

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

ABSTRACT
Strongly sorbing chemicals (e.g., heavy metals, radionuclides,
pharmaceuticals, and explosives) in porous media are associatedpredominantly with the solid phase, which is commonly assumed to bestationary. However, recent field- and laboratory-scale observations
have shown that in the presence of mobile colloidal particles (e.g.,
microbes, humic substances, clays, and metal oxides), colloids can act
as pollutant carriers and thus provide a rapid transport pathway forstrongly sorbing contaminants. To address this problem, we developed
a one-dimensional numerical model based on the HYDRUS-1D
software package that incorporates mechanisms associated withcolloid and colloid-facilitated solute transport in variably saturated
porous media. The model accounts for transient variably saturated
water f low, and for both colloid and solute movement due toadvection, diffusion, and dispersion, as well as for solute movement
facilitated by colloid transport. The colloid transport module additionally considers the processes of attachment/detachment to/fromthe solid phase and/or the air–water interface, straining, and/or size
exclusion. Various blocking and depth dependent functions can beused to modify the attachment and straining coefficients. The solutetransport module uses the concept of two-site sorption to describe
nonequilibrium adsorption–desorption reactions to the solid phase.
The module further assumes that contaminants can be sorbed onto
surfaces of both deposited and mobile colloids, fully accounting for the
dynamics of colloid movement between different phases. Applicationof the model is demonstrated using selected experimental data frompublished saturated column experiments, conducted to investigate thetransport of Cd in the presence of Bacillus subtilis spores in alluvial
gravel aquifer media. Numerical results simulating bacteria transport,
as well as the bacteria-facilitated Cd transport, are compared withexperimental results. A sensitivity analysis of the model to various
parameters is also presented.

C

OLLOIDAL PARTICLES (e.g., humic substances, clays,
metal oxides, and microorganisms) are commonly
found in subsurface environments (Kretzschmar et al.,1999). Many contaminants can sorb onto colloids insuspension, thereby increasing their concentrations insolution beyond thermodynamic solubilities (Kim et al.,1992). Experimental evidence now exists that many contaminants are transported not only in a dissolved state
by water, but also sorbed to moving colloids. This colloidfacilitated transport has been illustrated in the literature for numerous contaminants, including heavy metals(Grolimund et al., 1996), radionuclides (Von Gunten

J. Šimůnek and C. He, Dep. of Environmental Sciences, Univ. of California Riverside, CA 92521; L. Pang, Institute of Environmental
Science & Research, Christchurch, NZ; S.A. Bradford, George E. Brown
Jr. Salinity Laboratory, USDA, ARS, Riverside, CA 92521. Received 22
Dec. 2005. *Corresponding author (Jiri.Simunek@ucr.edu).
Published in Vadose Zone Journal 5:1035–1047 (2006).
Original Researchdoi:10.2136/vzj2005.0151ª Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

1035

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1036

VADOSE ZONE J., VOL. 5, AUGUST 2006

tween the contaminants and colloids. Various equilibriumand kinetic models have been used for this purpose.
Although already relatively complex, no existingmodel for colloid-facilitated solute transport to ourknowledge includes all of the major processes contributing to colloid and colloid-facilitated solute transport. For
example, the majority of models for colloid-facilitatedsolute transport consider flow and transport only in fully
saturated groundwater systems, usually for steady-state
flow, and thus do not account for colloid interactions withthe air–water interface. Also, no colloid-facilitatedtransport model has considered straining and size
exclusion as mechanisms of colloid retention and transport, respectively. These two processes, and especiallystraining, have recently received much attention sinceclassical colloid transport models are often unable todescribe simultaneously both breakthrough curves vs.
time and concentration profiles vs. depth (Bradford et al.,2003, 2006; Li et al., 2004; Tufenkji and Elimelech, 2005).
Straining involves the entrapment of colloids in downgradient pores that are too small to allow particle
passage. The critical pore size for straining will dependon the size of the colloid and the pore-size distribution of
the medium (McDowell-Boyer et al., 1986; Bradford
et al., 2002, 2003). Straining may have significant implications for colloid-facilitated solute transport, as illustrated by Bradford et al. (2006). Considering average
capillary pressure-saturation curves for the 12 major soil
textural groups given by Carsel and Parrish (1988), theycalculated that a 2 mm colloid, which is the size of a clay
particle, will be excluded or strained in 10 to 86% of thesoil pore space for the various soil textures (Bradford
et al., 2006). These percentages should significantly
increase if the soils become unsaturated.
Size exclusion, a process closely related to straining,
affects the mobility of colloids by constraining them to flow
domains and pore networks that are physically accessible(Ryan and Elimelech, 1996; Ginn, 2002). Electrostaticforces also play an important role in the distribution (andmobility) of colloids. Anionic colloids will be excluded
from locations adjacent to negatively charged solid
surfaces; similar to the much reported anion exclusion
process for anionic solutes (Krupp, 1972; Gvirtzman and
Gorelick, 1991; Ginn, 1995). In case of size or anion
exclusion, colloids will tend to reside in larger pores and in
more conductive parts of the flow domain. As a result, colloids will be transported faster than a conservative solute
tracer (Reimus, 1995; Cumbie and McKay, 1999; Harter
et al., 2000; Bradford et al., 2004). Differences in thedispersive flux for colloids and a conservative solute tracer
are also anticipated as a result of exclusion (Scheibe andWood, 2003). Bradford et al. (2002) observed that thedispersivity of 3.2 mm carboxyl latex colloids was up toseven times greater than bromide in saturated aquifer
sand. Conversely, Sinton et al. (2000) found in a fieldmicrobial transport experiment that the apparent colloid
dispersivity decreased with increasing particle size.
Colloid and colloid-facilitated contaminant transport inpartially saturated porous media is even more complex
than in water-saturated systems. In addition to all of theprocesses and difficulties discussed above, colloid trans-

port in partially saturated porous media is further complicated by the presence of an air phase and thin water
films, in addition to the solid and water phases present insaturated media. Wan and Wilson (1994) observed thatcolloidal particles deposit preferentially on the air–waterinterface via a capillary force acting on the particles, andthat particle transport was tremendously retarded since
the air–water interface acted as a strong sorption phase.
Another physical restriction on colloid transport inunsaturated systems is imposed by thin water films; thisprocess is often referred to as film straining (Wan andTokunaga, 1997; Saiers and Lenhart, 2003). Wan andTokunaga (1997) proposed that colloid transport in unsaturated systems depends on the ratio of the colloid size
to water film thickness. Corapcioglu and Choi (1996) developed a mathematical model describing colloid transport in unsaturated porous media, and also studied theeffect of colloids on volatile contaminant transport and
air–water partitioning in unsaturated porous media.
Most current models for colloid-facilitated solutetransport assume that the number of colloids with respect to contaminant is large and that kinetic reactions
coefficients are not dependent on the number of colloids
in the system. Although this may be true for some systems, the number of colloids (or concentrations) is often
highly variable, with colloids being mobilized (or immobilized) due to changing chemical or hydrologicalconditions. Thus the reaction coefficients need to be adjusted to the number of colloids in the system in different
phases (i.e., mobile, immobile, attached to the air–waterinterface). This adjustment needs to be performed also
for numerical stability reasons. For example, if the number of colloids in the system decreases dramatically andthe sorption constants for solute to colloids are assumed
to be constant, this may lead to large sorbed concentrations, and hence numerical instabilities.
In this study we developed a one-dimensional numericalmodel based on the HYDRUS-1D software package thatincorporates processes associated with colloid and colloidfacilitated solute transport in variably saturated porous
media. The model accounts for both colloid and solute
movement due to advection, diffusion, and dispersion invariably saturated soils, as well as for solute movement
facilitated by colloid transport. The colloid transportmodule additionally considers the processes of attachment/detachment to/from the solid phase and/or the air–water interface, straining, and/or size exclusion. The modelallows for different pore water velocities and dispersivities
for colloids and solute. In this paper we first present themodel itself, provide information about the numerical
techniques used to solve the governing equations, andprovide a simple application of the model to laboratory
data. We also present a brief sensitivity analysis of themodel to selected transport and reaction parameters.

CONCEPTUAL MODEL
The conceptual model of colloid and colloid facilitated solute transport (Fig. 1) is based on several assumptions. We assume that the porous medium consistsof three phases, that is, solid phase, air, and water. Ac-

1037

www.vadosezonejournal.org

Air

Table 1. Species considered in the colloid and colloid-facilitatedsolute transport model (L is the length unit, M is the mass unit,
and n is the number of colloids).

Air-Water Interface Colloids, Γc

ψaca

kdca

kaca

Species description

Mobile Colloids, Cc

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Water
kac

kstr

ψsstr

kdc

ψs

Strained Colloids, Scstr

Attached Colloids, Scatt

Solid

Air

Units

Cc
att
Sc
str
Sc
Gc

nL
21
nM
21
nM
22
nL

C
Se

ML
21
MM

Sk

MM

Smc

Mn

Sic

Mn

Sac

Mn

23

23

21

21
21
21

Contaminant sorbed
to mobile colloids, Smc

and kinetic solute sorption to the solid phase, respectively,
while the factors c account for nonlinearity of the process.

kdca

kaca
kaac

ψm

kdac
kamc

kstr

kdmc

Dissolved
Contaminant, C
Kd

Solid

Variable

Contaminant sorbed to colloids at air-water interface, Sac

ψg

Water

Colloid species
Mobile colloids
Colloids attached to the solid phase
Strained colloids
Colloids attached to theair–water interface
Solute (contaminant) species
Dissolved contaminantContaminant instantaneously sorbed
to the solid phase
Contaminant kinetically sorbed to thesolid phase
Contaminant sorbed tomobile colloids
Contaminant sorbed toimmobile colloids
Contaminant sorbed to colloids
attached to the air–water interface

Instantaneously Sorbed
Contaminant, Se

kac

kdc

kdic

ω

kaic

ψi
Kinetically Sorbed
Contaminant, Sk

Contaminant sorbed to immobile colloids, Sic

Fig. 1. Schematic of the (top) colloid transport and (bottom) colloidfacilitated solute transport model.

cording to van de Weerd and Leijnse (1997), it is assumed that (a) only one type of colloids exists in thesystem that these colloids can be described by theirmean behavior, (b) colloids are stable and that theirmutual interactions in the liquid phase may be neglected, (c) immobile colloids do not affect the flow andtransport properties of the porous medium because of
their clogging of small pores, and (d) the transport ofcolloids is not affected by sorption of contaminants on
colloids. We further assume that colloids exist in fourstates: mobile (suspended in water), attached to thesolid phase, strained by the solid phase, and attached tothe air–water interface (Fig. 1 top, Table 1). We further
assume that colloid phase changes can be describedusing nonlinear first-order processes.
We further assume that contaminants can be dissolved in
the liquid phase, as well as sorbed instantaneously andkinetically to the solid phase, and to colloids in all of theirstates. We thus distinguish six different states of contaminants (Fig. 1 bottom; Table 1): dissolved, instantaneously or
kinetically sorbed to the solid phase, sorbed to mobile orimmobile colloids (either strained or/attached to the solid
phase), and colloids attached to the air–water interface.
The coefficients k in Fig. 1 (Table 2) represent various firstorder attachment/detachment or sorption/desorption
rates, the coefficients Kd and v represent instantaneous

MATHEMATICAL MODEL
Variably Saturated Water Flow
Transient one-dimensional variably saturated waterflow in porous media is described using the Richards
equation (Richards, 1931):
 

]u
]
]h
5
K
1 cos a
2S
[1]
]t
]x
]x
where h is the water pressure head [L], u is the volumetric water content [L3 L23], t is time [T], x is thespatial coordinate [L] (positive upward), S is the sinkterm [L3 L23 T21], a is the angle between the flow
direction and the vertical axis, and K(h) is the unsaturated hydraulic conductivity function [L T21]. The soil
hydraulic properties, that is, the retention curve u(h) andthe unsaturated hydraulic conductivity function K(h)can be evaluated using, for example, the analytical modelof van Genuchten (1980).
Table 2. Kinetic reactions considered in the colloid and colloidfacilitated solute transport model (T is the time unit).
Reaction description
Colloid species
Colloids attachment coefficient
Colloids detachment coefficient
Colloids straining coefficient
Colloids attachment coefficient to theair–water interface
Colloids detachment coefficient fromthe air–water interface
Solute (contaminant) species
Solute adsorption rate tomobile colloids
Solute desorption rate frommobile colloids
Solute adsorption rate toimmobile colloids
Solute desorption rate from
immobile colloids
Solute adsorption rate to colloids at
the air–water interface
Solute desorption rate from colloids at
the air–water interface

Variable

Units

kac
kdc
kstr
kaca

T
21
T
21
T
21
T

kdca

T

21

kamc

T

21

kdmc

T

21

kaic

T

21

kdic

T

21

kaac

T

21

kdac

T

21

21

1038

VADOSE ZONE J., VOL. 5, AUGUST 2006

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Colloid Transport
Colloid fate and transport models are commonly
based on some form of the advection–dispersion equation, but modified to account for colloid filtration(Harvey and Garabedian, 1991; Hornberger et al., 1992;
Corapcioglu and Choi, 1996; Bradford et al., 2003). A
more comprehensive equation for colloid transport describing both colloid/matrix and colloid/air–water interface mass partitioning in one-dimensional form is given by]uw Cc]Sc]Aaw Gc
]
]Cc
1r
1
5
uw Dc
]t
]t
]t
]x
]x
]qc Cc
2
1 Rc
[2]
]x
where Cc is the colloid concentration in the aqueous phase[n L23], Sc and Gc are colloid concentrations adsorbed tothe solid phase [nm21] and air–water interface [n L22],
respectively; uw is the volumetric water content accessible
to colloids [L3 L23] (due to ion or size exclusion, uw may
be smaller than the total volumetric water content u), Dc isthe dispersion coefficient for colloids [L2 T21], r is thebulk density [M L23], Aaw is the air-water interfacial area
per unit volume [L2 L23], qc is the volumetric water flux
density for colloids [L T21], while Rc represents various
chemical and biological reactions [n L23 T21]. The second
and third terms on the left side of Eq. [2] represent colloid
mass-transfer terms between the aqueous phase and thesolid phase or the air–water interface [n L23 T21],
respectively. The first two terms on the right side of Eq.
[2] represent the dispersive and advective colloid fluxes,
respectively. The value of u w is defined as
uw 5 u 2 eg
[3]where e is porosity (equal to the saturated water content,us) [L3 L23] and g is the water saturation that is not accessible to mobile colloids [-] (Bradford et al., 2006). The
volumetric water flux density for colloids, qc, moving onlyin pores from which they are not excluded, can be calculated from the ratio of the relative hydraulic conductivity of the entire pore space (Krw) to the relative hydraulicconductivity of the colloid-accessible pores (Krc)
qc 5 q

Krc
Krw

[4]

where q is the volumetric water flux density for water [LT21]. The relative hydraulic conductivity of the colloidaccessible pores, Krc [-], can then be calculated using, forexample, the pore-size distribution model of Burdine(1953) by adjusting its limits of integration to reflect waterflow to straining sites (Bradford et al., 2003) as
Sw

#
2

dSw
[h(Sw )]2g

Krc 5 Sw 1

Colloid mass transfer between the aqueous and solid
phases is traditionally described using attachment–
detachment models of the form:]Sc
5 uw cs fs kac Cc 2 rkdc Sc 2 Rscr[6]]t
where Rsc represents various chemical and biologicalreactions acting on colloids that are kinetically attached tothe matrix [n L23 T21], kac and kdc are the first-order
colloid attachment and detachment coefficients [T21],
respectively, cs is a dimensionless colloid retentionfunction [-], and fs is the fraction of the solid surface
area that is available for attachment. The attachment
coefficient is generally calculated using filtration theory(Logan et al., 1995), a quasi-empirical formulation interms of the median grain diameter of the porous medium(often termed the collector), the pore-water velocity, andcollector and collision (or sticking) efficiencies accounting
for colloid removal due to diffusion, interception andgravitational sedimentation (Rajagopalan and Tien, 1976;
Logan et al., 1995; Tufenkji and Elimelech, 2004). The
first-order detachment coefficient in Eq. [6] accounts forcolloid mobilization, presumably as affected by changes inpore-water chemistry (ionic strength, ionic composition,and pH) and physical perturbations in flow, includingchanges in the flow rate and the water content.
To simulate reductions in the attachment coefficientdue to filling of favorable sorption sites, cs is sometimes
assumed to decrease with increasing colloid mass retention. Random sequential adsorption (Johnson and Elimelech, 1995) and Langmuirian dynamics (Adamczyk et al.,1994) equations have been proposed for cs to describe thisblocking phenomenon, with the latter equation given by:Sc
cs 5 1 2 max[7]Sc
in which Scmax is the maximum solid-phase colloid concentration [n M21]. Conversely, enhanced colloid retention during porous-medium ripening can theoretically bedescribed using a functional form of cs that increases withincreasing mass of retained colloids (Tien, 1989; Tien and
Chiang, 1985; Deshpande and Shonnard, 1999) We referto several recent studies (Ginn, 2002; DeNovio et al.,
2004; Rockhold et al., 2004) for more detailed discussionsof the attachment and detachment coefficients in Eq. [6].
The first term on the right side of Eq. [6] is multiplied by
the fraction of the solid surface area, fs, that is available forattachment. Since mobile colloids are due to size exclusiontransported in only a fraction of the pore space (Sw 2 g),
only a portion of the solid surface area is accessible forattachment. The fraction of the solid surface area that isavailable for attachment can be estimated as (Bradford
et al., 2006):
Sw

Sw . gdSw#[h(S )]

2w

0

where Sw is the soil water saturation [-].

[5]

# h(S )dS
w

w

g

[8]

fs 5 1#h(S )dS
w

0

w

1039

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

Notice that Eq. [6] lumps the effects of a variety ofphysical and chemical processes into a single attachmentcoefficient parameter. The attachment coefficient in
Eq. [6] has been found to depend strongly on water content, with attachment significantly increasing as the watercontent decreases. This has been attributed to a variety
of processes, including interfacial attachment and film
straining. Bradford et al. (2002, 2003, 2006) hypothesizedthat the influence of straining and attachment mechanisms on colloid retention should be separated into two
distinct components: attachment and detachment per se,
and straining (being the entrapment of colloids in porethroats that are too small to allow passage). They modeled the influence of straining using an irreversible firstorder expression as follows:
r

]Sc
]Sstr
]Satt
5 r c 1 r c 5 uw cstrs kstr Cc]t]t]t
1 (uw cs fs kac Cc 2 rkdc Sc ) 2 Rsc

for attachment [-], and kaca and kdca are the first-order
colloid attachment and detachment coefficients to/fromthe air–water interface [T21], respectively. The interfacial
area model of Bradford and Leij (1997) can be used toestimate the air–water interfacial area Aaw as:
Aaw 5 2

Sw

1

#

#

e
Paw (S)dS 5 erw g h(S)dS
saw 1
saw S

[12]

w

where e is the porosity [L3 L23], Paw and saw are the airwater capillary pressure [M L21 T22] and surface tension[M T22], respectively; h is the pressure head [L], rw is thedensity of water [M L23], and g is the gravitational
acceleration [L T22].

Colloid-Facilitated Solute Transport[9]21

where kstr is the first-order straining coefficient [T ], and
Scstr and Scatt are the solid-phase concentrations ofstrained and attached colloids [n M21], respectively. The
first term on the right-hand side of the above equation
now accounts for straining and the second term (in parentheses) for attachment–detachment.
Application of the first-order attachment–detachment
model given by Eq. [6] typically leads to exponential
colloid distributions vs. depth. Bradford et al. (2003)
showed that such exponential distributions are often inconsistent with experimental data. They obtained muchbetter results using a depth-dependent straining coefficient cstrs in Eq. [9] of the form
d50 1 x 2b
5
[10]
cstr
s
d50
where d50 is the median grain size of the porous media[L], b is a fitting parameter [-], and x is distance from theporous-medium inlet [L]. Data from Bradford et al.
(2002, 2003) for different colloid diameters (dp) andporous media median grain sizes showed an optimalvalue of 0.43 for b, while kstr could be described using aunique increasing function of the ratio of dp/d50.
Subsequent studies with layered soils (Bradford et al.,2004) showed the importance of straining at and close to
textural interfaces when flow occurs in the direction ofcoarse-textured to medium- and fine-textured media.
Straining was found to be a significant mechanism forcolloid retention for values of dp/d50 greater than 0.005.
Finally, a model similar to Eq. [6] may be used to describe
the partitioning of colloids to the air–water interface]Aaw Gc
5 uw caca fa kaca Cc 2 Aaw kdca Gc 2 Raca
[11]]twhere Gc is the colloid concentration adsorbed to the air–water interface [n L22], Raca represents various chemicaland biological reactions of attached colloids to the air–water interface [n L23 T21], Aaw is the air–waterinterfacial area per unit volume [L2 L23], caca is adimensionless colloid retention function for the air–waterinterface [-] similarly as used in Eq. [6], and fa is thefraction of the air–water interfacial area that is available

Colloid-facilitated solute transport model requires
knowledge of colloid transport, dissolved contaminanttransport, and colloid-facilitated contaminant transport,
and of various interactions between colloids, solute, soil,
and air phase. Transport and/or mass-balance equations
must therefore be formulated for the total contaminant,for contaminant sorbed kinetically or instantaneously tothe solid phase, and for contaminant sorbed to mobile
colloids, to colloids attached to the soil solid phase, andto colloids accumulated at the air–water interface, thatis, for contaminant in all different phases or pools.
Our colloid-facilitated transport model closely follows
development of many similar models (e.g., Corapcioglu
and Jiang, 1993; Jiang and Corapcioglu, 1993; Saiers andHornberger, 1996; van de Weerd and Leijnse, 1997;
Bekhit and Hassan, 2005). Major differences in ourmodel include consideration of water contents, water
fluxes, and air–water interfacial areas that may change
in time and space (i.e., transient variably-saturated
water flow), the presence of colloids on the air–waterinterface, and adjustment of all kinetic rates to thenumber of colloids present in the system.
Mass-Balance Equation for the Total Contaminant
The combined dissolved and colloid-facilitated contaminant transport equation (in one dimension) is given by:
]uC
]Se
]Sk
]uw Cc Smc]Sc Sic]Aaw Gc Sac
1r
1
1r
1
1r
]t 
]t 
]t
]t
]t

 ]t
]
]C
]qC
]
]Cc
5
uD
2
1
uw Smc Dc
]x
]x
]x
]x
]x
]qc Cc Smc
1R
[13]
]x
where C is the dissolved-contaminant concentration in
the aqueous phase [M L23], Se and Sk are contaminantconcentrations sorbed instantaneously and kinetically,
respectively, to the solid phase [M M21], Smc, Sic, and Sac
are contaminant concentrations sorbed to mobile andimmobile (attached to solid and air–water interface)colloids [M n21], respectively; D is the dispersion
coefficient for contaminants in solution [L2 T21], q is
the volumetric water density for the contaminant [L T21],
and R represents various chemical and biological reac2

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1040

VADOSE ZONE J., VOL. 5, AUGUST 2006

tions, such as degradation and production [M L23 T21],
discussed below. Note that the left-hand side sums themass of contaminant associated with the different phases(contaminant in the liquid phase, contaminant sorbed instantaneously or kinetically to the solid phase, and contaminant sorbed to mobile colloids and immobile colloids
attached to solid phase or air–water interface), while theright-hand side considers various spatial mass fluxes(dispersion and advective transport of the dissolved
contaminant, as well as dispersion and advective transport of contaminant sorbed to mobile colloids) and sink/
source reactions.
Mass-Balance Equation for Contaminant Sorbed to
the Solid Phase
Eq. [13] invokes the concept of two-site sorption formodeling non-equilibrium adsorption–desorption reactions (e.g., van Genuchten and Wagenet, 1989). The
two-site sorption concept assumes that total sorption, S,
can be divided into two fractions:S 5 Se 1 S k

[14]

21

with sorption Se [M M ] on one fraction of the sites(type-1 sites) assumed to be instantaneous, and sorption
Sk [M M21] on the remaining sites (type-2 sites) being
time dependent according to]Sk
5 v[(1 2 f )C(C) 2 Sk ] 1 Rsk
r
[15]
]t
where v is the first-order rate constant [T21], f is thefraction of exchange sites assumed to be in equilibriumwith the solution phase [-], C(C) is the adsorption isotherm [M M21] that can be expressed using Freundlich,
Langmuir, Freundlich-Langmuir, or linear adsorption
models (Šimůnek et al., 1998), and Rsk represents various chemical and biological reactions of the kineticallysorbed contaminant [M L23 T21].
Mass-Balance Equation for Contaminant Sorbed to
Mobile Colloids
The mass-balance equation for contaminant sorbed tomobile colloids can be written as]uw Cc Smc]]Cc]qc Cc Smc
52
1 ucm kamc C
uw Smc Dc
]t
]x
]x
]x
2 uw kdmc Cc Smc 2 uw (cs fs kac 1 cstr kstr )Cc Smc
1 rkdc Sc Sic 2 uw caca fa kaca Cc Smc
1 Aaw kdca Gc Sac 1 Rmc[16]

where kamc is the adsorption rate to mobile colloids[T21], kdmc is the desorption rate from mobile colloids[T21], and Rmc represents various chemical and biological reactions for contaminant sorbed to mobile colloids[M L23 T21]. The parameter cm adjusts the sorption rate
to the number of mobile colloids present, that is,
cm 5 Cc /Ccref
ref

[17]

where Cc is the reference concentration of colloids forwhich the sorption rate kamc is valid [M L23]. In Eq. [16],the first two terms on the right-hand side represent
dispersive and advective transport, respectively, of contaminant sorbed to mobile colloids; the third and fourthterms account for sorption and desorption of contaminants to/from mobile colloids, respectively; the fifth andsixth terms account for the attachment (including straining) and detachment of mobile colloids containing sorbed
contaminants, respectively; the seventh and eighth terms
account for the attachment and detachment of mobile colloids containing sorbed contaminants to/from air–waterinterface, respectively; while the ninth term represents
degradation or other reactions involving contaminantsorbed to mobile colloids.
Mass-Balance Equation for Contaminant Sorbed to
Immobile Colloids
The mass-balance equation for contaminant sorbed toimmobile colloids can be written as follows]Sc Sicr
5 uci kaic C 2 rkdic Sc Sic]t
1 uw (cs fs kac 1 cstr kstr )Cc Smc
2 rkdc Sc Sic 1 Ric[18]

where kaic is the adsorption rate to immobile colloids[T21], kdic is the desorption rate from immobile colloids[T21], and Ric represents various reactions for contaminant sorbed to immobile colloids [M L23 T21]. The
parameter ci adjusts the sorption rate to the number ofimmobile colloids present:
ci 5 Sc /Srefc[19]

where Scref is the reference concentration of immobile
colloids for which the sorption rate kaic is valid [nm21].
In Eq. [18] the first two terms on the right-hand side
represent adsorption and desorption of contaminant to/from immobile colloids, respectively; the third andfourth terms describe the attachment (including straining) and detachment of immobile colloids with sorbed
contaminant, respectively; and the fifth term represents
reactions of contaminant sorbed to immobile colloids.
Mass-Balance Equation for Contaminant Sorbed to
Colloids Attached to the Air–Water Interface
The mass-balance equation for contaminant sorbed to
colloids attached to the air–water interface may be
written as]Aaw Gc Sac
5 ucg kaac C 2 Aaw kdac Gc Sac]t
1 uw caca fa kaca Cc Smc 2 Aaw kdca Gc Sac 1Rac [20]where kaac is the adsorption rate to colloids at the air–water interface [T21], kdac is the desorption rate fromcolloids at the air–water interface [T21], and Rac
represents various reactions for contaminant sorbed to
colloids at the air–water interface [M L23 T21]. Similarly
as in Eq. [19], the parameter cg adjusts the sorption rate
to the number of colloids at the air–water interface:cg 5 Gc /Grefc[21]

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

where Gcref is the reference concentration of immobile
colloids for which sorption rate kaac is valid [nL22]. In
Eq. [20] the first two terms on the right-hand side
represent the sorption and desorption, respectively, ofcontaminant to/from colloids at the air–water interface;
the third and fourth terms account for the attachment
and detachment, respectively, of colloids with sorbed
contaminant to/from the air–water interface; whereas
the fifth term represents degradation and other reactions of contaminant sorbed to colloids accumulated at
the air–water interface.
Reaction Term
The reaction term R in Eq. [13] may be used to account for a variety of chemical and biological reactions and transformations, including degradation andproduction, not already explicitly incorporated in themain total contaminant mass transport equation. In the
newly developed colloid-facilitated solute transportmodule we do not fully support current capabilities ofthe HYDRUS software package to simulate sequential
first-order decay chains and can simulate transport ofonly one single compound. The reaction term R may
thus include provisions for only one first-order degradation reaction that may have different rate coefficientsin each phase (i.e., in the liquid, sorbed, and colloid
phase). The reaction term R in Eq. [13] for colloidfacilitated transport scenarios is now given by:
R 5 2mw uC 2 ms r(Sk þ Se )2 mc (uw Cc Smc þ rSc Sic þ Aaw Gc Sac )[22]

where mw, ms, and mc are first-order rate constants [T21]for solutes in the liquid, solid, and colloid phases, respectively. Note that we assume that the solute degradation rate is the same for solute sorbed to all colloid,that is, mobile and immobile.
The terms Rsk, Rmc, Ric, and Rac for reactions in the
kinetically sorbed phase, on mobile colloids and on
colloids associated either with the solid phase or the air–water interface, respectively, are as follows:
Rsk 5 2ms rSk
Rmc 5 2mc uw Cc Smc
Ric 5 2mc rSc Sic
Rac 5 2mc Aaw Gc Sac[23]

SOLUTION OF THE
GOVERNING EQUATIONS
The above mathematical development shows that acomplete description of colloid-facilitated solute transport requires, in addition to variably-saturated flow
equation, a total of nine coupled partial differential
equations involving nine unknown variables (Cc, Scstr,
Scatt, Gc, C, Sk, Smc, Sic, and Sac). Four partial-differential
equations pertain to four forms of colloids, that is,mobile colloids (Cc), strained colloids (Scstr), attachedcolloids to the solid phase (Scatt), and colloids attached tothe air–water interface (Gc). Additional five partialdifferential equations pertain to five forms of contami-

1041

nants, that is, dissolved contaminant (C), contaminant
kinetically sorbed to the solid phase (Sk), contaminantsorbed to mobile colloids (Smc), contaminant sorbed to
colloids associated with the solid phase (Sic), andcontaminant sorbed to colloids attached to the air–water interface (Sac). Note that the instantaneously
sorbed contaminant does not need a special partial
differential equation, since it is described by the algebraic equation (adsorption isotherm) and treated simultaneously with the dissolved contaminant.
To find a solution for the flow Eq. [1] and the ninetransport equations (Eq. [2], [6], [9], [11], [13], [15], [16],[18], and [20]) analytically is difficult, if not impossible.
Numerical methods hence must be used to solve the flow
and transport equations. Our numerical solution was
implemented into the HYDRUS-1D software package(Šimůnek et al., 1998). The numerical solution is performed at each time step in three sequential steps.
First, because of the assumption that colloids do notaffect the flow and transport properties of the porous
medium, the variably-saturated flow Eq. [1] can besolved independently of the colloid and solute transport
equations. Had we considered the effect of attached andstrained colloids on the hydraulic conductivity (i.e.,
clogging), the soil hydraulic properties would have had
to be updated at each time step depending on thenumber of strained and attached colloids. HYDRUS-1D
uses a mass-lumped linear finite elements scheme for the
spatial discretization and a fully implicit finite difference
scheme for the temporal discretization (Šimůnek et al.,1998). The mixed form of the nonlinear Richards Eq. [1]was solved with a Picard iterative solution scheme. After
achieving convergence, the nodal values of the fluid flux
are determined from nodal values of the pressure headby applying Darcy’s law, and together with nodal valuesof the water content subsequently used as input to thenumerical solution of the colloid transport and colloidfacilitated solute transport equations.
Second, because of the assumption that the transportof colloids is not affected by sorption of contaminants on
colloids, the four partial differential equations describing colloid transport can be solved independently of
the five solute transport equations. The four governing
equations (2, 6, 9, and 11) form a linear system when thedimensionless colloid retention functions cs, csstr, andcaca are equal to one (i.e., when blocking and depthdependent phenomena are not considered). In that case,
Eq. [6], [9], and [11], none of which contains spatial
derivatives, are discretized using finite differences,similarly as for physical nonequilibrium solute transportin Šimůnek et al. (1998), and incorporated directly into
the discretized form of Eq. [2]. The Galerkin finite
element method (spatial discretization) coupled with the
Crank–Nicholson finite difference method (temporal
discretization) was used to solve the colloid transport
Eq. [2] subject to appropriate initial and boundary
conditions. When the dimensionless colloid retentionfunctions cs, csstr, and caca are not equal to one (i.e.,
when blocking or depth-dependent processes are
considered), then the system of equations is nonlinear.
For those cases we used a Picard iterative solution

VADOSE ZONE J., VOL. 5, AUGUST 2006

scheme to solve the system of nonlinear equations.
Again, after the solution of the colloid transport
equations is obtained, nodal values of colloid concentrations in the different states (e.g., mobile, strained, and/or
attached) are used as input to the numerical solution of
the colloid-facilitated solute transport equations.
Finally, a Picard iterative solution scheme was used tosolve the system of five equations describing colloidfacilitated solute transport (13, 15, 16, 18, and 20). We
used again the finite element method for spatial discretization (of 13 and 15) and finite difference method
for temporal discretization (of all five equations). After
achieving convergence, and hence the solute concentrations of the different phases, the numerical solution canproceed to the next time step, starting again with theRichards equation.

MODEL VERIFICATIONWe selected experimental data of Pang et al. (2005)
to test our colloid- and colloid-facilitated solute transport model. In Pang et al. (2005), a column experiment
was conducted to study the movement of Cd in thepresence of Bacillus subtilis spores under saturated
flow conditions. The column (18 cm long, 10 cm i.d.) was
uniformly packed with coarse gravel aquifer material(d50 5 16.69 mm, d60/d10 5 53), which gave bulk densityof 1.9 g cm23 and an effective porosity of 0.27. A solution (at fixed pH) containing approximately 4 mg L21
of Cd and 2 mg L21 of Br was first injected for about
1 h (5 pore-volumes). Bacillus subtilis spores were subsequently introduced to the column, along with Cd andBr, for another 40 min (3.4 pore-volumes). After that,the column was flushed with fresh water. A constantflow rate was applied during the entire experiment.
Untreated groundwater extracted from alluvial gravelaquifers was used as the background electrolyte.
Detailed information is given in Pang et al. (2005).
Figure 2 shows the observed breakthrough curves(BTCs) of Br, B. subtilis spores and Cd from the column
experiment. Although Br and Cd solutes were injected
simultaneously in the beginning of the experiment, the
BTC of Cd showed a slow increase of concentrations com-pared to the BTC of Br, which suggests a higher adsorption capability of Cd onto the solid matrix. After injectionof the spores, a simultaneous steep rise and flat plateau ofboth B. subtilis spores and Cd occurred. This observation
implies that Cd was co-transported with the spores.
The column experiment was simulated in three steps(Pang and Šimůnek, 2006). First, the Br breakthroughcurve was used to optimize the pore-water velocity andlongitudinal dispersivity. Next, the B. subtilis breakthrough curve was used to optimize the first-orderattachment and detachment coefficients, while the porewater velocity and longitudinal dispersivity were fixed atthe values determined for Br. Finally, the Cd breakthrough curve was analyzed to obtain the distribution
coefficient, the first-order rate constant for exchanging
between two sorption sites, and the adsorption and desorption rate coefficients to/from immobile bacteria. Pangand Šimůnek (2006) independently determined adsorption and desorption rate coefficients to/from mobile
bacteria from a batch kinetic study in the absence ofaquifer media. Results of preliminary simulations suggestthat almost all Cd sorption sites of the solid phase werekinetic. These rate coefficients, together with parameters
for the bacteria and flow, were fixed during simulationsof the Cd data. Parameters were optimized using theinternal HYDRUS optimization routine (Marquardt,
1963) for the Br and bacteria data and the PEST software(Doherty, 1994) for the Cd data. The optimized parameters are given in Table 3, while a comparison of modelpredicted results and observed experimental data ispresented in Fig. 3. Notice an excellent correspondence
between optimized and measured breakthrough curves
for all three compounds. Since the column experiment was
conducted under saturated condition, the simulation did
not consider the attachment-detachment of colloids to/from the air–water interface. As we assumed the samepore-water velocity and dispersion for Br and thebacteria, size exclusion was also not considered here.
Pang and Šimůnek (2006) provide more details about
Table 3. Transport and reaction parameters used to simulate theBr, Bacillus subtilis spores, and Cd breakthrough curves presented in Fig. 3.
Description

Parameter

Value

Br tracer parameters

1.2

Concentration (C/C0)

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1042

Observed Br
Observed B.subtilis

1.0

Observed Cd
0.8
0.6

Saturated water content
us21
Ks
Saturated hydraulic conductivity, cm h23rBulk density, g cm
Longitudinal dispersivity, cmaL
Saturation
Sw
Colloid (Bacillus subtilis spores) parameters
21

First-order colloid attachment coefficient, h21First-order colloid detachment coefficient, hWater saturation that is not accessible to mobile
colloids
Solute (Cd) parameters

0.4
0.2

3

0.0
0.0

0.5

1.0

1.5

2.0

2.5

3.0

Time (hour)Fig. 2. Observed Br, Bacillus subtilis spores, and Cd breakthroughcurves for the column experiment of Pang et al. (2005).

21

Adsorption isotherm coefficient of Cd, cm g
Cd adsorption rate to mobile colloids, h21
21
Cd desorption rate from mobile colloids, h
21
Cd adsorption rate to immobile colloids, h
21
Cd desorption rate from immobile colloids, h21
Cd first-order rate constant for kinetic sorption, h

0.27
16.7
1.9
0.9
1.0

kac
kdc
g

1.697
0.0001
0

Kd
kamc
kdmc
kaic
kdic
v

2.29
2.30†
1.28†
18.7
0.00
0.25

† Independently determined from a batch kinetic study in the absence ofaquifer media (Pang and Šimůnek, 2006).

1043

www.vadosezonejournal.org

Concentration (C/C0)

a.

1.0

Observed Br
Simulated Br

0.8
0.6
0.4
0.2
0.0

Concentration (C/C0)

0.0

0.5

1.0

1.5

2.5

3.0

Observed B.subtilis

b.

0.8

2.0

Simulated B.subtilis

0.6
0.4
0.2
0.0
0.0
0.6

Concentration (C/C0)

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1.2

0.5

1.0

1.5

2.0

c.

2.5

3.0

Observed Cd

0.5

Simulated Cd
0.4
0.3
0.2
0.1
0.0
0.0

0.5

1.0

1.5

2.0

2.5

3.0

Time (hour)Fig. 3. (a) Observed and simulated Br, (b) Bacillus subtilis spores,and (c) Cd breakthrough curves for the column experiment of Panget al. (2005).

the experimental set up, batch studies, as well as
more examples of model experimental validation involving Cd transport with either B. subtilis spores orEscherichia coli.

In the sensitivity analysis here we assume that allphysical properties, such as porosity, bulk density, water
flux, dispersivity, and pore-water velocity are the same as
for the column experiment discussed above for B. subtilis
spores and Cd. Since the sensitivity of colloid BTCs tovarious parameters and processes, such as straining andsize exclusion, was presented elsewhere (e.g., Bradford
et al., 2003), we evaluate here only the sensitivity ofmodel output to the various solute transport parameters.
In this section we evaluate the effects of the solute
sorption coefficients on the solute BTC. The analysis
considers the same properties for colloids as before(Table 4). We assume that the column is fully saturated
and that all pores are accessible to colloids. We further
assume that solute sorption to almost all sorption sites ofthe solid phase is kinetic ( f 5 0.01) and that the
distribution coefficient Kd is equal to 2.5 cm3g21 (and 0,
1, and 2.5 in the sensitivity analysis). The adsorption rateof contaminant to mobile colloids is considered to beequal to the adsorption rate to immobile colloids (i.e.,
kamc 5 kaic) and, similarly, the desorption rate frommobile colloids is assumed to be equal to the desorption
rate from immobile colloids (i.e., kdmc 5 kdic). While
colloids were applied during the time interval between
1 and 1.9 h, solute was applied between 0 and 1.9 h.
In the first analysis we assumed that both theadsorption and desorption coefficients were constant(i.e., kamc 5 kdmc 5 10 h21) and varied the solute
distribution coefficient Kd from 0 (no retardation) to 1,
and 2.5 cm3g21. In the second analysis we kept constant
Kd (5 2.5 cm3g21) and kdmc (5 1 h21) and changed theadsorption coefficient kamc from 1 to 10, and 50 h21.
Finally, in the last analysis, we kept constant Kd (5 2.5 cm3g21) and kamc (5 10 h21) values, but varied the value of
kdmc from 0.1 to 1, and 10 h21.
The colloid BTCs were the same in all three runs
(Fig. 4). Figure 4 (top) shows additionally solute BTCs
for three different values of Kd. Recall that almost all
sorption sites are kinetic and thus that the distribution
coefficient does not lead to a retardation of the solute
Table 4. Parameter values used in the sensitivity analysis of colloid
and solute breakthrough curves evaluating effects of various
solute sorption coefficients. Parameter values in parenthesis
were used in alternative runs in the sensitivity analyses.
Description

Parameter

Value

Colloid parameters

SENSITIVITY ANALYSIS
A sensitivity analysis can be one of the most powerfuluses of complex transport models. Such an analysis
permits a systematic evaluation of the impact of various
model parameters (both transport and reactive), andinitial and boundary conditions on the model output.
The results of a sensitivity analysis provide insight into
the relative importance of individual reactions within acomplex geochemical system. These results can thus
help one to identify the most important parameters andprocesses, and thus help guide allocation of resources
for laboratory and field investigations (Šimůnek andValocchi, 2002).

kac
First-order colloid attachment
coefficient, h21
kdc
First-order colloid detachment
coefficient, h21
Water saturation that is not accessiblegto mobile colloids
Water saturation
Sw
Solute parameters
3 21

Distribution coefficient, cm g
Mass transfer coefficient, h21
Fraction of exchange sites assumed
to be in equilibrium with thesolution phase
First-order attachment coefficient to21
mobile and immobile colloids, hFirst-order detachment coefficient frommobile and immobile colloids, h21

1.970.00010
1Kdv
f

2.5 (0, 1, 2.5)
0.3
0.01

kamc, kaic

10 (1, 10, 50)

kdmc, kdic

10 (0.1, 1, 10)1044

VADOSE ZONE J., VOL. 5, AUGUST 2006

Colloid
Solute (10, 10, 0)
Solute (10, 10, 1)
Solute (10, 10, 2.5)

Concentration (C/C0)

1.0
0.8
0.6
0.4
0.2
0.0
0.0

0.5

1.0

1.5

2.0

2.5

3.0

Colloid
Solute ( 1,
Solute (10,Solute (50,

1, 2.5)1, 2.5)1, 2.5)Time (Hour)

Concentration (C/C0)

1.2
1.0
0.8
0.6
0.4
0.2
0.0
0.0

0.5

1.0

1.5

2.0

2.5

3.0

Time (Hour)
Colloid
Solute (10, 0.1, 2.5)
Solute (10, 1, 2.5)
Solute (10, 10, 2.5)

1.2

Concentration (C/C0)

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1.2

1.0
0.8
0.6
0.4
0.2
0.0
0.0

0.5

1.0

1.5

2.0

2.5

3.0

Time (Hour)Fig. 4. Sensitivity of simulated breakthrough curves to the (top)distribution coefficient Kd; (middle) adsorption coefficients to mobile and immobile colloids kamc and kaic, respectively; and (bottom)
desorption coefficients from mobile and immobile colloids kdmc andkdic, respectively. The first position in parentheses in the legend
represents kamc and kaic, the second kdmc and kdic, and the third Kd.

front. Kinetic adsorption to the solid phase leads to a
gradual increase in solute concentration after arrival of
the solute front, and an additional gradual increase inconcentration after arrival of the colloid front. For thecase of no adsorption to the solid phase (i.e., Kd 5 0),
solute concentrations quickly reached unit relative concentration after arrival of the solute front, and stayedat this level until the end of the solute pulse. Colloids
did not have any effect on the solute breakthroughfor this set of solute transport parameters since solute
was not retarded and thus its transport cannot be ac-celerated by colloids unless they are excluded fromsome part of the pore space and travel faster than unretarded solute. The effects of colloids on solute transport become apparent when solute is kinetically sorbed.
Increasing the kinetic sorption coefficient of the solute
produced a more pronounced effect of colloids on solute transport (Fig. 4, top).
A different picture is obtained when sorption to thesolid phase is assumed constant (Kd 5 2.5 cm3g21), but
the adsorption coefficient to the colloids (kamc 5 kaic 5 1,
10, or 50 h21) is varied. The initial phase of the soluteconcentration BTC was the same for all three runs (Fig. 4,
middle) since all solute transport parameters characterizing interactions with the soil were the same. Dramaticchanges occurred only after the arrival of colloids since
sorption rates to the colloids were different for the threeruns. The higher the adsorption rate to the colloids, thehigher the peaks of the solute concentration BTC as aconsequence of more solute arriving at the end of thecolumn sorbed to mobile colloids.
Similar results were obtained when varying the desorption rates from the mobile and the immobile colloids (kdmc 5 kdic 5 0.1, 1, or 10 h21). We obtained again
the same initial phase of the solute concentration BTC
for all three runs (Fig. 4, bottom). Then due to relatively
high sorption to both mobile and immobile colloids,
concentrations relatively quickly increased with arriving
colloids. Due to different desorption rates from the colloids, the plateaus of the solute concentration BTCs
were at different levels.

SUMMARY
We developed a one-dimensional numerical model
based on the HYDRUS-1D software package that incorporates various processes associated with colloid andcolloid-facilitated solute transport in variably saturated
porous media. The numerical model accounts for transient variably saturated water flow and for both colloid
and solute movement due to convection, diffusion, anddispersion, as well as for solute movement facilitated bycolloid transport. The colloid transport module additionally considers the processes of attachment/detachment to/from the solid phase and/or the air–water interface. A
unique feature of our colloid-facilitated solute transportmodel is that it accounts for size exclusion and depthdependent straining.
The solute transport module uses the concept of twosite sorption to describe nonequilibrium adsorption–
desorption reactions to the solid phase. The module
further assumes that contaminants can be sorbed onto
surfaces of both deposited and mobile colloids, thus fully
accounting for the dynamics of colloid transfer between
different phases. While a majority of models for colloidfacilitated solute transport assume that the number ofcolloids with respect to the contaminant is large and thatkinetic reaction coefficients are not dependent on thenumber of colloids in the system, our newly developedmodel accounts for a variable number of colloids in thesystem; this since colloids may be mobilized (or immobilized) due to changing chemical or hydrological con-

1045

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

ditions. Our model can also account for different transport
velocities of colloids and solute due to the size exclusion ofcolloids from small pores.
The use of the model was demonstrated using
experimental data from a saturated column experiment,
conducted to investigate the transport of Cd in thepresence of Bacillus subtilis spores. We also presented abrief sensitivity analysis of the model to selected reaction
parameters characterizing various sorption and desorption coefficients.
Our future work will involve further verification of
the code using experimental data derived for more
complex conditions, such as transient water flow withmobilization and immobilization of colloids, with colloids interacting with the air–water interface and/or
being excluded from part of the pore space.
APPENDIX: LIST OF VARIABLES
Aaw
C
Cc
Ccref
d50
D
Dc
f
fa
fs
g
h
kac
kaac
kaca
kaic
kamc
kdc
kdac
kdca
kdic
kdmc
kstr
K
Krc
Krw
Paw
R
Raca
Rac
Rc

air–water interfacial area per unit volume [L2 L23]
dissolved-contaminant concentration in the aqueous
phase [M L23]
colloid concentration in the aqueous phase [n L23]
reference concentration of colloids for which the sorption rate kamc is valid [M L23]
median grain size of the porous media [L]dispersion coefficient for contaminants in solution[L2 T21]dispersion coefficient for colloids [L2 T21]
fraction of exchange sites assumed to be in equilibriumwith the solution phase [-]fraction of the air–water interfacial area that is available for attachment [-]fraction of the solid surface area that is available forattachment [-]gravitational acceleration [L T22]
pressure head [L]first-order colloid attachment coefficient [T21]
adsorption rate to colloids at the air–water interface[T21]
first-order colloid attachment coefficient to the air–water interface [T21]
adsorption rate to immobile colloids [T21]
adsorption rate to mobile colloids [T21]first-order colloid detachment coefficients [T21]
desorption rate from colloids at the air–water interface[T21]
first-order colloid detachment coefficient from theair–water interface [T21]
desorption rate from immobile colloids [T21]
desorption rate from mobile colloids [T21]first-order straining coefficient [T21]
unsaturated hydraulic conductivity [LT21]
relative hydraulic conductivity of the colloid-accessible
pores [-]
relative hydraulic conductivity of the entire pore space [-]air–water capillary pressure [M L21 T22]
various chemical and biological reactions, such as degradation and production, for contaminant [M L23 T21]
various chemical and biological reactions of attachedcolloids to the air–water interface [n L23 T21]
various reactions for contaminant sorbed to colloids
at the air–water interface [M L23 T21]
various chemical and biological reactions for colloids[n L23 T21]

Ric
Rmc
Rsc
Rsk
q
qc
S
Sac
Sc
Scref
Scatt
Scmax
Scstr
Se
Sic
Sk
Smc
Sw
t
x
b
g
Gc
Gcref
e
mc
ms
mw
caca
cg
ci
cm
cs
C
v
saw
u
uw
r
rw

various reactions for contaminant sorbed to immobile
colloids [M L23 T21]
various chemical and biological reactions for contaminant sorbed to mobile colloids [M L23 T21]
various chemical and biological reactions acting on kinetically attached colloids to the matrix [n L23 T21]
various chemical and biological reactions of the
kinetically sorbed contaminant [M L23 T21]
volumetric water density for the contaminant [L T21]
volumetric water flux density for colloids [L T21]
sink term [L3 L23 T21]
contaminant concentration sorbed to immobile (attached to the air–water interface) colloids [M n21]
colloid concentrations adsorbed to the solid phase [n M21]
reference concentration of immobile colloids for which
the sorption rate kaic is valid [n M21]
solid-phase concentrations of attached colloids [n M21]
maximum solid-phase colloid concentration [n M21]
solid-phase concentration of strained colloids [n M21]
contaminant concentration sorbed instantaneously tothe solid phase [M M21]
contaminant concentration sorbed to immobile (attached to the solid phase) colloids [M N21]
contaminant concentration sorbed kinetically to thesolid phase [M M21]
contaminant concentration sorbed to mobile colloids[M n21]soil water saturation [-]time [T]spatial coordinate [L] (positive upward)
fitting parameter [-]water saturation that is not accessible to mobile colloids [-]colloid concentrations adsorbed to the air–water interface [n L22]
reference concentration of immobile colloids for which
sorption rate kaac is valid [n L22]
porosity [L3 L23]
first-order rate constant for solutes in the colloid phase[T21]
first-order rate constant for solutes in the solid phase[T21]
first-order rate constant for solutes in the liquid phase[T21]
dimensionless colloid retention function for the air–water interface [-]parameter adjusting the sorption rate to the number ofcolloids at the air–water interface [-]parameter adjusting the sorption rate to the number ofimmobile colloids present
parameter adjusting the sorption rate to the number ofmobile colloids present
dimensionless colloid retention function [-](C) adsorption isotherm [M M21]
first-order rate constant [T21]air–water surface tension [M T22]
volumetric water content [L3 L23]
volumetric water content accessible to colloids [L3 L23]
bulk density [M L23]
density of water [M L23]
ACKNOWLEDGMENTS

The work was supported by the Terrestrial Sciences Program of the Army Research Office (Terrestrial Processes andLandscape Dynamics And Terrestrial System Modeling andModel Integration), project number 45690EV.

1046

VADOSE ZONE J., VOL. 5, AUGUST 2006

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

REFERENCES
Adamczyk, Z., B. Siwek, M. Zembala, and P. Belouschek. 1994.
Kinetics of localized adsorption of colloid particles. Adv. Colloid
Interface Sci. 48:151–280.
Bekhit, H.M., and A.E. Hassan. 2005. Two-dimensional modeling ofcontaminant transport in porous media in the presence of colloids.
Adv. Water Resourc. 28:1320–1335.
Bradford, S.A., M. Bettahar, J. Šimůnek, and M.Th. van Genuchten.
2004. Straining and attachment of colloids in physically heterogeneous porous media. Vadose Zone J. 3:384–394.
Bradford, S.A., and F.J. Leij. 1997. Estimating interfacial areas for
multi-fluid soil systems. J. Contam. Hydrol. 27:83–105.
Bradford, S.A., J. Šimůnek, M. Bettahar, M.Th. van Genuchten, and S.R.
Yates. 2003. Modeling colloid attachment, straining, and exclusion
in saturated porous media. Environ. Sci. Technol. 37:2242–2250.
Bradford, S.A., J. Šimůnek, M. Bettahar, M.Th. van Genuchten, and
S.R. Yates. 2006. Significance of straining in colloid deposition:
Evidence and implications, Water Resour. Res. (in review).
Bradford, S.A., S.R. Yates, M. Bettahar, and J. Šimůnek. 2002. Physicalfactors affecting the transport and fate of colloids in saturated
porous media. Water Resour. Res. 38(12):1327 doi:10.1029/
2002WR001340.
Burdine, N.T. 1953. Relative permeability calculations from pore-size
distributiondata. Petr. Trans. Am. Inst. MiningMetall.Eng. 198: 71–77.
Carsel, R.F., and R.S. Parrish. 1988. Developing joint probability
distributions of soil water retention characteristics. Water Resour.
Res. 24:755–769.
Corapcioglu, M.Y., and H. Choi. 1996. Modeling colloid transport inunsaturated porous media and validation with laboratory column
data. Water Resour. Res. 32:3437–3449.
Corapcioglu, M.Y., and S. Jiang. 1993. Colloid-facilitated groundwater contaminant transport. Water Resour. Res. 29(7):2215–2226.
Cumbie, D.H., and L.D. McKay. 1999. Influence of diameter on particle
transport in a fractured shale saprolite. J. Contam. Hydrol. 37:139–157.
DeNovio, N.M., J.E. Saiers, and J.N. Ryan. 2004. Colloid movement in
unsaturated porous media: Recent advances and future directions.
Vadose Zone J. 3:338–351.
Deshpande, P.A., and D.R. Shonnard. 1999. Modeling the effects ofsystematic variation in ionic strength on the attachment kinetics ofPseudomonas fluorescens UPER-I in saturated sandy media. Water
Resour. Res. 35(5):1619–1627.
Doherty, J. 1994. PEST. Watermark Computing, Corinda, Australia.
Dunnivant, F.M., P.M. Jardine, D.L. Taylor, and J.F. McCarthy. 1992.
Cotransport of cadmium and hexachlorobiphenyl by dissolved
organic carbon through columns containing aquifer material.
Environ. Sci. Technol. 26:360–368.
Ginn, T. R. 1995. A brief review of bacterial transport in natural
porous media. PNL-10876, UC-402. Pacific Northwest Natl. Lab.,
Richland, WA.
Ginn, T.R. 2002. A travel time approach to exclusion on transport in porous media. Water Resour. Res. 38:1041 doi:10.1029/2001WR000865.
Grolimund, D., M. Borkovec, K. Barmettler, and H. Sticher. 1996.
Colloid-facilitated transport of strongly sorbing contaminants innatural porous media; a laboratory column study. Environ. Sci.
Technol. 30:3118–3123.
Gvirtzman, H., and S.M. Gorelick. 1991. Dispersion and advection in
unsaturated porous media enhanced by anion exclusion. Nature(London) 352:793–795.
Hanselman, T.A., D.A. Graetz, and A.C. Wilkie. 2003. Manure-borne
estrogens as potential environmental contaminants: A review.
Environ. Sci. Technol. 37:5471–5478.
Harter, T., S. Wagner, and E.R. Atwill. 2000. Colloid transport andfiltration of Cryptosporidium parvum in sandy soils and aquifer
sediments. Environ. Sci. Technol. 34:62–70.
Harvey, R.W., and S.P. Garabedian. 1991. Use of colloid filtration
theory in modeling movement of bacteria through a contaminatedsandy aquifer. Environ. Sci. Technol. 25:178–185.
Hornberger, G.M., A.L. Mills, and J.S. Herman. 1992. Bacterialtransport in porous media: Evaluation of a model using laboratory
observations. Water Resour. Res. 28:915–938.
Jiang, S., and M.Y. Corapcioglu. 1993. A hybrid equilibrium model forsolute transport in porous medium in the presence of colloids.
Colloids Surf. A. 73:275–286.

Johnson, P.R., and M. Elimelech. 1995. Dynamics of colloid depositionin porous media: Blocking based on random sequential adsorption.
Langmuir 11:801–812.
Kan, A.T., and M.B. Tomson. 1990. Ground water transport ofhydrophobic organic compounds in the presence of dissolved
organic matter. Environ. Toxicol. Chem. 9:253–263.
Kim, J.I., P. Zeh, and B. Delakowitz. 1992. Chemical interactions ofactinide ions with groundwater colloids in Gorleben aquifer
systems. Radiochim. Acta 58–9:147–154.
Kretzschmar, R., M. Borkovec, D. Grolimund, and M. Elimelech.
1999. Mobile subsurface colloids and their role in contaminanttransport. Adv. Water Resour. 66:121–193.
Krupp, H.K., J.W. Biggar, and D.R. Nielsen. 1972. Relative flow rates
of salt and water in soil. Soil Sci. Soc. Am. Proc. 36:412–417.
Li, X., T.D. Scheibe, and W.P. Johnson. 2004. Apparent decreases incolloid deposition rate coefficients with distance of transport underunfavorable deposition conditions: A general phenomenon. Environ. Sci. Technol. 38:5616–5625.
Lindqvist, R., and C.G. Enfield. 1992. Biosorption of dichlorodiphenyltri-chloroethane and hexachlorobenzene in groundwater and its
implications for facilitated transport. Appl. Environ. Microbiol.
58:2211–2218.
Logan, B.E., D.G. Jewett, R.G. Arnold, E.J. Bouwer, and C.R.
O’Melia. 1995. Clarification of clean-bed filtration models. J.
Environ. Eng. 121:869–873.
Magee, B.R., L.W. Lion, and A.T. Lemley. 1991. Transport of dissolved
organic macromolecules and their effect on the transport ofphenanthrene in porous media. Environ. Sci. Technol. 25:323–331.
Mansfeldt, T., H. Leyer, K. Barmettler, and R. Kretzschmar. 2004.
Cyanide leaching from soil developed from coking plant purifier
waste as influenced by citrate. Vadose Zone J. 3:471–479.
Marquardt, D.W. 1963. An algorithm for least-squares estimation of
nonlinear parameters. SIAM J. Appl. Math. 11:431–441.
McCarthy, J.F., and J.M. Zachara. 1989. Subsurface transport ofcontaminants. Environ. Sci. Technol. 23(5):496–502.
McDowell-Boyer, L.M., J.R. Hunt, and N. Sitar. 1986. Particle
transport through porous media. Water Resour. Res. 22:1901–1921.
Mills, W.M., S. Liu, and F.K. Fong. 1991. Literature review and model(COMET) for colloid/metals transport in porous media. Ground
Water 29:199–208.
Noell, A.L., J.L. Thompson, M.Y. Corapcioglu, and I.R. Triay. 1998. The
role of silica colloids on facilitated cesium transport through glass
bead columns and modeling. J. Contam. Hydrol. 31(1–2):23–56.
Pang, L., M. Close, M. Noonan, M. Flintoft, and P. van den Brink.
2005. A laboratory study of bacteria-facilitated cadmium transport
in alluvial gravel aquifer media. J. Environ. Qual. 34:237–247.
Pang, L., and J. Šimůnek. 2006. Evaluation of bacteria-facilitatedsolute transport model. Water Resour. Res. (in press).
Rajagopalan, R., and C. Tien. 1976. Trajectory analysis of deep-bed
filtration with the sphere-in-cell porous media model. AIChE J.
22:523–533.
Reimus, P.W. 1995. The use of synthetic colloids in tracer transport
experiments in saturated rock fractures, LA-13004-T. Los AlamosNatl. Lab., Los Alamos, NM.
Richards, L.A. 1931. Capillary conduction of fluid through porous
mediums. Physics 1:318–333.
Rockhold, M.L., R.R. Yarwood, and J.S. Selker. 2004. Coupledmicrobial and transport processes in soils. Vadose Zone J. 3:368–383.
Ryan, J.N., and M. Elimelech. 1996. Colloid mobilization and transport ingroundwater. Colloids Surf. A: Physicochem. Eng. Aspects 107:1–56.
Saiers, J.E., and G.M. Hornberger. 1996. The role of colloidal kaolinite
in the transport of cesium through laboratory sand columns. Water
Resour. Res. 32(1):33–41.
Saiers, J.E., and J.J. Lenhart. 2003. Ionic-strength effects on colloidtransport and interfacial reactions in partially saturated porous
media. Water Resour. Res. 39:1256 doi:10.1029/2002WR001887.
Scheibe, T.D., and B.D. Wood. 2003. A particle-based model of size oranion exclusion with application to microbial transport in porous
media. Water Resour. Res. 39(4):1080 doi:10.1029/200IWRO01223.
Šimůnek, J., M. Šejna, and M.Th. van Genuchten. 1998. The
HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variablysaturated media, Version 2.0. IGWMC- TPS- 70. Int., Ground
Water Modeling Center, Colorado School of Mines, Golden.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

Šimůnek, J., and A.J. Valocchi. 2002. Geochemical transport.
p. 1511–1536 In J.H. Dane and G.C. Topp (ed.) Methods of soil
analysis. Part 1. Physical methods. 3rd ed. SSSA, Madison, WI.
Sinton, L.W., M.J. Noonan, R.K. Finlay, L. Pang, and M.E. Close.
2000. Transport and attenuation of bacteria and bacteriophages
in an alluvial gravel aquifer. N. Z. J. Mar. Freshwater Res. 34(1):
175–186.
Thiele-Bruhn, S. 2003. Pharmaceutical antibiotic compounds in soils–
A review. J. Plant Nutr. Soil Sci. 166:145–167.
Thompson, S.S., and M.V. Yates. 1999. Bacteriophage inactivation at
the air-water-solid interface in dynamic batch systems. Appl.
Environ. Microbiol. 65:1186–1190.
Tien, C. 1989. Granular filtration of aerosols and hydrosols. Butterworth, London.
Tien, C., and H. Chiang. 1985. Dynamics of deep-bed filtration. II.
Experiment. AIChE J. 31(8):1360–1371.
Tolls, J. 2001. Sorption of veterinary pharmaceuticals in soils: A review.
Environ. Sci. Technol. 35:3397–3406.
Tufenkji, N., and M. Elimelech. 2004. Correlation equation forpredicting single-collector efficiency in physiochemical filtrationin saturated porous media. Environ. Sci. Technol. 38:529–536.
Tufenkji, N., and M. Elimelech. 2005. Breakdown of colloid filtration
theory: Role of the secondary energy minimum and surface charge
heterogeneities. Langmuir 21:841–852.
van de Weerd, H., and A. Leijnse. 1997. Assessment of the effect ofkinetics on colloid facilitated radionuclide transport in porous
media. J. Contam. Hydrol. 26:245–256.

1047

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., and J. Šimůnek. 2004. Integrated modeling ofvadose zone flow and transport processes. p. 37–69, x–xi. In R.A.
Feddes et al. (ed.) Proc. Unsaturated Zone Modelling: Progress,
Challenges and Applications, Wageningen UR Frontis Ser., Vol., 6,
Wageningen, the Netherlands. 3–5 Oct. 2004. Kluwer Academic
Publ., Dordrecht, the Netherlands.
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.
Vinten, A.J.A., B. Yaron, and P.H. Nye. 1983. Vertical transport ofpesticides into soil when adsorbed on suspended particles. J. Agric.
Food Chem. 31:662–664.
Von Gunten, H.R., V.E. Weber, and U. Krahenbuhl. 1988. The reactor
accident at Chernobyl: A possibility to test colloid-controlled transportof radionuclides in a shallow aquifer. J. Contam. Hydrol. 2:237–247.
Wan, J., and T.K. Tokunaga. 1997. Film straining of colloids in
unsaturated porous media: Conceptual model and experimental
testing. Environ. Sci. Technol. 31:2413–2420.
Wan, J.M., and T.K. Tokunaga. 2002. Partitioning of clay colloids at airwater interfaces. J. Colloid Interface Sci. 247:54–61.
Wan, J., and J.L. Wilson. 1994. Visualization of the role of the gaswater interface on the fate and transport of colloids in porous
media. Water Resour. Res. 30:11–23.