WATER RESOURCES RESEARCH, VOL. 47, W10503, doi:10.1029/2011WR010812, 2011

Modeling colloid transport and retention in saturated porous
media under unfavorable attachment conditions
Scott A. Bradford,1 Saeed Torkzaban,2 and Jiri Simunek2
Received 18 April 2011; revised 15 August 2011; accepted 20 August 2011; published 6 October 2011.

[1] A mathematical model is presented for colloid transport and retention in saturated

porous media under unfavorable attachment conditions. The model accounts for colloidtransport in the bulk aqueous phase and adjacent to the solid surface, and rates of colloid
collision, interaction, release, and immobilization on the solid phase. Model parameters
were estimated using (1) filtration theory; (2) calculated interaction energies in conjunction
with the Maxwellian kinetic energy model of diffusion; (3) information about the velocity
magnitude and distribution adjacent to the solid phase that was obtained from pore scalewater flow simulations; (4) colloid and collector sizes ; (5) the balance of appliedhydrodynamic and resisting adhesive torques; and (6) time dependent filling of retention
locations using a Langmuirian approach. The presented theory constrains the model
parameters and output to physically realistic values in many instances, and minimizes theneed for parameter optimization. Example simulations demonstrate that our modelingformulation is qualitatively consistent with observed trends for retention with colloid sizeand concentration, grain size, and velocity for many systems. The model provides a clear
conceptual explanation for the causes of hyperexponential, exponential, uniform, andnonmonotonic retention profiles without invoking hypotheses with regard to colloid
heterogeneity, aggregation, or multiple deposition rates. Furthermore, the model
formulation and research presented herein helps to identify areas where additional researchand theory development are still needed.

Citation: Bradford, S. A., S. Torkzaban, and J. Simunek (2011), Modeling colloid transport and retention in saturated porous media
under unfavorable attachment conditions, Water Resour. Res., 47, W10503, doi:10.1029/2011WR010812.

1.

Introduction[2] Unfavorable attachment conditions occur in many
natural environments because most colloids and solid surfaces are negatively charged at the prevailing pH conditions[Wan and Tokunaga, 2002]. In this case, colloids may interact with the solid-water interface (SWI) at a separation distance as a result of the summation of repulsive electrostaticand attractive van der Waals interactions, which results in a
weak interaction in the secondary minimum, and an energy
barrier to attachment in the primary minimum [Franchi andO’Melia, 2003; Hahn et al., 2004; Tufenkji and Elimelech,
2005a; Kuznar and Elimelech, 2007]. Microscopic physicalor chemical heterogeneity can also lead to local reductionsin the energy barrier height and weak adhesive interactions
at a separation distance due to repulsive interactions fromneighboring regions [Elimelech et al., 2003; Hoek andAgarwal, 2006; Kozlova and Santore, 2006, 2007; Kalasin
and Santore, 2008; Duffadar and Davis, 2007, 2008; Duffadar et al., 2009]. Weak adhesive interactions at a separation distance may also occur for colloids under otherwise
1

US Salinity Laboratory, USDA, ARS, Riverside, California, USA.
Department of Environmental Sciences, University of California,
Riverside, California, USA.
2

Copyright 2011 by the American Geophysical Union.
0043-1397/11/2011WR010812

favorable attachment conditions as a result of forces that are
neglected in Derjaguin, Landau, Verwey, and Overbeek(DLVO) theory [Derjaguin and Landau, 1941; Verwey andOverbeek, 1948], such as electrosteric repulsion, hydration
forces, and Born repulsion [Elimelech et al., 1995; Israelachvili, 1992].
[3] There are many important consequences for weak adhesive interactions of colloids with the SWI at a separation
distance. First, only a fraction of the colloids that collidewith the SWI may interact because a fraction will possess
sufficient kinetic energy to overcome the weak adhesiveinteraction [Yao et al., 1971; Ryan and Elimelech, 1996;
Simoni et al., 1998; Dong et al., 2002; Shen et al., 2007].
Second, a fraction of the interacting colloids may detachfrom the SWI due to random variations in kinetic energy(diffusion) [Simoni et al., 1998; Dong et al., 2002; Shenet al., 2007]. Third, the interacting colloids may be susceptible to removal due to hydrodynamic forces during fluid
flow that cause colloids to roll on the SWI to regions where
the hydrodynamic force is weaker (grain-grain contacts or
surface roughness locations) or the adhesive force is stronger (chemical heterogeneity) [Hubbe, 1984; Bergendahl
and Grasso, 1998, 1999, 2000; Elimelech et al., 2003;
Bradford et al., 2007, 2011; Torkzaban et al., 2007, 2008].
Fourth, only a fraction of the solid surface area may contribute to colloid immobilization [Song and Elimelech,
1994; Bradford et al., 2007, 2011; Torkzaban et al., 2007,

W10503

1 of 12

W10503

BRADFORD ET AL.: MODELING COLLOID TRANSPORT AND RETENTION

2008]. Fifth, this fraction may fill up with time, and the rateof filling will depend on the colloid concentration [Adamczyk et al., 1994; Song and Elimelech, 1994; Kim et al.,2009; Bradford et al., 2009a].
[4] Mathematical models have been developed in an
attempt to describe colloid transport and retention underunfavorable attachment conditions. The most common
models employ the advective dispersion equation with firstorder colloid attachment and detachment terms [Harvey
and Garabedian, 1991; Corapcioglu and Choi, 1996;
Bolster et al., 1999; Schijven and Hassanizadeh, 2000].
The attachment coefficient is typically determined using filtration theory [Yao et al., 1971], with a modification toaccount for the fact that only a fraction of the colliding colloids attach to the SWI [Ryan and Elimelech, 1996]; i.e.,the so called sticking efficiency  (–). The value of  is frequently considered to be an empirical parameter determined by fitting to experimental data, although theory has
also been developed to predict  [Simoni et al., 1998;
Dong et al., 2002; Shen et al., 2007]. It should be noted,
however, that first-order attachment/detachment models
have been unable to account for many experimental observations. This includes the observed colloid transport andretention dependence on solution chemistry, grain size,
water velocity, and colloid concentration [Song and Elimelech, 1994; Bradford et al., 2003, 2006a, 2006b, 2007;
Tufenkji and Elimelech, 2005a, 2005b; Tong et al., 2005;
Li and Johnson, 2005; Li et al., 2005; Johnson et al.,2007; Torkzaban et al., 2007, 2008]. In particular, theshape of the colloid retention profiles (RPs) is frequently
not exponential with depth [Albinger et al., 1994; Baygentset al., 1998; Simoni et al., 1998; Bolster et al., 2000;
DeFlaun et al., 1997; Zhang et al., 2001; Redman et al.,
2001; Bradford et al., 2002, 2006a, 2006b, 2007, 2009a;
Li et al., 2004; Bradford and Bettahar, 2005, 2006; Tonget al., 2005], in contrast to predictions of first-order attachment models. In an attempt to account for these discrepancies a variety of two-site chemical nonequilibrium [Schijven
et al., 2002; Bradford et al., 2003; Tufenkji and Elimelech,
2005b], physical nonequilibrium [Cherrey et al., 2003],
chemical and physical nonequilibrium [Leij and Bradford,
2009], dual permeability [Bradford et al., 2009b; Yuan andShapiro, 2011], and stochastic models [Bradford and Toride, 2007; Shapiro and Bedrikovetsky, 2010] have beendeveloped. These models have provided an improved
description of collected experimental data, but model parameters are largely empirical and obtained by optimization
to experimental data. To date, no model exists that provide
reliable predictions of colloid transport and retention underunfavorable conditions, even under relatively simple, welldefined conditions.
[5] The objective of this work is to present a mathematical model for colloid transport and retention under unfavorable conditions that accounts for observed trends invelocity, adhesive interaction, concentration, and grain
size. Our approach considers the five factors describedabove for weak adhesive interactions. Below we describe
the mathematic model, the approach used to estimate
model parameters, and provide example simulations todemonstrate the predicted sensitivity to colloid size, grainsize, velocity, adhesive interaction, and concentration. We
then provide illustrative comparisons to experimental data.

2.

W10503

Theory[6] The dual permeability model formulation is quite
flexible and has been used to describe hyperexponential[Bradford et al., 2009b] and nonmonotonic [Yuan and Shapiro, 2011] RPs in homogeneous porous media. These simulation results suggest that the dual-permeability model
provides a promising physical interpretation for colloid
retention under unfavorable attachment conditions. Unfortunately, the dual permeability model contains a large numberof coupled parameters that have not yet been reliably measured or predicted. Consequently, the dual permeability
model has been of limited utility because its parameters
have to be fitted to experimental data and the optimization
is frequently nonunique. Below we present an improved
conceptual and mathematical description of colloid transport and retention in saturated porous media, and provide an
initial estimate of the associated dual permeability model
parameters. Our approach constrains the model parameters
and output to physically realistic values in many instances,and thereby minimizes the need for parameter optimization.
In addition, our model framework and theory helps to identify gaps in knowledge and provides a starting point for further model improvement.
2.1. Mass Balance Equations[7] Colloid transport and retention is conceptualized asshown in Figure 1. Colloids are transported through thebulk aqueous phase by advection and dispersion in region 1.
Region 2 is associated with the zone of colloid interaction
with the SWI. The thickness of region 2 is very small (onFigure 1. A conceptual picture of the flow and colloidtransport processes that are implemented in our model formulation. The subscripts 1 and 2 on parameters denote therespective regions, C is the colloid concentration in the
aqueous phase, J is the colloid flux, w is the volumetricwater content, s is the colloid concentration on the solid
phase, k2s is the first order colloid immobilization rate coefficient from region 2 to the solid phase, k12 is the mass transfer coefficient for colloids from region 1 to 2, k21 is the masstransfer coefficient for colloids from region 2 to 1, 2s
accounts for time and concentration dependent blocking,and  is the colloid sticking efficiency.

2 of 12

W10503

BRADFORD ET AL.: MODELING COLLOID TRANSPORT AND RETENTION

the order of the colloid size). Colloids may be transported byadvection and dispersion in region 2, but the velocity is much
lower than in region 1. Mass transfer of colloids to and fromregions 1 to 2 is quantified using first-order kinetic expressions. Only a fraction of the colloids in region 2 will interact
with the solid phase at any given time, and this fraction is
subject to kinetic retention and release. Immobilized colloidson the solid phase may fill up retention locations over time.
[8] The following aqueous and solid phase mass balance
equations apply to this conceptual model:@ ðw1 C1 Þ@J1 w1 k12 C1 þ w2 ð1  Þk21 C2 ;
¼
@z
@t

ð1Þ

ð2Þ

@ ðb sÞ
¼ w2 2s k2s C2 ;@t

ð3Þwhere subscripts 1 and 2 denote the respective regions, t(T; T denotes units of time) is time, z (L; L denotes unitsof length) is the depth, C (Nc L3 ; Nc denotes the numberof colloids) is the colloid concentration in the aqueous
phase, J (Nc L2 T1) is the colloid flux (the sum of the advective and dispersive fluxes), w (–) is the volumetricwater content, s (Nc M1 ; M denotes units of mass) is thecolloid concentration on the solid phase, b (M L3) is thebulk density, k2s (T1) is the first order colloid immobilization rate coefficient from region 2 to the solid phase, k12(T1) is the mass transfer coefficient for colloids fromregion 1 to 2, and k21 (T1) is the mass transfer coefficientfor colloids from region 2 to 1. The parameter 2s (–)
accounts for time and concentration dependent blocking
using a Langmuirian approach as [Adamczyk, 1994]:2s ¼ 1 

s
;
smax

ð4Þ

where smax (Nc M1) is the maximum solid phase concentration of retained colloids. The total water flux (qt ; L T1),
volumetric water content (wt ), and flux concentration of
colloids (Ct, Nc L3) are given in the model as [Simůnek
and van Genuchten, 2008]:
qt ¼ q1 þ q2 ;ð5Þ

wt ¼ w1 þ w2 ;ð6Þ

q1 C1 þ q2 C2;
q1 þ q2

ð7Þ

Ct ¼

algorithm [Marquardt, 1963] to determine model parameters from measured breakthrough curves (BTCs) and/or
RPs. For the simulations discussed below, a third-typeboundary condition was used at the inlet, and a concentration gradient of zero was fixed at z equal to the outlet depth.
2.2. Determination of Parameters[10] Our conceptual model (Figure 1) assumes that
region 2 is associated with the zone of colloid interaction
with the SWI. Consequently, an initial estimate for thevalue of w2 was obtained as the product of the geometricsolid surface per unit volume (As, L1) and a boundary
thickness (L2, L) that is given asw2 ¼ As L2 ¼ As ð2rc þ hÞ;@ ðw2 C2 Þ@J2þ w1 k12 C1  w2 ð1  Þk21 C2¼
@z@t w2 2s k2s C2 ;where q1 and q2 (L T1) are the Darcy velocities for regions
1 and 2, respectively.
[9] The model outlined above has been implemented into
the COMSOL software package (COMSOL, Inc., Palo
Alto, CA 94,301) and into HYDRUS1D [Simůnek
et al.,
2008]. The model is coupled with a nonlinear least squares
optimization routine based upon the Levenberg-Marquardt

W10503

ð8Þ

where rc (L) is the colloid radius and h (L) is the interactionenergy separation distance. The value of h was initially set
equal to 20 nm in this work to be consistent with a secondary minimum interaction. The value of w1 follows directlyfrom equation (6) and the measured value of wt .
[11] The value of q2 can be determined from w2 and the
median velocity in region 2 (v2, L T1) which occurs at a
distance of 0.5L2 from the SWI as
q2 ¼ w2 v2 :

ð9ÞBradford et al. [2011] presented pore scale simulations ofwater flow in sphere packs, and determined the cumulative
density function (CDF) of water velocity adjacent to the
SWI. This information was further extended by these authors
using scaling and interpolation techniques to predict the
CDF of v2 for a range of grain sizes and distributions, watervelocities, and colloid sizes. The value of v2 can be estimated from this information as
v2 ¼

L2
LV

     
d50
qt
"
v ;
d50
qt
" 50

ð10Þ

where the superscript  denotes parameters in the reference
simulation, LV (L) is the reference simulation voxel length
equal to 1.077  106 m, d50 (L) is the median grain diameter, " ðÞ is the porosity, and v50 (L T1) is the medianboundary water velocity in the reference simulation equal to
7.12  106 m s1. The values d50
and qt =" were 100 lm
5
1
and 2.8  10 m s , respectively. The value of q1 follows
directly from equation (5) and the selected value of qt. Values of the hydrodynamic dispersion coefficients D1 and D2(L2 T1), were set equal to the product of their respective
pore water velocity and a constant dispersivity (e.g., 0.1
times the length of the simulation domain).
[12] The rate of colloid mass transfer to the SWI undersaturated conditions depends on diffusion, sedimentation,and interception. Colloid filtration theory [Yao et al., 1971] iscommonly used to estimate this mass transfer coefficient as
k12 ¼

3ð1  w Þ
vavg :
2d50

ð11Þ

Here vavg (L T1) is the average pore water velocity and ðÞ is the collector efficiency. Correlation equations topredict  as a function of system variables have been

3 of 12

W10503

BRADFORD ET AL.: MODELING COLLOID TRANSPORT AND RETENTION

developed from pore scale simulations studying colloid
mass transfer to the collector surface in the sphere-in-cellmodel [Rajagopalan and Tien, 1976; Tufenkji and Elimelech, 2004]. In this work we employ the correlation of
Tufenkji and Elimelech [2004].
[13] Values of k12 in equations (1) and (2) are multipliedby  because only a fraction of the colloids will interact
with the solid surface under unfavorable attachment conditions. Several theoretical approaches have been developed
to estimate  [Simoni et al., 1998; Dong et al., 2002; Shenet al., 2007]. The kinetic energy method assumes that the
distribution of kinetic energies of diffusing colloids can be
described by Maxwell’s function [Chandrasekhar, 1943],
and that DLVO or extended DLVO energy profiles accurately predict the strength of the interaction between colloids and collectors. The value of  is then related to thefraction of diffusing colloid particles that possess kinetic
energy less than a given dimensionless (divided by kbTk,
where kb¼ 1.38  1023 J K1 is the Boltzmann constantand Tk is the temperature in degree Kevin, K) depth of theinteraction energy minimum (min ) asZmin pﬃﬃﬃﬃ
2 
pﬃﬃﬃ expðÞd;
¼


duf ¼



d50 1  Sf
:
2Nf

ð14Þ

[16] Here Sf (–) is the fraction of the collector surface
where the torque balance is favorable for immobilization,
and Nf (–) is the number of favorable locations of equal
size and distribution on the collector surface. On a smooth,chemically homogeneous collector the value of Nf may
reflect the average number of grain-grain contacts. Unlessotherwise noted we assume in this work a value of Nf ¼ 5
to be consistent with this condition. Alternatively, if surface roughness or chemical heterogeneity is controlling thevalue of Sf then Nf is given as
d50 Sf
:
2df

ð15Þ

ð12Þ

where  (–) is a dummy variable of integration. The kinetic
energy model implies that a complementary fraction of colloids in region 2, 1  , may diffuse away from the solid
surface at any given time and values of k21 are therefore
multiplied by (1  ) in equations (1) and (2).
[14] It should be noted that filtration theory assumes thatinteracting colloids are immobilized instantaneously. In
this work, we assume that interacting colloids are mobile
and may roll along the solid surface until they encounter alocation where the torque balance is favorable for immobilization. The rate of colloid immobilization is given by k2s.
Recent experimental information under unfavorable attachment conditions supports the assumption of colloid rollingin region 2 [Kuznar and Elimelech, 2007], and colloid
immobilization near grain-grain contacts [Bradford et al.,2005, 2006a; Xu et al., 2006; Li et al., 2006; Yoon et al.,2006; Gaillard et al., 2007; Tong et al., 2008] and locations associated with surface roughness [Choi et al., 2007]
and/or chemical heterogeneity [Hoek and Agarwal, 2006;
Kozlova and Santore, 2006, 2007; Kalasin and Santore,
2008; Duffadar and Davis, 2007, 2008; Duffadar et al.,2009]. However, we are unaware of any experimental ortheoretical information to directly calculate k2s for the complex geometries found in porous media. In the absence of
direct measurements, we initially assume that k2s is inversely related to the time it takes colloids in region 2 to travel
to a retention location on the solid surface as
v2
:
duf

not instantaneously immobilized on the solid surface when
duf > 0 but continue to be transported in region 2 until
arriving at a retention location.
[15] The value of duf may be considered as an empirical
optimization parameter or a physical estimate may be
given as

Nf ¼

0

k2s ¼

W10503

ð13Þ

Here duf (L) is the average distance that colloids must travel
in region 2 to find a retention location. As duf approaches 0
the value of k2s becomes large so that colloids in region 2
are almost instantaneously immobilized on the solid phaseas assumed in filtration theory. Conversely, colloids are

Here df (L) is equal to the average size of the heterogeneity
parallel to the collector surface. If surface molecules arehindering the immobilization of colloids with the solid surface, then it is logical to anticipate that other functional
forms for Nf would need to be developed.
[17] Bradford et al. [2011] presented a detailed approachto predict the value of Sf based on the balance of appliedhydrodynamic (Tapplied, ML2 T2) and resisting adhesive(Tadhesion, ML2 T2) torques. In summary, the value of Sf on
a smooth, chemically homogeneous collector (Sf ) is given as
1 1lnðTadhesion Þ  
pﬃﬃﬃ
Sf ¼ þ erf
;
2 2
 2

ð16Þ

where  and  are the mean and variance of the lognormal
CDF of Tapplied determined from results of pore scale simulations and scaling approaches. The value of  for the lognormal distribution is defined as ¼ lnðT50 Þ:

ð17Þ

The value of T50 (ML2 T2) corresponds to the medianvalue of Tapplied on the SWI. Bradford et al. [2011] provide
a detailed explanation of how to predict this value for various colloids, porous media, and water flow conditions. In
equation (16) the CDF of Tapplied is evaluated at Tadhesion todetermine Sf . The value of Tadhesion may be determinedusing DLVO and JKR theories [Derjaguin and Landau,
1941; Verwey and Overbeek, 1948; Johnson et al., 1971] orfrom a coefficient of rolling friction [Duffadar and Davis,
2008]. Details are again given by Bradford et al. [2011]. It
should be mentioned that Sf does not account for the potential influence of surface roughness and/or nanoscale chemical heterogeneity which are expected to produce a minimum
value of Sf (Sfmin ). We therefore determine Sf as

4 of 12

Sf ¼ Sfmin þ Sf :ð18Þ

W10503

BRADFORD ET AL.: MODELING COLLOID TRANSPORT AND RETENTION

Unless otherwise noted we set the value of Sfmin in thiswork equal to 0.03 based on experimental data presented
by Bradford et al. [2009a].
[18] It should be mentioned that the value of smax inequation (4) is related to Sf as [Kim et al., 2009; Bradford
et al., 2009a]
smax ¼

ð1  ÞAs Sf
;
Ac b

ð19Þ

where Ac (L2 N1) is the cross section area per colloid, and(–) is the porosity of a monolayer packing of colloids on
the solid surface. In this work we assume a value of¼ 0:5 in all simulations based on information presented
by Johnson and Elimelech [1995].
[19] Equations (1)–(3) are mathematically equivalent tothe conventional first-order attachment/detachment model
when k2s ¼ 0 and q2 ¼ 0. In this case, the attachment anddetachment rate coefficients are equal to k12 and ð1  Þ
k21 , respectively. The value of ð1  Þk21 has important
implications for the shape of the BTCs in this model. When
ð1  Þk21 is low then concentration tailing occurs after recovery of the unretarded breakthrough curve. However, a
low value of detachment cannot account for many experimental observations of colloid retention under unfavorable
attachment conditions, such as hyperexponential RPs[Albinger et al., 1994; Baygents et al., 1998; Simoni et al.,
1998; Bolster et al., 2000; DeFlaun et al., 1997; Zhang
et al., 2001; Redman et al., 2001; Bradford et al., 2002; Liet al., 2004; Bradford and Bettahar, 2005]. Conversely,
when ð1  Þk21 is high then the first-order attachment/
detachment model approaches linear equilibrium conditions with a retardation coefficient (R) equal to 1 þ k12 =ð1  Þk21 [Schijven and Hassanizadeh, 2000]. As a first
approximation we assume in this work that k21 ¼ k12 andthis yields a value of R ¼ 1 þ =ð1  Þ that is controlled
by chemistry as expected. This assumption is certainly justifiable under highly unfavorable attachment conditions
with small  and R close to 1 because k12  ð1  Þk21and is consistent with results reported in the literature[Gargiulo et al., 2007, 2008].
[20] It should be mentioned that the outlined model formulation approaches filtration theory predictions underfavorable attachment conditions. In particular, the values of, Sf, and 2s go to 1, ð1  Þk21 goes to 0, k12 is determinedusing filtration theory [equation (11)], and k2s becomes verylarge such that colloids that enter region 2 are instantaneously retained.

3.

Results and Discussion

3.1. Predictions[21] Under unfavorable attachment conditions the model
predicts a complex coupling between many physical andchemical factors. Below we present representative simulations to demonstrate the sensitivity of the model output tospecific input parameters. The simulation domain length
was 0.1 m to be representative of typical packed column
studies, the initial concentration in the simulation domain
was zero, and the input pulse duration was 3 pore volumes(PV). BTCs are presented herein by plotting the relative
effluent concentrations [Ct/Ci ; where Ci (Nc L3) is theW10503

input concentration] as a function of PV. RPs are shown on
a semilog plot of normalized total concentration in region 2
and on the solid phase [ðb s þ wt C2 Þ=Ci ] as a function of
distance from the inlet.
[22] We first examine the influence of grain size andmin on colloid transport and retention. Figure 2 presents
simulated BTCs (Figures 2a and 2c) and RPs (Figures 2band 2d) when rc ¼ 500 nm, qt ¼ 0.1 cm min1, wt ¼ 0:35,
Ci ¼ 106 Nc mL1, and d50 ¼ 50, 100, 250, 500, and 1000
lm. The value of min ¼ 0:15 in Figures 2a and 2b, andmin ¼ 0:5 in Figures 2c and 2d. A systematic trend ofincreasing colloid retention occurs with decreasing d50 for
both min equal to 0.15 and 0.5. This occurs in part because
filtration theory predicts that decreasing d50 will increase
k12 in equation (11). In contrast to the predicted filtration
theory, however, the RPs are not always log linear with distance. Smaller values of d50 and min tend to produce more
hyperexponential RPs (a greater amount of retention near
the inlet). Hyperexponential RPs have been experimentallyobserved by many researchers under highly unfavorable
attachment conditions [Albinger et al., 1994; Baygentset al., 1998; Simoni et al., 1998; Bolster et al., 2000;
DeFlaun et al., 1997; Zhang et al., 2001; Redman et al.,
2001; Bradford et al., 2002; Li et al., 2004; Bradford andBettahar, 2005]. In addition, RPs have also been observed
to become more hyperexponential with decreasing d50[Bradford et al., 2002, 2003]. The implemented mathematical model provides a clear explanation for these observations as will be discussed below.
[23] Hyperexponential RPs occur because w2 [equation(8)], v2 [equation (10)], q2 [equation (9)], and k2s [see equations (13) and (14)] all increase with decreasing d50. Masstransfer to region 2 is controlled by exchange with region
1 (k12 ) and J2 at the inlet (J2 ). The value of J2 increases
with decreasing d50 and this leads to a greater amount of
colloid retention near the inlet. As J2 is depleted by colloid
retention with distance, then the mass transfer to region 2 iscontrolled by exchange with region 1 (k12 ). Consequently,
the relative importance of J2 on colloid retention decreases
with increasing min [due to its influence on  given byequation (12)] and with distance from the inlet. Consistent
with this explanation, Li et al. [2004], Tufenkji and Elimelech [2005a], and Bradford et al. [2007] observed that colloid RPs became less hyperexponential with increasing
min . The above information demonstrates that accounting
for the median hydrodynamics [equations (8)–(10)] of colloids near the solid surface provides a viable explanation forthe dependence of hyperexponential RPs on d50 and min .
[24] It should be emphasized that hyperexponential RPs
are not an experimental artifact of the boundary conditions
at the soil surface. This behavior has also been observed tooccur at textural interfaces when water flows from a coarser
to a finer textured material [Bradford et al., 2005], especially for larger colloids. As discussed above, the model
predicts that w2 and q2 increases with decreasing d50. Consequently, J2 at the textural interfaces changes and this difference will increase with the contrast in d50. However, it ispresently unclear how concentrations in regions 1 and 2 mix
at textural interfaces (e.g., mixing may be complete, partial,
or absent). Additional research is warranted on this topic.
[25] Figure 3 presents simulated BTCs and RPs when d50 ¼
250 lm, qt ¼ 0.1 cm min1, wt ¼ 0:35, Ci ¼ 106 Nc mL1,

5 of 12

W10503

BRADFORD ET AL.: MODELING COLLOID TRANSPORT AND RETENTION

W10503Figure 2. Simulated (a and c) BTCs and (b and d) RPs when rc ¼ 500 nm, qt ¼ 0.1 cm min1,
wt ¼ 0:35, Ci ¼ 106 Nc mL1, min ¼ ða and bÞ 0.15 and (c and d) 0.5, and d50 ¼ 50, 100, 250, 500,
and 1000 lm. The BTCs are presented herein by plotting the relative effluent concentrations (Ct/Ci) as afunction of pore volumes. The RPs are shown on a semilog plot of normalized total concentration inregion 2 and on the solid phase [ðb s þ wt C2 Þ=Ci ] as a function of distance from the inlet.
and rc ¼ 250, 500, 1000, and 2000 nm. Since min (the secondary minimum) is also a function of rc, values of minwere determined from DLVO calculations assuming a zeta
potential of 30 mV for both the colloid and collector, anda ionic strength (IS) of 2 mM of a monovalent electrolyte
solution. Under these assumptions the value of min was
0.039, 0.079, 0.158, and 0.316 when rc equals 250, 500,
1000, and 2000 nm, respectively. A systematic increase incolloid retention occurs with increasing rc. In addition, allof the RPs are hyperexponential and the amount of retention near the inlet increases with rc (profiles become more
hyperexponential). These results are consistent with experimental observations reported in the literature [Bradford
et al., 2002, 2003, 2007; Tufenkji and Elimelech, 2005a]. It
should be mentioned that if Figure 3b was plotted using anonlog scale then the hyperexponential profiles for thesmallest colloids would not have been observable due totheir low values. Similarly, the sensitivity of the analytic
approach (spectrophotometer or fluorometer) to determine
RPs may not always be sufficient to quantify such subtle
differences in shape.
[26] A detailed explanation for the observations in Figure 3 can be obtained from the implemented mathematicalmodel. Differences in the amount and shape of the RPs are
due to the influence of rc on J2 and on k12 . Values of w2[equation (8)], v2 [equation (10)], q2 [equation (9)], and k2s[equation (13)] all increase with increasing rc. This produces a nonlinear increase in J2 with increasing rc and consequently more hyperexponential RPs. The amount ofretention away from the inlet also increases with increasing
rc because of a corresponding increase in k2s and k12
[equations (11) and (12)]. Slight differences in the apparent
pulse duration are due to increases in w2 and corresponding decreases in w1 with increasing rc.
[27] Additional simulations were run using the same
conditions as in Figure 3, but at a IS ¼ 20 mM (data notshown). Observed trends with rc were similar to that shownin Figure 3. Similar to differences in Figure 2 with min ,increasing the solution IS produced greater amounts ofretention and profiles that are less hyperexponential
because the contribution of k12 to retention increases relative to that of J2 . It should be mentioned that if a constantvalue of min and  were assumed with rc in Figure 3 thendifferences in the retention away from the inlet would havebeen controlled by k12 (instead of k12 ). Filtration theory
predicts a nonlinear dependence of k12 on rc, with a minimum around rc ¼ 1000 nm. Consequently, a constant valueof min and  produces a more complex dependence of colloid retention on rc than that shown Figure 3.
[28] Figure 4 presents simulated BTCs and RPs whenrc ¼ 500 nm, d50 ¼ 250 lm, qt ¼ 0.1 cm min1,
wt ¼ 0:35, min ¼ 0:5, and Ci ¼ 106, 107, 108, 109, and

6 of 12

W10503

BRADFORD ET AL.: MODELING COLLOID TRANSPORT AND RETENTION

W10503Figure 3. Simulated (a) BTCs and (b) RPs when d50 ¼ 250 lm, qt ¼ 0.1 cm min1, wt ¼ 0:35, Ci ¼
106 Nc mL1, and rc ¼ 250, 500, 1000, and 2000 nm. Values of min were determined from DLVO calculations when the zeta potential of the colloid and collector were 30 mV and the IS of a monovalent
electrolyte solution was 2 mM.
1010 Nc mL1. This data was normalized by Ci so that differences in transport and retention behavior were apparent.
Results indicate that above a threshold in Ci (or input pulse
duration) the normalized BTCs become sensitive to Ci. In
this case, the threshold value is greater than Ci ¼ 108 Nc
mL1. Values of Ci above this threshold influence theshape of the BTCs and the RPs, with the BTCs exhibiting
blocking behavior (a decreasing rate of retention over time)
and the RPs transitioning from exponential to uniform inshape with depth. Our mathematical formulation indicatesthat this occurs as a result of decreases in 2s [equations (4)and (19)] with filling of retention sites. Furthermore, themodel predicts a sensitivity of colloid transport and retention to Ci that is a function of colloid size, grain size, chemistry, and velocity. For example, similar simulations tothose shown in Figure 4 were conducted with min ¼ 0:15
(data not shown). In this case, analogous concentration
effects on the BTCs and RPs occurred to those shown in
Figure 4 but over a smaller range in retention values because
the value of min determines  and therefore the maximum
amount of retention. However, when min ¼ 0:15 the RPstransitioned from hyperexponential to uniform in shape withdepth as Ci varied from 106 to 1010 Nc mL1. Simulatedtransport and retention trends with Ci are qualitatively

consistent with experimental observations reported in the literature [Tan et al., 1994; Lindqvist et al., 1994; Liu et al.,
1995; Bradford and Bettahar, 2006; Bradford et al.,2009a]. However, additional research is warranted to test
the predicted concentration dependency over a wider rangeof experimental conditions and to improve the conceptualmodeling framework. For example, we postulate that thevalue of  may decrease with higher values of Ci because
the number of colloid collisions increases with Ci and the kinetic energy distribution may therefore not be accurately
described by the Maxwellian distribution [Chapman and
Cowling, 1991; Bradford and Bettahar, 2006].
[29] Figure 5 presents simulated BTCs and RPs whenrc ¼ 500 nm, d50 ¼ 250 lm, wt ¼ 0:35, min ¼ 0:15, Ci ¼
106 Nc mL1, and qt equals 0.01, 0.05, 0.1, 0.2, and 0.4 cm
min1. A systematic decrease in colloid retention occurs
with increasing qt. This result is consistent with experimental observations reported in the literature [Kretzschmar andSticher, 1998; Schijven and Hassanizadeh, 2000]. Filtration theory predicts that k12 is approximately proportional
to vavg raised to the 1/3 power [Schijven and Hassanizadeh,2000]. However, the overall rate of advection is proportional to vavg. Consequently, filtration theory predicts thatcolloid retention will decrease with increasing vavg.

Figure 4. Simulated (a) BTCs and (b) RPs when rc ¼ 500 nm, d50 ¼ 250 lm, qt ¼ 0.1 cm min1,
wt ¼ 0:35, min ¼ 0:5, and Ci ¼ 106, 107, 108, 109, and 1010 Nc mL1.
7 of 12

W10503

BRADFORD ET AL.: MODELING COLLOID TRANSPORT AND RETENTION

W10503

1
(a)

CT\Ci

0.8
qT (cm min-1)

0.6

0.01
0.05
0.1
0.2
0.4

0.4
0.2
0
0

2

4
6
Pore Volumes

8

10

Figure 5. Simulated (a) BTCs and (b) RPs when rc ¼ 500 nm, d50 ¼ 250 lm, wt ¼ 0:35, min ¼ 0:15,
Ci ¼ 106 Nc mL1, and qt equals 0.01, 0.05, 0.1, 0.2, and 0.4 cm min1.
[30] The retention profile in Figure 5b is almost exponential at the lowest value of qt ¼ 0.01 cm min1, but
becomes increasingly hyperexponential for increasing qt.
Values of v2, q2, and (in these simulations) k2s increase
with increasing qt. This produces an increase in J2 withincreasing qt and consequently more hyperexponential RPs.
This predicted trend has not yet been experimentally verified. In fact, Li et al. [2004] and Bradford et al. [2007]
reported that RPs became less hyperexponential withincreasing velocity. This discrepancy could arise from a velocity dependency on  [Johnson et al., 2007] that is notaccounted for in equation (12), analytical insensitivity tolower concentrations of retained colloids, or the dependency of RPs on Ci as previously demonstrated. In addition,
the simulated value of Sf was not influenced much bychanges in qt for the considered min and rc. However, variations in min and rc over a wider range of conditions is
expected to influence Sf [Torkzaban et al., 2007; Bradford
et al., 2011] and a decreasing value of Sf with qt would
increase the sensitivity of RPs to Ci. Additional research iswarranted on all of these velocity related topics. It shouldbe mentioned that similar simulations to those shown in
Figure 5 were conducted using min ¼ 0:5 (data notshown). Similar to Figure 2, increasing min to 0.5 produced a greater amount of retention and profiles that wereless hyperexponential than those shown in Figure 5. However, overall trends of decreasing retention with increasing
qt were preserved.
3.2. Analysis of Experimental Data
[31] The outlined model was used to describe transportand retention data for monodispersed suspensions of fluorescent 1 lm carboxyl latex microspheres [Bradford andBettahar, 2006; Bradford et al., 2007, 2009a] and E. coli
O157:H7/pGFP [Bradford et al., 2006b] in saturated, packcolumn experiments under unfavorable attachment conditions. The selected data for analysis exhibits a wide rangeof transport and retention behavior that is intended to illustrate the capabilities and limitations of this model.
[32] Experimental details are given in the indicated publications, but will be briefly highlighted below. Selected
sieve sizes of sand were wet packed into 4.8 cm diameterby 15 cm long columns that were equipped with an adjustable length flow adapter. A salt cleaning method [Bradford

et al., 2002] was used to remove trace amounts of clay
from the packed sand. An electrolyte solution with aselected chemistry was prepared for the resident, tracer,and eluting solutions. The sand in the columns was equilibrated with the resident solution, and then a tracer solutionof microspheres or E. coli O157:H7 was pumped throughthe packed columns at a steady flow rate for several PV followed by continued flushing with the eluting solution. The
sand in the columns was incrementally excavated into vials
containing excess eluting solution, the vials were shaken,
an aliquot was collected for analysis, and the mass of sandin the vials was determined after drying. The concentrations of fluorescent microspheres or E. coli O157:H7 were
determined in the solution from the column effluent andvials using a fluorometer. A mass balance was conductedbased on measured concentrations in the influent, effluent,and sand.
[33] Figure 6 presents observed and simulated BTCs andRPs for carboxyl latex microspheres under several different
conditions of solution chemistry (curves 1, 2, and 3 were for
pH ¼ 7 and IS ¼ 1 mM, pH ¼ 10 and IS ¼ 31 mM, andpH ¼ 10 and IS ¼ 81 mM, respectively), colloid hydrophobicity (curves 1, 2, and 3 were hydrophobic, hydrophilic,and hydrophilic, respectively), and Ci (curves 1, 2, and3 were for 3.9  107, 3.9  109, and 2.7  107 Nc mL1,
respectively), but for similar values of rc ¼ 500 nm, d50 ¼
150 lm, qt ¼ 0.1 cm min1, wt ¼ 0:35, and pulse duration(1.5–2.1 PV). Curves 1, 2, and 3 approached hyperexponential, uniform, and exponential distributions with depth. The
simulations were obtained by optimizing values of  and L2
to the curves, and the dispersivities were set equal to 0.01
cm be consistent with the data. The coefficient of linear
regression (r2) between observed and simulated data ranged
from 0.85 to 0.94, and this indicates that the optimized modelgave a reasonable description to a wide range of experimental data and retention profile shapes. It should be mentionedthat the time dependent retention behavior observed incurves 1 and 2 was incompletely characterized by the model
because values of smax were predicted instead of optimized.
[34] Optimized values of  to curves 1, 2, and 3 wereequal to 0.025, 0.045, and 0.055, respectively. These valuesof  account for differences in solution chemistry and colloid hydrophobicity, and increased with IS due to compression of the double layer thickness and a corresponding

8 of 12

W10503

BRADFORD ET AL.: MODELING COLLOID TRANSPORT AND RETENTION

W10503Figure 6. Observed and simulated (a) BTCs and (b) RPs for carboxyl latex microspheres under several
different conditions of solution chemistry (curves 1, 2, and 3 were for pH ¼ 7 and IS ¼ 1 mM, pH ¼ 10
and IS ¼ 31 mM, and pH ¼ 10 and IS ¼ 81 mM, respectively), colloid hydrophobicity (curves 1, 2, and3 were hydrophobic, hydrophilic, and hydrophilic, respectively), and Ci (curves 1, 2, and 3 were for3.9  107, 3.9  109, and 2.7  107 Nc mL1, respectively), but for similar values of rc ¼ 500 nm,
d50 ¼ 150 lm, qt ¼ 0.1 cm min1, wt ¼ 0:35, and pulse durations (1.5–2.1 PV). Values of  were optimized to each data set [Bradford and Bettahar, 2006; Bradford et al., 2007, 2009a], as well as a single
value of L2.
increase in the depth of the secondary minimum. To
improve the description of the RPs it was also necessary
to optimize a single value of L2 ¼ 4 lm. The value of L2
influences J2 by changing w2 , v2, and q2 that were basedon information from pore scale water flow simulations onsmooth, spherical collector surfaces. Differences in predicted values of L2 are therefore likely with grain roughness, angularity, and size distribution.
[35] In addition, our modeling approach is based on anumber of assumptions with regard to the chemical interaction that may be violated in some instances. For example,
our model predictions assumed that the colloid was muchlarger than any physical or chemical heterogeneity on thecolloid and/or collector surface, such that colloids on thecollector surface experienced an effective secondary minimum interaction. It should be noted that DLVO theory will
not be adequate to characterize the interaction energy when
the size of the colloid approaches that of the heterogeneity
on the collector surface or when surface macromolecules
produce significant non-DLVO forces, and this will lead tocorresponding errors in predictions for , Sf, smax, k2s, andk21. To overcome these limitations additional model
assumptions, theory, and/or parameter optimization are
required. An illustrative example is given below.
[36] Some bacteria species have been reported in the literature to exhibit nonmonotonic RPs with depth [Tong et al.,2005; Bradford et al., 2006b]. Observed and simulated
BTCs and nonmonotonic RPs for E. coli O157:H7 are shownin Figure 7 in several different pore volumes (times). Simulations employed values of rc ¼ 250 nm, d50 ¼ 150 lm, qt ¼
0.1 cm min1, wt ¼ 0:34, Ci ¼ 1.72  108 Nc mL1, and apulse duration of 75 min that were based on measurements.
The values of  ¼ 0:05, Sf ¼ 0.007, D2 ¼ 0.01  D1, k21 ¼
0, and Nf ¼ 0.01 had to be optimized by trial and error
because of the coupling between parameters and nonuniqueness in the fit. The simulation for E. coli O157:H7 clearly

demonstrates that our modeling framework is capable of producing nonmonotonic RPs (r2 ¼ 0.76). Furthermore, predictions for E. coli O157:H7 based on these fitted parametervalues were able to accurately capture the general transportand retention behavior in several different sands at different
velocities [Bradford et al., 2006b] (data not shown), and aslowly moving and spreading nonmonotonic retention profile
with time [Tong et al., 2005].
[37] In these simulations we assumed a value of rc ¼ 250
nm based on the width of the rod shaped cells, because
they will likely orient horizontally on the SWI. The valueof D2 ¼ 0.01  D1 is consistent with the much lower velocities of rod shaped cells next to the SWI than in region 1.
Values of k21 and Nf were much lower than those predictedpreviously, which suggests hindered cell release and immobilization by surface macromolecules. Additional researchis needed to predict all of these fitted parameters over awide range of conditions.

4.

Summary and Conclusions[38] A mathematical modeling formulation was presented to describe colloid transport in saturated porous
media under unfavorable attachment conditions. The modelaccounts for colloid transport in the bulk aqueous phaseand adjacent to the solid phase, and rates of colloid collision, interaction, release, and immobilization on the solid
phase. Colloid immobilization depends on all these factors.
However, the relative importance of a given process canchange with transport distance and filling of retention locations. The exact location and mechanism of retention is not
explicitly considered in this work, just the rate of immobilization which depends on the velocity adjacent to the solid
phase and the amount and distribution of retention locations (e.g., grain-grain contacts, surface roughness, andchemical heterogeneity).

9 of 12

W10503

BRADFORD ET AL.: MODELING COLLOID TRANSPORT AND RETENTION

W10503Figure 7. Observed and simulated (a) BTCs and (b) RPs for E. coli O157:H7. The value of rc ¼ 250 nm,
d50 ¼ 150 lm, qt ¼ 0.1 cm min1 and wt ¼ 0:35 were based on measurements. The values of  ¼ 0:05,
Sf ¼ 0.007, D2 ¼ 0.01 D1, k21 ¼ 0, and k2s ¼ 0.01 min1 were optimized to the E. coli O157:H7 data
[Bradford et al., 2006b].
[39] Simulations presented herein indicate that our modeling formulation is at least qualitatively consistent withobserved trends for retention with d50, rc, Ci, and qt for
many systems. Furthermore, the model provides a clear
conceptual explanation for the causes of hyperexponential,
exponential, uniform, and nonmonotonic RPs without
invoking hypotheses with regard to colloid heterogeneity,
aggregation, or multiple deposition rates. Hyperexponential
RPs occur when J2 is larger than the net colloid exchange
rate from the bulk aqueous phase to the solid surface. Thiseffect increases for large colloid sizes, smaller grain sizes,
and higher velocities under highly unfavorable attachment
conditions. In contrast, exponential RPs occur when the net
colloid exchange rate from the bulk aqueous phase to thesolid surface dominates such as for smaller colloid sizes,
larger grain sizes, lower velocities, and conditions that are
more favorable for attachment. Uniform RPs occur whensmall amounts of retention sites are filled at higher input
concentration and/or longer pulse durations. Nonmonotonic
profiles reflect an increase in the residence time of mobile
colloids on the solid phase due to a decrease in the rates of
release and immobilization, presumably due to surface
macromolecules.
[40] The research presented herein also helps to identify
areas where additional studies are needed. For example,
what are the implications of J2 at textural interfaces and inheterogeneous soils under conditions that produce hyperexponential RPs? Is the value of  also a function of velocity
and concentration ? How do we account for transients in solution chemistry in the above model formulation? How do
we predict model parameters when DLVO theory is not
valid and on rough surfaces ?
[41] Acknowledgments. This research was supported by the Agricultural and Industrial Byproducts project (NP 214) of the USDA-ARS. Mention of trade names and company names in this manuscript does not imply
any endorsement or preferential treatment by the USDA.

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.

Albinger, O., B. K. Biesemeyer, R. G. Arnold, and B. E. Logan (1994),
Effect of bacterial heterogeneity on adhesion to uniform collectors by
monoclonal populations, FEMS Microbiol. Lett., 124(3), 321–326.
Baygents, J. C., J. R. Glynn Jr., O. Albinger, B. K. Biesemeyer, K. L.
Ogden, and R. G. Arnold (1998), Variation of surface charge density inmonoclonal bacterial populations: Implications for transport through porous media, Environ. Sci. Technol., 32, 1596–1603.
Bergendahl, J., and D. Grasso (1998), Colloid generation during batchleaching tests: Mechanics of disaggregation, Colloid Surf. A: Physicochem. Eng. Aspects, 135, 193–205.
Bergendahl, J., and D. Grasso (1999), Prediction of colloid detachment in amodel porous media: Thermodynamics, AIChE J., 45, 475–484.
Bergendahl, J., and D. Grasso (2000), Prediction of colloid detachment in amodel porous media: Hydrodynamics, Chem. Eng. Sci., 55, 1523–1532.
Bolster, C. H., A. L. Mills, G. M. Hornberger, and J. S. Herman (1999),
Spatial distribution of deposited bacteria following miscible displacement experiments in intact cores, Water Resour. Res., 35, 1797–1807,
doi:10.1029/1999WR900031.
Bolster, C. H., A. L. Mills, G. Hornberger, and J. Herman (2000), Effect of
intra-population variability on the long-distance transport of bacteria,
Ground Water, 38(3), 370–375.
Bradford, S. A., and M. Bettahar (2005), Straining, attachment, and detachment of Cryptosporidium oocysts in saturated porous media, J. Environ.
Qual., 34, 469–478.
Bradford, S. A., and M. Bettahar (2006), Concentration dependent colloidtransport in saturated porous media, J. Contamin. Hydrol., 82, 99–117.
Bradford, S. A., and N. Toride (2007), A stochastic model for colloid transport and deposition, J. Environ. Qual., 36, 1346–1356.
Bradford, S. A., S. R. Yates, M. Bettahar, and J. Simůnek(2002), Physicalfactors affecting the transport and fate of colloids in saturated porous
media, Water Resour. Res., 38, 1327(12), doi:10.1029/2002WR001340.
Bradford, S. A., J. Simůnek,
M. Bettahar, M. Th. van Genuchten, and S. R.
Yates (2003), Modeling colloid attachment, straining, and exclusion insaturated porous media, Environ. Sci. Technol., 37, 2242–2250.
Bradford, S. A., J. Simůnek,
M. Bettahar, Y. F. Tadassa, M. Th. van Genuchten, and S. R. Yates (2005), Straining of colloids at textural interfaces, Water Resour. Res., 41, W10404, doi:10.1029/2004WR003675.
Bradford, S. A., J. Simůnek,
M. Bettahar, M. Th. van Genuchten, and S. R.
Yates (2006a), Significance of straining in colloid deposition: Evidence
and implications, Water Resour. Res., 42, W12S15, doi:10.1029/2005
WR004791.
Bradford, S. A., J. Simůnek,
and S. L. Walker (2006b), Transport andstraining of E. coli O157:H7 in saturated porous media, Water Resour.
Res., 42, W12S12, doi:10.1029/2005WR004805.
Bradford, S. A., S. Torkzaban, and S. L. Walker (2007), Coupling of physical and chemical mechanisms of colloid straining in saturated porous
media, Water Res., 41, 3012–3024.
Bradford, S. A., H. N. Kim, B. Z. Haznedaroglu, S. Torkzaban, and S. L.
Walker (2009a), Coupled factors influencing concentration dependent
colloid transport and retention in saturated porous media, Environ. Sci.
Technol., 43, 6996–7002.

10 of 12

W10503

BRADFORD ET AL.: MODELING COLLOID TRANSPORT AND RETENTION
Bradford, S. A., S. Torkzaban, F. Leij, J. Simůnek,
and M. Th. van Genuchten (2009b), Modeling the coupled effects of pore space geometry
and velocity on colloid transport and retention, Water Resour. Res., 45,
W02414, doi:10.1029/2008WR007096.
Bradford, S. A., S. Torkzaban, and A. Wiegmann (2011), Pore-scale simulations to determine the applied hydrodynamic torque and colloid immobilization, Vadose Zone J., 10, 252–261.
Chandrasekhar, S. (1943), Stochastic problems in physics and astronomy,Rev. Mod. Phys., 15, 1–89.
Chapman, S., and T. G. Cowling (1991), The Mathematical Theory of Nonuniform Gases, Cambridge University Press, Cambridge, UK, p. 447.
Cherrey, K. D., M. Flury, and J. B. Harsh (2003), Nitrate and colloid transport through coarse Hanford sediments under steady state, variably saturated flow, Water Resour. Res., 39(6), 1165, doi:10.1029/2002WR
001944.
Choi, N. C., D. J. Kim, and S. B. Kim (2007), Quantification of bacterial
mass recovery as a function of pore-water velocity and ionic strength,Res. Microbiol., 158, 70–78.
Corapcioglu, M. Y., and H. Choi (1996), Modeling colloid transport in unsaturated porous media and validation with laboratory column data,
Water Resour. Res., 32, 3437–3449, doi:10.1029/95WR02693.
DeFlaun, M. F., C. J. Murray, M. Holben, T. Scheibe, A. Mills, T. Ginn,
T. Griffin, E. Majer, and J. L. Wilson (1997), Preliminary observations
on bacterial transport in a coastal plain aquifer, FEMS Microbiol. Rev.,20(3–4), 473–487.
Derjaguin, B. V., and L. D. Landau (1941), Theory of the stability of strongly
charged lyophobic sols and of the adhesion of strongly charged particles
in solutions of electrolytes, Acta Physicochim. USSR, 14, 733–762.
Dong, H., T. C. Onstott, C.-H. Ko, A. D. Hollingsworth, D. G. Brown, andB. J. Mailloux (2002), Theoretical prediction of collision efficiency
between adhesion-deficient bacteria and sediment grain surface, Colloid
Surf. B: Biointerfaces, 24, 229–245.
Duffadar, R. D., and J. M. Davis (2007), Interaction of micrometer-scaleparticles with nanotextured surfaces in shear flow, J. Colloid Interface
Sci., 308, 20–29.
Duffadar, R. D., and J. M. Davis (2008), Dynamic adhesion behavior of micrometer-scale particles flowing over patchy surfaces with nanoscale
electrostatic heterogeneity, J. Colloid Interface Sci., 326, 18–27.
Duffadar, R. D., S. Kalasin, J. M. Davis, and M. M. Santore (2009), The
impact of nanoscale chemical features on micron-scale adhesion: Crossover from heterogeneity-dominated to mean field behavior, J. Colloid
Interface Sci., 337, 396–407.
Elimelech, M., J. Gregory, X. Jia, and R. A. Williams (1995), in Particle
Deposition & Aggregation Measurement, Modeling and Simulation, Butterworth-Heinemann, Woburn, MA.
Elimelech, M., J. Y. Chen, and Z. A. Kuznar (2003), Particle deposition
onto solid surfaces with micropatterned charge heterogeneity: The‘‘hydrodynamic bump’’ effect, Langmuir, 19, 6594–6597.
Franchi, A., and C. R. O’Melia (2003), Effects of natural organic matterand solution chemistry on the deposition and reentrainment of colloids inporous media, Environ. Sci. Technol., 37(6), 1122–1129.
Gaillard, J. F., C. Chen, S. H. Stonedahl, B. L. T. Lau, D. T. Keane, and A.
I. Packman (2007), Imaging of colloidal deposits in granular porous
media by X-ray difference micro-tomography, Geophys. Res. Lett., 34,
L18404, doi:10.1029/2007GL030514.
Gargiulo, G., S. A. Bradford, J. Simunek, P. Ustohal, H. Vereecken, and E.
Klumpp (2007), Bacteria transport and deposition under unsaturated 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. Simunek, P. Ustohal, H. Vereecken, and E.
Klumpp (2008), Bacteria transport and deposition under unsaturated flow
conditions: The role of water content and bacteria surface hydrophobicity, Vadose Zone J., 7, 406–419.
Hahn, M. W., D. Abadzic, and C. R. O’Melia (2004), Aquasols: On the roleof secondary minima, Environ. Sci. Technol., 38, 5915–5924.
Harvey, R. W., and S. P. Garabedian (1991), Use of colloid filtration theory
in modeling movement of bacteria through a contaminated sandy aquifer,
Environ. Sci. Technol., 25, 178–185.
Hoek, E. M. V., and G. K. Agarwal (2006), Extended DLVO interactions
between spherical particles and rough surfaces, J. Colloid Interface Sci.,
298, 50–58.
Hubbe, M. A. (1984), Theory of detachment of colloidal particles from platsurfaces exposed to flow, Colloid Surf., 12, 151–178.
Israelachvili, J. N. (1992), in Intermolecular and Surface Forces, 2nd ed.,Academic, London.

W10503

Johnson, K. L., K. Kendall, and A. D. Roberts (1971), Surface energy andthe contact of elastic solids, Proc. R. Soc. London Ser. A, 11, 301–313.
Johnson, P. R., and M. Elimelech (1995), Dynamics of colloid deposition in
porous media: blocking based on random sequential adsorption, Langmuir, 11, 801–812.
Johnson, W. P., X. Li, and S. Assemi (2007), Deposition and re-entrainment
dynamics of microbes and non-biological colloids during non-perturbed
transport in porous media in the presence of an energy barrier to deposition, Adv. Water Resour., 30, 1432–1454.
Kalasin, S., and M. M. Santore (2008), Hydrodynamic crossover indynamic microparticle adhesion on surfaces of controlled nanoscale heterogeneity, Langmuir, 24, 4435–4438.
Kim, H. N., S. A. Bradford, and S. L. Walker (2009), Escherichia coli
O157:H7 transport in saturated porous media: Role of solution chemistry
and surface macromolecules, Environ. Sci. Technol., 43(12), 4340–4347.
Kozlova, N., and M. M. Santore (2006), Micron-scale adhesion dynamics
mediated by nanometer-scale surface features, Langmuir, 22, 1135–1142.
Kozlova, N., and M. M. Santore (2007), Micrometer scale adhesion on
nanometer-scale patchy surfaces: Adhesion rates, adhesion thresholds,
and curvature-based selectivity, Langmuir, 23, 4782–4791.
Kretzschmar, R., and H. Sticher (1998), Colloid transport in natural porous
media: Influence of surface chemistry and flow velocity, Phys. Chem.
Earth, 23, 133–139.
Kuznar, Z. A., and M. Elimelech (2007), Direct microscopic observation ofparticle deposition in porous media: Role of the secondary energy minimum, Colloid Surf. A: Physicochem. Eng. Aspects, 294, 156–162.
Leij, F. J., and S. A. Bradford (2009), Combined physical and chemical
nonequilibrium transport model: Analytical solution, moments, andapplication to colloids, J. Contamin. Hydrol., 110, 87–99.
Li, X., and W. P. Johnson (2005), Non-monotonic variations in removal
rate coefficients of microspheres in porous media under unfavorable deposition conditions, Environ. Sci. Technol., 39, 1658–1665.
Li, X., T. D. Scheibe, and W. P. Johnson (2004), Apparent decreases in colloid deposition rate coefficient with distance of transport under unfavorable deposition conditions: A general phenomenon, Environ. Sci. Technol.,38, 5616–5625.
Li, X., P. Zhang, C. L. Lin, and W. P. Johnson (2005), Role of hydrodynamic
drag on microsphere deposition and re-entrainment in porous media underunfavorable conditions, Environ. Sci. Technol., 39, 4012–4020.
Li, X., C.-L. Lin, J. D. Miller, and W. P. Johnson (2006), Pore-scale observation of microsphere deposition at grain-to-grain contacts over assemblage-scale porous media domains using X-ray microtomography,
Environ. Sci. Technol., 40, 3762–3768.
Lindqvist, R., J. S. Cho, and C. G. Enfield (1994), A kinetic model for celldensity dependent bacterial transport in porous media, Water Resour.
Res., 30, 3291–3299, doi:10.1029/94WR01725.
Liu, D., P. R. Johnson, and M. Elimelech (1995), Colloid deposition dynamics in flow-through porous media: role of electrolyte concentration,
Environ. Sci. Technol., 29, 2963–2973.
Marquardt, D. W. (1963), An algorithm for least squares estimation of nonlinear parameters, SIAM J., 11, 431–441.
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.
Redman, J. A., S. B. Grant, T. M. Olson, and M. K. Estes (2001), Pathogen
filtration, heterogeneity, and the potable reuse of wastewater, Environ.
Sci. Technol., 35, 1798–1805.
Ryan, J. N., and M. Elimelech (1996), Colloid mobilization and transport ingroundwater, Colloid Surf. A: Physicochem. Eng. Aspects, 107, 1–56.
Schijven, J. F., and S. M. Hassanizadeh (2000), Removal of viruses by soil
passage: Overview of modeling, processes, and parameters, Crit. Rev.
Environ. Sci. Technol., 30, 49–127.
Schijven, J. F., S. M. Hassanizadeh, and H. de Bruin (2002), Two-site kinetic modeling of bacteriophages transport through columns of saturated
dune sand, J. Contamin. Hydrol., 57, 259–279.
Shapiro, A. A., and P. G. Bedrikovetsky (2010), A stochastic theory fordeep bed filtration accounting for dispersion and size distributions, Physica A, 389, 2473–2494.
Shen, C., B. Li, Y. Huang, and Y. Jin (2007), Kinetics of coupled primaryand secondary-minimum deposition of colloids under unfavorable chemical conditions, Environ. Sci. Technol., 41, 6976–6982.
Simoni, S. F., H. Harms, T. N. P. Bosma, and A. J. B. Zehnder (1998),
Population heterogeneity affects transport of bacteria through sand columns at low flow rates, Environ. Sci. Technol., 32, 2100–2105.
Simůnek,J., and M. Th. van Genuchten (2008), Modeling nonequilibrium
flow and transport using HYDRUS, Vadose Zone J., 7, 782–797.

11 of 12

W10503

BRADFORD ET AL.: MODELING COLLOID TRANSPORT AND RETENTION
Simůnek,J., M. Th. van Genuchten, and M. Sejna(2008), The HYDRUS1D software package for simulating the one-dimensional movement ofwater, heat, and multiple solutes in variably saturated media, Version
4.0, HYDRUS Software Series 3, Department of Environmental Sciences, University of California Riverside, Riverside, CA, p. 315.
Song, L., and M. Elimelech (1994), Transient deposition of colloidal
particles in heterogeneous porous media, J. Colloid Interface Sci., 167,
301–313.
Tan, Y., J. T. Gannon, P. Baveye, and M. Alexander (1994), Transport ofbacteria in an aquifer sand: Experiments and model simulations, Water
Resour. Res., 30, 3243–3252, doi:10.1029/94WR02032.
Tong, M., X. Li, C. N. Brow, and W. P. Johnson (2005), Detachmentinfluenced transport of an adhension-deficient bacterial strain withinwater-reactive porous media, Environ. Sci. Technol., 39, 2500–2508.
Tong, M., H. Ma, and W. P. Johnson (2008), Funneling of flow into grainto-grain contacts drives colloid-colloid aggregation in the presence of anenergy barrier, Environ. Sci. Technol., 42, 2826–2832.
Torkzaban, S., S. A. Bradford, and S. L. Walker (2007), Resolving thecoupled effects of hydrodynamics and DLVO forces on colloid attachment to porous media, Langmuir, 23, 9652–9660.
Torkzaban, S., S. S. Tazehkand, S. L. Walker, and S. A. Bradford (2008),
Transport and fate of bacteria in porous media: Coupled effects of chemical conditions and pore space geometry, Water Resour. Res., 44,
W04403, doi:10.1029/2007WR006541.
Tufenkji, N., and M. Elimelech (2004), Correlation equation for predicting
single-collector efficiency in physiochemical filtration in saturated porous media, Environ. Sci. Technol., 38, 529–536.
Tufenkji, N., and M. Elimelech (2005a), Breakdown of colloid filtration
theory: Role of the secondary energy minimum and surface charge heterogeneities, Langmuir, 21, 841–852.

W10503Tufenkji, N., and M. Elimelech (2005b), Spatial distributions of Cryptosporidium oocysts in porous media: Evidence for dual mode deposition,
Environ. Sci. Technol., 39, 3620–3629.
Verwey, E. J. W., and J. Th. G. Overbeek (1948), Theory of the Stability of
Lyophobic Colloids, Elsevier, Amsterdam.
Wan, J. M., and T. K. Tokunaga (2002), Partitioning of clay colloids at air–water interfaces, J. Colloid Interface Sci., 247, 54–61.
Xu, S., B. Gao, and J. E. Saiers (2006), Straining of colloidal particles insaturated porous media, Water Resour. Res., 42, W12S16, doi:10.1029/
2006WR004948.
Yao, K. M., M. T. Habibian, and C. R. O’Melia (1971), Water and waste
water filtration—Concepts and applications, Environ. Sci. Technol., 5,1105–1112.
Yoon, J. S., J. T. Germaine, and P. J. Culligan (2006), Visualization of particle behavior with a porous medium: Mechanisms for particle filtration
and retardation during downward transport, Water Resour. Res., 42,
W06417, doi:10.1029/2004WR003660.
Yuan, H., and A. A. Shapiro (2011), A mathematical model for non-monotonic deposition profiles in deep bed filtration systems, Chem. Eng. J.,
166, 105–115.
Zhang, P., W. P. Johnson, T. D. Scheibe, K. Choi, F. C. Dobbs, and B. J.
Mailloux (2001), Extended tailing of bacteria following breakthrough at
the Narrow Channel Focus Area, Oyster, Virginia, Water Resour. Res.,37, 2687–2698, doi:10.1029/2000WR000151.

S. A. Bradford, US Salinity Laboratory, USDA, ARS, 450 West Big
Springs Road, Riverside, CA 92507, USA. (Scott.Bradford@ars.usda.gov)
J. Simunek and S. Torkzaban, Department of Environmental Sciences,
University of California, 2320 Geology, Riverside, CA 92521, USA.

12 of 12