Special Section: Reactive
Transport Modeling

Qiulan Zhang
S. Majid Hassanizadeh*
Amir Raoof
M.Th. van Genuchten
Saskia M. Roels

Remobilization of colloids occurs
under transient flow conditions.
Three different approaches were
employed to simulate observed
transient flow eﬀects on the remobiliza on of a ached viruses in sandycolumns. The approach taking into
account specific air-water interfacial
area showed be er agreement thanthe previously proposed models.
Q. Zhang, S.M. Hassanizadeh, and A. Raoof,
Dep. of Earth Sciences, Utrecht Univ., Budapestlaan 4, 3584CD Utrecht, the Netherlands;
S.M. Hassanizadeh, Soil and Groundwater
Systems, Deltares, the Netherlands; M.Th.
van Genuchten, Dep. of Mechanical Engineering, COPPE/LTTC, Federal Univ. of Rio deJaneiro, Rio de Janeiro, RJ 21945-970, Brazil; and S.M. Roels, Dep. of Geotechnology,
Del Univ. of Technology, Stevinweg 1, 2628
CN Del , the Netherlands. *Corresponding
author (S.M.Hassanizadeh@uu.nl).
Vadose Zone J.
doi:10.2136/vzj2011.0090
Received 20 July 2011.

© Soil Science Society of America
5585 Guilford Rd., Madison, WI 53711 USA.

All rights reserved. No part of this periodical may
be reproduced or transmi ed in any form or by anymeans, electronic or mechanical, including photocopying, recording, or any informa on storage
and retrieval system, without permission in wri ng
from the publisher.

Modeling Virus Transport andRemobilizaƟon during TransientParƟally Saturated Flow
Virus transport in porous media is aﬀected by the water flow regime. During transient,variably saturated flow, fluctua ng flow regimes can enhance virus detachment fromboth solid–water interfaces (SWIs) and air–water interfaces (AWIs). The objec ve of thisstudy was to simulate the influence of drainage and imbibi on events on the remobilizaon of a ached viruses. Three diﬀerent modeling approaches were examined. In the first
approach, all a achment and detachment coeﬃcients were assumed to be constant, but
the values of the detachment coeﬃcients were increased dras cally for the dura on oftransient, unsaturated flow. The second and third modeling approaches involved extensions of the model of Cheng and Saiers, who assumed enhanced detachment of viruses tobe propor onal to the me rate of change in the water content. Their model was extended
to include separate terms for virus a achment–detachment on SWIs and AWIs. In our
second approach, we assumed kine c sorp on onto the AWI, with the desorp on rate
being described as a func on of temporal changes in the air content. This approach did
not explicitly account for the specific air–water interfacial area. Thus, in our third approach
we explicitly included the presence and varia on of air–water interfaces and assumed AWI
a achment–detachment to be an equilibrium sorp on process. The available air–waterinterfacial area was assumed to be a func on of fluid satura on. The models were used
to simulate a series of saturated–unsaturated virus transport experiments reported in theliterature for condi ons of both drainage and imbibi on. The most promising results wereobtained with the third approach, which explicitly accounts for adsorp on to air–waterinterfaces and assumes equilibrium sorp on on the available air–water interfacial area.
Abbrevia ons: AWI, air–water interface; SWI, solid–water interface.

Studies of virus and colloid transport during transient flow in variably saturated porous media are important for devising optimal management and remediation
practices for polluted water resources, especially for contaminants that may adsorb onto
and migrate with colloids, including the transport of biological colloidal particles (Goyal
et al., 1984; Ouyang et al., 1996). Transient-flow regimes triggered by rainfall, irrigation,
snowmelt, or water table fluctuations predominate in the shallow vadose zone. Much evidence exists that fluctuating flows may remobilize deposited pathogenic microorganisms
and can lead to very high peak concentrations in drinking water resources (DeNovio etal., 2004; Zhuang et al., 2007). The presence of such peak events can be a major threat to
public health.
Fate and transport processes of viruses and colloids in porous media are governed byadvection, dispersion, and inactivation, as well as by interactions with various interfaces(Keswick and Gerba, 1980; Yates et al., 1987; Schijven and Hassanizadeh, 2000). In unsaturated soils, the main removal processes are inactivation and attachment to solid surfaces
and air–water interfaces (Powelson et al., 1990; Schijven and Hassanizadeh, 2000; Chu
et al., 2001; Crist et al., 2005). These processes very much depend on prevailing hydraulicand chemical conditions that are prone to spatial and temporal changes (e.g., see Sadeghi
et al. [2011] for the effect of pH and ionic strength on virus removal). Of importance tovirus and colloid transport in the unsaturated zone is also the air–water interfacial area
between the wetting (water) and nonwetting (air) phases, which is a function of the degreeof saturation and the capillary pressure (Hassanizadeh and Gray, 1993). Much theoretical,
experimental, and modeling research has recently been undertaken to understand andquantify the relationships between air–water interfacial area, capillary pressure, and fluid
saturation (e.g., Pyrak-Nolte et al., 2008; Joekar-Niasar et al., 2010b).

www.VadoseZoneJournal.org

Several studies of virus and colloid transport during variably saturated flow have shown that attachment and release are enhanced atlower degrees of saturation (e.g., Torkzaban et al., 2006a; Chu etal., 2001); however, most of these studies were performed assuming
steady-state flow with uniform degrees of saturation in time andspace. Commonly, a steady-state flow regime is established first
before the virus transport experiments are initiated. The results ofsuch studies are valuable in that they show the dependence of sorption or attachment–detachment parameters on fluid saturation;
however, natural systems usually involve transient-flow conditions,which have been shown to often cause mobilization of attachedviruses and colloids.
Gómez-Suárez et al. (1999a, 1999b, 2001b) studied the detachment of microsized particles from solid surfaces by passing air
bubbles through a parallel-plate flow chamber. They varied fluid
properties and hydraulic conditions to study their effects on colloid detachment. Detachment was found to decrease linearly withincreasing air bubble velocity and decreasing air–water interfacial
tension. The percentages of detachment by passing air bubbles
also varied greatly depending on the characteristics of the bacterial strain and the solid surface (Leenaars and O’Brien, 1989;
Gómez-Suárez et al., 2001a; Sharma et al., 2005; Boks et al., 2008).
Results indicated that at low air bubble velocities, the hydrophobicity of the substratum had no influence on detachment. At high air
bubble velocities, however, all bacterial strains were more efficiently
detached from hydrophilic glass substrata. Particle shape (GómezSuárez et al., 2001a) and size (Gómez-Suárez et al., 2001b) werealso found to be an important factor in air-bubble-induced detachment. Sharma et al. (2008a, 2008b) performed experiments similar
to Gómez-Suárez et al. (1999a). They observed that an infi ltratingor draining water front removed colloids that were deposited on
the solid surfaces whenever the colloids came into contact withthe moving AWI.
Several pore-scale visualization studies on the role of moving
AWIs on colloid transport have been performed. For example,
Sirivithayapakorn and Keller (2003) conducted pore-scale imbibition experiments in an unsaturated physical micromodel. They
found that colloids strongly attached to the AWI. In subsequentstudies, Keller and Auset (2007) showed how imbibition eventsthickened the water fi lms where colloids were originally trapped
and mobilized colloids attached to the AWI. They found that
resaturation eventually remobilized all colloids attached to the
AWI or trapped within thin water fi lms. Crist et al. (2004, 2005)
performed three-dimensional micromodel visualization experiments using a packed-sand infi ltration chamber. They found thatcolloids not only attached to the AWI but also were trapped inthin regions where the AWI and the solid surface meet, often
referred to as the air–water–solid contact line. Other pore-scale
experiments on colloid transport, including retention and mobilization in unsaturated media during imbibition, were performed
by Gao et al. (2006). As the water content increased and immobile

water was converted to mobile water, particles were found to bereleased from the immobile water zones.
Several drainage and imbibition column experiments have been
performed to study the remobilization of colloids and viruses asa result of transient, unsaturated flow. Results showed that theconcentration of viruses at the outlet increased when a wetting
front arrived (Ryan et al., 1998; Powelson and Mills, 2001; Saiersand Lenhart, 2003; Torkzaban et al., 2006a, 2006b; Shang et al.,2008; Zhuang et al., 2009; Cheng and Saiers, 2009). These studies suggest that temporal changes in saturation, whether drainage
or imbibition, promote the remobilization of attached viruses
and colloids.
To account for colloid remobilization during transient flow, Cheng
and Saiers (2009) proposed a semi-empirical model in which thedetachment coefficient is assumed to be a function of the temporal change in fluid saturation. They used the model to simulate
the results of their colloid remobilization experiments. Unfortunately, the model has two shortcomings: (i) the dependency of theattachment and detachment coefficients on fluid saturation was
not taken into account, and (ii) attachment to and detachmentfrom AWIs was not specifically included in the model.
The objective of this study was to investigate the utility of
the Cheng and Saiers model, and several modifications andextensions thereof, to describe experimental data reported by
Torkzaban et al. (2006a, 2006b). We modified the model ofCheng and Saiers (2009) in an attempt to remove the two shortcomings mentioned above.

6 Column Experiments
The virus transport experiments of Torkzaban et al. (2006a, 2006b)
were designed to determine the effect of water content on the attachment of viruses in soil during steady-state, unsaturated flow. They
focused especially on the interaction of viruses with the SWI and
AWI. In their studies, they used cylindrical polyvinyl chloride columns with an internal diameter of 10 cm and length of 23 cm. Sand
with a median grain size (d50) of 140 μm and uniformity coefficient (d60/d10) of 1.6 was packed into the columns following theprocedures of Robinson and Friedman (2001) to produce as homogeneous a packing as possible. Both saturated and unsaturated flow
experiments at constant flow rates and constant saturations were
performed. The pH and ionic strength were also varied among theexperiments. Table 1 lists the relevant parameters of the experiments.
The experiments were designed such that constant saturation levels
and capillary pressures existed along the column and vs. time. Thiswas done by establishing zero pressure gradients in the column suchthat gravity was the only driving force for water flow.
Once steady-state flow was established, a virus suspension wasadded to the inlet solution. Effluent samples were taken from thewww.VadoseZoneJournal.org

Table 1. Experimental condition of the column experiments and modeling parameters (average ± standard deviations) as reported by Torkzaban et al.
(2006a, 2006b).
ParameterExp. LpHi100

Exp. HpLi65

Imbibition–drainage

pH

7

9

6.2

Ionic strength, mmol L−1

19

0.6

19

Initial saturation, %

100

65

50

Pore water velocity, cm min−1

0.68

0.56

0.46

Measured dispersivity, cm

0.06

0.44

0.56

Seeding duration, pore volumes

4.2

8.0

10.4

Column length, cm

23

23

23

Porosity (θ s), –

0.41

0.41

0.41

Saturated hydraulic conductivity, cm min−1

1.13

1.13

1.13

Soil dry bulk density (ρ b), kg cm−3

1650

1650

1650

0.055

0.055

0.055

8.2 (±1.7) × 10−30.82 (±0.34) × 10−3

0.2 (±0.06) × 10−3

1.4 (±0.2) × 10−3

7.0 (±1.4) × 10−3

9.5 (±1.6) × 10−3–

8.5 (±1.4) × 10−3

1.5 (±0.06) × 10−3–0.0052 (±0.001) × 10−3

0.0063 (±0.001) × 10−3

Inactivation rate (μ w, μ s, and μ a), min−1

1.3 × 10−5

1.3 × 10−5

1.5 × 10−5

Transient condition

draineddrainedresaturated, then drained

Residual water content (θ r)
Solid–water interface (SWI) attachment coefficient (kattSWI detachment coefficient (kdets), min−1s), min−1

Air–water interface (AWI) attachment coefficient (katt
AWI detachment coefficient (kdeta), min−1a), min−1

outlet of the column at regular intervals and analyzed for virus
concentrations. The virus seeding was continued until almost
steady-state breakthrough curves were obtained. At the end of theseeding period, the flow of clean water and measurement of the
breakthrough curve continued. Some of the steady-state experiments were followed by a transient drainage or imbibition flow
phase. For transient drainage, the inflow of water at the top of thecolumns was terminated and the columns were allowed to drain
by gravity to nearly residual water saturation. The unsaturated
columns in several other transient experiments were graduallyresaturated to almost full saturation by flooding them frombelow, after which they were drained to the residual water content. Both types of transient experiments caused steep increases inthe breakthrough concentrations of the viruses, as shown by typical breakthrough curves in Fig. 1 and 2. Torkzaban et al. (2006a,
2006b) used the HYDRUS-1D code (Šimůnek et al., 1998) tomodel the virus experimental data during steady-state flow. They
did not try to simulate the transient part of their experiments.

6 MathemaƟcal FormulaƟons
Flow in Unsaturated Porous Media

Transient flow of water through the unsaturated columns wasmodeled using the standard Richards equation:⎛ ∂h⎞⎤∂θ ∂ ⎡⎢ K ( h )⎜⎜ +1⎟⎟⎥=⎝⎜ ∂ x⎠⎟⎦⎥∂ t ∂ x ⎣⎢[1]

where h [L] is the water pressure head, θ (dimensionless) is thevolumetric water content, K [L T−1] is the hydraulic conductivity,
t [T] is time, and x [L] is the vertical coordinate, taken positive inthe upward direction. We used the expressions of van Genuchten(1980) to describe the nonlinear relationships between θ and h andbetween K and h:
h =−

1/ n
1 − 1/m
Se
−1)
(
α

2
⎡
1/ m m ⎤K ( h ) = K s S 1/2⎥)e ⎢1−(1− S e
⎣⎢
⎦⎥

Se =

θ−θ r
θ s −θ r

[2a]

[2b]

[3]

where θ s and θ r (both dimensionless) are the saturated and residual
water contents, respectively, Se (dimensionless) is the effective saturation, and α [L−1], n (dimensionless), and m (dimensionless) areempirical parameters, with m = 1 − 1/n.

Virus Transport during Transient,
Unsaturated Flow
We investigated three alternative models of virus transport duringtransient, variably saturated flow. The general equations governing virus transport in variably saturated media, including terms
accounting for sorption onto SWIs and AWIs, as stated by Torkzaban et al. (2006b), among others, arewww.VadoseZoneJournal.org

Fig. 2. Measured and fitted φX174 normalized virus concentration
breakthrough curve for the steady-state part of the imbibition–drainage
experiment (at 50% saturation) of Torkzaban et al. (2006a, 2006b).

rate coefficients, respectively, of viruses to and from the SWI; andkatta and kdeta [T−1] are similar attachment and detachment rate
coefficients, respectively, associated with the AWI. According to
Schijven and Hassanizadeh (2000), the ratio of μ w and μ s is generally close to unity.
Fig. 1. Normalized virus concentration breakthrough curves at pH 7
and 100% fluid saturation (top) and pH 9 and 65% fluid saturation(bottom). Figures constructed based on data and parameter values
reported by Torkzaban et al. (2006a).

∂θ C w=∂t∂ ⎛⎜ ∂ C w ⎞⎟ ∂ qC w
−μ w θ C w −γ s −γ a⎟−⎜θD
∂ x ⎜⎝∂ x ⎠⎟∂x[4a]∂ρ bC s= γ s −μ s ρ bC s∂tss
= θ k att
C w − k det
ρ bC s −μ s ρ bC s

[4b]∂θ aC a= γ a −μ a θ aC a∂taa
= θ k att
C w − k det
θ aC a −μ a θ aC a

Torkzaban et al. (2006a, 2006b) determined values of the attachment and detachment coefficients at different constant saturation
values. Those data cannot be used to provide a dependence of
the attachment/detachment coefficients on saturation duringtransient flow, however. As explained above, various researchers
have reported that during both drainage and imbibition a remobilization of attached viruses may occur. One way to model thisremobilization, using Eq. [4a–4c], is to assign a much larger, yetconstant, value to kdet for the duration of a drainage or imbibition
process. We shall refer to this as the constant-detachment model. A
more physically based model was proposed by Cheng and Saiers(2009). This model was extended in two different ways to include
the effect of AWIs.

Model 1: The Cheng and Saiers Model[4c]where Cw [pfu L−3] is the concentration of viruses in water(number of plaque-forming units per volume of water); Cs [pfuM−1] is the number of viruses adsorbed to the SWI per unit massof dry soil; Ca [pfu L−3] is the concentration of viruses adsorbed
to the AWI in terms of the number of viruses per unit volume
of air; ρ b [M L−3] is the soil bulk density; θ a (dimensionless) isthe air content (volume of air per unit volume of the soil); D [L2
T−1] is the dispersion coefficient, q [L T−1] is the Darcy–Buckingham flow rate; γ s and γ a [pfu L−3T−1] are rates of adsorption to
SWIs and AWIs, respectively; μ w, μ s, and μ a [T−1] are inactivation
rate coefficients associated with the water, solid, and air phases,
respectively; katts and kdets [T−1] are attachment and detachment

Cheng and Saiers (2009) hypothesized that mobilization ofadsorbed colloids occurs only in pores that are being emptied
during drainage or are fi lled during imbibition. During drainage,
air enters the larger pores first. As the capillary pressure increases(or water pressure decreases), smaller pores subsequently undergo
drainage. The variability in pore sizes leads to a distribution ofcompartments (pores) that empty at a distribution of entry pressures, referred to as snap-off pressures by Cheng and Saiers (2009).
Once a compartment empties as its entry pressure is reached,adsorbed particles are removed from the pore walls of that compartment and enter the water phase.
If the concentration of viruses that are adsorbed within the ith
compartment is denoted by Csi (number of adsorbed viruses perunit mass of dry soil), then the total concentration adsorbed to
SWIs at any given time and position is given by

www.VadoseZoneJournal.org

Nc

C s = ∑C si

[5]

i =1

where Nc is the number of compartments. Cheng and Saiers (2009)
assumed that the desorption of viruses occurs as a consequence of
two separate processes: a steady-state part, modeled by the standard desorption coefficient kdets, and an additional part accounting
for mobilization of viruses from a drained compartment, modeled
as kdeti(t)Csi. Thus, under transient conditions, Eq. [4b] is extended
to the following form:∂ρ bC s
=
∂t

Nc

s
s
i
(t )C s i
θ k att
C w −ρ b k det
C s −ρ b μ sC s −ρ b ∑ k det

[6]

i =1

i [T−1] is the transient detachment coefficient, standing

where kdet
for detachment from the solid grains induced by a passing AWI.
This coefficient is nonzero only if the corresponding compartment
is being drained. Thus, Cheng and Saiers (2009) proposed the following dependency for kdeti:
⎪⎧⎪0

⎪
i
k det
=⎨

∂θ
⎪⎪ N dr⎪⎪⎩∂t

for h > h s ifor h ≤ h s i[7]where Ndr (dimensionless) is an empirical coefficient that quantifies the kinetics of virus mobilization during a drainage event andhsi is the entry pressure head for the ith compartment. According
to Eq. [7], detachment of viruses within a given compartment takes
place only when the water pressure head is less than the criticalwater pressure head hsi of that compartment. Once air enters apore, detachment occurs as long as flow is locally unsteady (i.e.,∂ θ/∂t ≠ 0).
The same formula is applied during imbibition, but with different
restrictions because a pore fi lls with water when the water pressure
head exceeds the critical entry (snap-off ) pressure head:
⎪⎧⎪0
i
k det
= ⎪⎨
∂θ
⎪⎪N imb⎪⎪⎩∂t

for h < hs ifor h ≥ hs i[8]where Nimb (dimensionless) is an empirical coefficient that quantifies the kinetics of virus mobilization during an imbibition event.
The coefficients Ndr and Nimb are assumed to be the same for allcompartments.
The entry pressure head distribution hsi for different compartments
can be obtained from capillary pressure head–water content curves
such as those given by Eq. [2a]. Our choice was to divide the watercontent range between [θ r, θ ini], where θ ini is the initial water content in each transient experiment, evenly into Nc compartments,

with each ith compartment then providing a corresponding valueof the critical water pressure head, hsi. This approach assumes thatviruses attach uniformly throughout the entire pore size distribution range.

Model 2: Kine c Detachment from Air–Water
Interfaces as a Func on of Air Content Changes
The model of Cheng and Saiers (2009) does not specifically accountfor attachment to or detachment from AWIs. During imbibition,
AWIs are destroyed, and viruses adsorbed onto the AWI enter thewater phase. This desorption from the air phase to the water phase
may be accounted for by redefining the sorption rate γ a in Eq. [4a]and [4c] as⎪⎧⎪0⎪γ a = ⎨ 0 ∂θ a
⎪⎪C a
⎪⎪⎩
∂t

for h < h s ifor h ≥ h s i[9]where Ca 0 is the last known concentration of AWI-adsorbed
viruses before imbibition starts. Note that Ca0 is constant in time
but varies in space.

Model 3: Equilibrium Air–Water Interface
A achment–Detachment as a Func on of Area
An alternative and somewhat simpler formulation may be obtainedby assuming that attachment–detachment at the AWI occurs as afunction of the available air–water interfacial area. The governing
equation for virus transport in the liquid phase, including adsorption to SWIs and AWIs, is then obtained by adding Eq. [4a] and[4c] to eliminate γ a, resulting in∂θ C w ∂ aC aw
+
=
∂t
∂t
∂ ⎛⎜ ∂ C w ⎞⎟ ∂ qC ws−θ k att
Cw⎟⎟−⎜⎜ θ D
⎝
⎠
∂x
∂x
∂x

[10]

s
+ρ b k det
C s −μ w θ C w −μ a aC aw

where Caw [pfu L−2] is the concentration of viruses adsorbed to the
AWI in terms of number of viruses per unit air–water interfacial
area and a [L2 L−3] is the air–water interfacial area per unit volumeof porous medium. This equation should be solved in conjunction
with Eq. [6] for sorption onto the solid phase. Next we assume
linear equilibrium partitioning of viruses between the water andair–water interface:
C aw = K Da C w[11]where KDa [L] is the equilibrium distribution coefficient for sorption onto the AWI, expressed as KDa = θkatta/akdeta . Equation[10] then reduces towww.VadoseZoneJournal.org

∂θ RC w∂ ⎛ ∂ C w ⎞⎟ ∂ qC ws= ⎜⎜ θ D−θ k att
Cw⎟−∂t∂ x ⎜⎝∂ x ⎠⎟∂x[12]s
C s −μ t θ C w+ρ b k det

where
R =1+

aK Da
θ

μt = μw +[13a]

μ a aK Daθ[13b]

Equation [12], in conjunction with Eq. [6], is applicable to bothtransient and steady-state saturated and unsaturated flow.
The specific air–water interfacial area (a) in the above equations
is very much a transient parameter. Its value may depend not only
on saturation but also on capillary pressure (Hassanizadeh andGray, 1993; Marle, 1981; Joekar-Niasar et al., 2010a). For the purpose of this study, we assumed that a depends on water saturation,
being zero when the medium is fully saturated or completely dry.
We considered several specific relationships between the interfacial area, a, and the degree of fluid saturation, Sw = θ/θ s. The best
results were obtained using the following equation:a = a0 (1− S w ) S w
3

do so, the sum of squared residuals between calculated and measured effluent concentrations was used as the objective (or fitness)function (Raoof and Hassanizadeh, 2010). The objective function was minimized using a genetic algorithm to allow for bothconstrained and unconstrained optimization problems (Houck etal., 1995). Genetic algorithms differ from more traditional searchalgorithms in that they work with a number of candidate solutions(a population) rather than only one solution. The algorithm begins
by creating a random initial population. Then at each step, individuals are selected at random from the current population andused to produce the next generation. Over successive generations,the population “evolves” toward an optimal solution (for detailsof the method, see Houck et al., 1995). We used in our study the
Genetic Algorithm Toolbox included in Matlab (The Mathworks).
Our main objective was to estimate values of the detachment model
parameters that would fit the measurements satisfactorily. Steadystate model parameter values, such as D, katts, kdets, katta, kdeta, μ w,
μ s, and μ a, obtained by fitting the virus breakthrough curves withthe two-site kinetic attachment–detachment model are presented
in Table 1. The parameter Nc was assigned a fi xed value of 5 following a sensitivity analysis (larger values did not materially change
the optimization results). The hydraulic parameters α and n had
to be estimated also because no water retention data were available.
For this reason, we estimated these parameters as well as Ndr fromthe virus concentrations of the drainage experiment.

[14][L2 L−3] is the specific interfacial area corresponding towhere a0
residual saturation of the medium. Equation [14] assumes thatthe effective air–water interfacial area for virus sorption increases
when saturation of the medium decreases, except at very low fluid
saturations when the liquid water films around particles may
become too thin to allow unrestricted movement of viruses in theliquid phase to and from AWIs. The latter effect is accounted forby multiplication with Sw in Eq. [14]. This equation, without thelast term (Sw), is consistent with both theoretical (e.g., Bradford
and Leij, 1997; Oostrom et al., 2001) and experimental (e.g., Kim
et al., 1999) analyses of air–water interfacial areas as a functionof fluid saturation. The data of Kim et al. (1999) for a coarse-textured sand, similar to the one used in this study, indicate valuesof approximately 3 for the exponent in Eq. [14]. The value of 400
cm−1 for a0, which we used in our study, was determined by using
the simple power function a = 400(1 − Sw)3, which is physicallyrealistic to match the data of Kim et al. (1999).

6 Results

The van Genuchten hydraulic parameters in Eq. [2a] and [2b]for a very coarse-textured soil, in our case a medium sand witha median grain size of 140 μm and a uniformity coefficient of
1.6, are generally between 2 and 15 for n and between 14 and 35
m−1 for α (e.g., Schaap et al., 2001). These values were specified
as limiting values in the optimization process to ensure that theparameter values remained physically realistic. The genetic algorithm repeatedly modified a population of individual solutions.
In our case, each individual solution consisted of three parametervalues. We used a population size of 20 individual solutions anda repetition of 10 generations.
Optimized values of the hydraulic parameters α and n and the drainage detachment constant Ndr are listed in Table 2. For imbibition,
we used the same value of n as for drainage but shifted the drainage
curve by increasing the value of α by a factor of 2 (e.g., Luckner et al.,1989). With known unsaturated soil hydraulic parameters, only thetransient imbibition detachment constant Nimb had to be optimized.
The estimated value for Nimb is included in Table 2.

SimulaƟons of the Virus Transport Experiments

Parameter OpƟmizaƟon

Simula ons of the Drainage Experiments

The three different transport models presented above were solved
by means of COMSOL Multiphysics 3.5a. Model parameters were
estimated by fitting calculated concentration breakthrough curves
to the experimental data of Torkzaban et al. (2006a, 2006b). To

Figure 1 presents normalized virus concentration breakthroughcurves for the column experiments conducted by Torkzaban et al.
(2006a) at two different fluid saturations. The solid line is the simulated breakthrough curve obtained with Eq. [4a–4b] as fitted towww.VadoseZoneJournal.org

Table 2. Best-fit soil hydraulic parameter values, including retention curve parameters a and n, empirical coefficients that quantify the kinetics of virus
mobilization during imbibition and drainage (Nimb and Ndr, respectively), the specific interfacial area corresponding to residual saturation (a0), theequilibrium distribution coefficient for sorption onto the air–water interface (KDa), and the sum of squared residual (SSQ) for different experiments.

Experiment

Retention curve parameters

Cheng–Saiers extended model

α

Nimb

n

Ndr

m−1

a0

K Da

cm−1

m

SSQ

LpHi100, drainage

14.29

12.6

–

152.5 (±8.54)†––

0.0137

HpLi65, drainage

14.29

12.6

–

171.5 (±7.72)

400 (±5)0.04 (±0.01)

0.0632

Imbibition‡

28.5

12.6

0.15 (±0.008)

400 (±5)0.02 (±0.01)

0.0478Drainage

14.29

12.6

–

400 (±5)0.02 (±0.003)

Imbibition–drainage–
25.59 (±1.41)† Values are plus or minus the standard error of the fitted parameter estimate.
‡ Estimated from the drainage curve.

the steady-state flow part of the curve. Corresponding parameters
for the simulation are presented in Table 1. The two plots show
significant increases in the virus concentration due to drainage ofthe column at the end of the experiments.

better fits with the data than the variable detachment model ofCheng and Saiers (2009), as reflected by lower values of the sumof squared residuals between measured and calculated virus concentrations in the optimization process. Fitted parameters using
this approach are listed in Table 2. The equilibrium distribution
coefficients (2–4 cm) in our study were found to be several ordersof magnitude larger than values (about 0.015 cm) reported by Wan
and Tokunaga (1998). We do not have a good explanation for thesedifferences, except perhaps the different experimental protocols
used in the two studies (transient imbibition–drainage experiments in our study and a steady-state bubble column method of
Wan and Tokunaga, 1998). Further studies about these differences
are required.

As discussed above, virus remobilization can be modeled in several
ways: (i) increasing the otherwise constant value of the detachment coefficient, kdet, by some constant factor as soon as drainage
or imbibition starts; (ii) changing the value of kdet as a functionof changes in fluid saturation or air content, or (iii) changing
attachment–detachment as a function of the available air–waterinterfacial area. Results of the three modeling approaches are compared with the experimental virus drainage data in Fig. 3. Noticethat only data for the transient-flow period (and the ensuing remobilization) are shown, starting from the time when transient flow
was initiated. Parameter values for the three models are listed inTables 1 and 2. For the constant-coefficient model, the value of kdet
was increased from 1.38 × 10−3 to 3.9 × 10−2 min−1 in the caseof 100% fluid saturation (Exp. LpHi100) and from 7.02 × 10−3
to 6 × 10−2 min−1 in the case of 65% saturation (Exp. HpLi65) as
determined from the optimizations. The results are shown in the
top graphs of Fig. 3. It is clear that the fits are poor; even the trends
in the data are not properly simulated with this model.

The measured breakthrough curve, including remobilization, of
the imbibition–drainage experiments of Torkzaban et al. (2006a,
2006b) is shown in Fig. 2. The solid line is again the simulated
breakthrough curve using the conventional virus attachment–
detachment model based on Eq. [4a–4c] fitted to the steady-stateflow period. The column in this case was resaturated to nearly fullsaturation by raising the hanging water tube and then was drained
to residual water saturation.

Results for the variable-coefficient model of Cheng and Saiers(2009) are shown in the middle graphs of Fig. 3. This model fitted
the data much better than the constant-coefficient formulation.
The variations in the detachment coefficient kdeti with time at a
position 10 cm below the column inlet are shown in Fig. 4. The
maximum value occurred when the drainage front reached thisposition. It is evident that this maximum value will be reached atlater times at positions farther down the column.

The transient part of this experiment was simulated with thevariable-detachment model as well as the formulation based on
equilibrium partitioning between the water phase and air–waterinterfacial area during both resaturation and drainage. Optimized
values of the detachment parameters Nimb and Ndr, as well as otherparameters, are listed in Table 2. Figure 5 shows the measured andfitted breakthrough curves for the two approaches. The time isagain presented from the start of resaturation.

Results of the modeling approach assuming linear equilibrium
virus partitioning between the water phase and the available (transient) air–water interfacial area are shown in the bottom graphsof Fig. 3. Simulations using this approach provided somewhat

As can be seen from Fig. 5, the simulations assuming equilibrium
sorption onto the air–water interfacial area agree very well withthe observed data. It is encouraging to see that the simpler equilibrium sorption model works as well, or better, that the numerically

Simula ons of the Imbibi on–Drainage Experiment

www.VadoseZoneJournal.org

Fig. 3. Observed (open circles) and fitted (lines) virus remobilization breakthrough curves for two drainage experiments: Exp. LpHi100 (low pH, highionic strength, 100% fluid saturation) (left) and HpLi65 (high pH, low ionic strength, 65% fluid saturation) (right) simulations assuming constantdetachment rate (kdet) (top), variable kdet (middle), and equilibrium partitioning on the air–water interface (bottom).

more intricate variable detachment formulation. We note here thatthe formulation based on sorption onto the air–water interfacial
area requires an expression relating the specific interfacial area tothe degree of fluid saturation, or alternatively the air content. In
this study, we used Eq. [14] for this purpose, which performed
better than several other expressions we investigated for a, including simple power functions of relative fluid saturation such as a =a0(1 − Sw)p. Equation [14] suggests that, as expected, the effective
AWI for virus sorption increases when the medium desaturates,
except in very dry systems. Alternative expressions for Eq. [14] may

need to be investigated in future work, especially for medium- andfine-textured sand.

6 Concluding Remarks
Virus remobilization during draining column experiments andduring the processes of resaturation conducted by Torkzaban et al.
(2006a, 2006b) was modeled. First, we used the constant-detachment coefficient model, which just increased the detachment
coefficient by a constant factor. Then we modified the model of

www.VadoseZoneJournal.org

Fig. 4. Variable detachment rates(kdet) during drainage at 10 cm
below the column inlet: drainage
starting from a saturated column(left) and drainage starting froma fluid saturation of 65% (right).
Vertical bars indicate the kdet
where the compartments definedin Cheng and Saiers (2009)
model are drained.

Fig. 5. Observed (open circles) and fitted (lines) virus remobilization breakthrough curves for the imbibition–drainage experiment of Torkzaban etal. (2006a, 2006b): breakthrough results using the model with air–water interface virus attachment–detachment as a function of air content (left) andas a function of the air–water interfacial area.

Cheng and Saiers (2009) by including the effect of virus attachment on AWIs. Two alternative formulations were considered; one
was virus attachment–detachment on the AWI modeled as a function of transient changes in the air content, and the other was virus
sorption on the AWI modeled as an apparent linear equilibrium
partitioning process subject to the available (transient) air–waterinterfacial area (being zero at full saturation). The latter model
formulation was found to fit the observed data satisfactorily. Separate SWI and AWI virus mass balances (see Fig. 6 and 7) showed
that, in total, more viruses accumulated at, and released from, theSWI than the AWI.
One point needing emphasis is that the experiments of Torkzaban
et al. (2006a, 2006b) showed that steady-state attachment coefficients obtained from experiments performed at different fluid
saturations are different (see Table 1). This means that in drainage–
imbibition experiments, the steady-state attachment coefficientskatts and katta should be some function of the fluid saturation. As
such, we included in our transient simulations a linear dependenceof the attachment coefficient on saturation. Comparison of theresults with the constant-attachment model simulations, however,showed very little difference. This may have been due to the large
values of the transient detachment coefficients, which may have
obscured any variations in the attachment coefficients (causing
them to have little effect on the results).
We studied the effect of fluid saturation changes on colloid
remobilization at the column scale. As noted by Bradford et al.
(2009), among others, column scale analysis, especially if limited
to breakthrough curves, provides only limited information about
the exact processes involved. While important, such column-scalestudies may need to be augmented with more intricate pore-scale
theoretical and experimental studies to more precisely identify theunderlying processes responsible for colloid remobilization duringtransient imbibition and drainage events. Additional research isalso needed to investigate the applicability of the model of Cheng
and Saiers (2009), or the formulations developed in this study,
for virus transport in structured field soils as well as for two- or
three-dimensional flow conditions. It is further worth noting thatthe recent study of Sadeghi et al. (2011) has shown that transient

www.VadoseZoneJournal.org

Fig. 7. Concentration profiles of viruses at the solid–water interface(SWI) and the air–water interface (AWI). Dashed lines represent
concentration profiles at the very beginning of the remobilization
simulation. Solid lines represent concentration profiles at the end of
the remobilization simulation.

Acknowledgments

Fig. 6. Profiles of the amount of viruses being retained at the solid–
water interface (SWI) (ρbCs), in water (θCw), and at the air–waterinterface (AWI) (aCaw).

changes in water chemistry also result in the remobilization ofadsorbed viruses. It is important to investigate which other transient effects may have similar consequences for virus transport insoil systems.

Saeed Torkzaban (CSIRO, Adelaide, Australia) is gratefully acknowledged for providing the experimental data. We are thankful to Jack Schijven (National Instituteof Public Health and the Environment, the Netherlands) and Vahid Joekar-Niasar(Utrecht University, the Netherlands) for very helpful discussions. Financial help for
Q. Zhang by the China Scholarship Council is greatly appreciated. We are thankful
to Sabine Goldberg and three anonymous reviewers for their critical comments andvaluable recommendations that helped to improve our manuscript. The first three authors are members of the International Research Training Group NUPUS, financedby the German Research Foundation (DFG) and Netherlands Organization for Scientific Research (NWO).

References
Boks, N.P., W. Norde, H.C. van der Mei, and H.J. Busscher. 2008. Forces involvedin bacterial adhesion to hydrophilic and hydrophobic surfaces. Microbiology 154:3122–3133. doi:10.1099/mic.0.2008/018622-0
Bradford, S.A., and F.J. Leij. 1997. Es ma ng interfacial areas for mul fluid soil systems. J. Contam. Hydrol. 27:83–105. doi:10.1016/S01697722(96)00048-4

www.VadoseZoneJournal.org

Bradford, S.A., S. Torkzaban, F.J. Leij, J. Šimůnek, and M.Th. van Genuchten.
2009. Modeling the coupled eﬀects of pore space geometry and velocity on colloid transport and reten on. Water Resour. Res. 45:W02414.
doi:10.1029/2008WR007096
Cheng, T., and J.E. Saiers. 2009. Mobiliza on and transport of in situ colloids
during drainage and imbibi on of par ally saturated porous media. Water
Resour. Res. 45:W08414. doi:10.1029/2008WR007494
Chu, Y., Y. Jin, M. Flury, and M.V. Yates. 2001. Mechanisms of virus removal during transport in unsaturated porous media. Water Resour. Res. 37:253–263.
doi:10.1029/2000WR900308
Crist, J.T., J.F. McCarthy, Y. Zevi, P. Baveye, J.A. Throop, and T.S. Steenhuis. 2004.
Pore-scale visualiza on of colloid transport and reten on in partly saturated porous media. Vadose Zone J. 3:444–450.
Crist, J.T., Y. Zevi, J.F. McCarthy, J.A. Troop, and T.S. Steenhuis. 2005. Transportand reten on mechanisms of colloids in par ally saturated porous media.
Vadose Zone J. 4:184–195.
DeNovio, N.M., J.E. Saiers, and J.N. Ryan. 2004. Colloid movement in unsaturated porous media: Recent advances and future direc ons. Vadose Zone
J. 3:338–351.
Gao, B., J.E. Saiers, and J.N. Ryan. 2006. Pore-scale mechanisms of colloid deposi on and mobiliza on during steady and transient flow
through unsaturated granular media. Water Resour. Res. 42:W01410.
doi:10.1029/2005WR004233
Gómez-Suárez, C., H.J. Busscher, and H.C. van der Mei. 2001a. Analysis ofbacterial detachment from substratum surfaces by the passage of air–liquid interfaces. Appl. Environ. Microbiol. 67:2531–2537. doi:10.1128/
AEM.67.6.2531-2537.2001
Gómez-Suárez, C., J. Noordmans, H.C. van der Mei, and H.J. Busscher. 1999a.
Detachment of colloidal par cles from collector surfaces with diﬀerent electrosta c charge and hydrophobicity by a achment to air bubbles in a parallel plate flow chamber. Physiol. Chem. Phys. 1:4423–4427.
doi:10.1039/a905156b
Gómez-Suárez, C., J. Noordmans, H.C. van der Mei, and H.J. Busscher. 1999b.
Removal of colloidal par cles from quartz collector surfaces as s mulated by the passage of liquid–air interfaces. Langmuir 15:5123–5127.
doi:10.1021/la981608c
Gómez-Suárez, C., H.C. van der Mei, and H.J. Busscher. 2001b. Air bubble-induced detachment of polystyrene par cles with diﬀerent sizes from collector surfaces in a parallel plate flow chamber. Colloids Surf. A 186:211–219.
doi:10.1016/S0927-7757(00)00799-8
Goyal, S.M., B.H. Keswick, and C.P. Gerba. 1984. Viruses in groundwater beneathsewage irrigated cropland. Water Res. 18:299–302. doi:10.1016/00431354(84)90103-9
Hassanizadeh, S.M., and W.G. Gray. 1993. Thermodynamic basis of capillary pressure in porous media, Water Resour. Res. 29:3389–3405.
doi:10.1029/93WR01495
Houck, C.R., J.A. Joines, and M.G. Kay. 1995. A gene c algorithm for func on
op miza on: A Matlab implementa on. Tech. Rep. NCSU-IE-TR-95–09.
North Carolina State Univ., Raleigh.
Joekar-Niasar, V., S.M. Hassanizadeh, and H.K. Dahle. 2010a. Non-equilibrium eﬀects in capillarity and interfacial area in two-phase flow: Dynamic pore-network modelling. J. Fluid Mech. 655:38–71. doi:10.1017/
S0022112010000704
Joekar-Niasar, V., M. Prodanovic, D. Wildenschild, and S.M. Hassanizadeh.
2010b. Network model inves ga on of interfacial area, capillary pressure
and satura on rela onships in granular porous media. Water Resour. Res.
46:W06526. doi:10.1029/2009WR008585
Keller, A.A., and M. Auset. 2007. A review of visualiza on techniques of biocolloid transport processes at the pore scale under saturated and unsaturated condi ons. Adv. Water Resour. 30:1392–1407. doi:10.1016/j.advwatres.2006.05.013
Keswick, B.H., and C.P. Gerba. 1980. Viruses in groundwater. Environ. Sci. Technol. 14:1290–1297. doi:10.1021/es60171a602
Kim, H., P.S.C. Rao, and M.D. Annable. 1999. Gaseous tracer technique for es ma ng air–water interfacial areas and interface mobility. Soil Sci. Soc. Am.
J. 63:1554–1560. doi:10.2136/sssaj1999.6361554x
Leenaars, A.F.M., and S.B.G. O’Brien. 1989. Par cle removal from silicon substrates using surface-tension forces. Philips J. Res. 44:183–209.
Luckner, L., M.Th. van Genuchten, and D.R. Nielsen. 1989. A consistent set ofparametric models for the two-phase flow of immiscible fluids in the subsurface. Water Resour. Res. 25:2187–2193. doi:10.1029/WR025i010p02187
Marle, C.M. 1981. From the pore scale to the macroscopic scale: Equa ons
governing mul phase fluid flow through porous media. In: A. Verruijt and
F.B.J. Barends, editors, Flow and Transport in Porous Media: Proceedings ofEuromech 143, Del , the Netherlands. 2–4 Sept. 1981. A.A. Balkema, Rotterdam, the Netherlands. p. 57–61.
Oostrom, M., M.D. White, and M.L. Brusseau. 2001. Theore cal es ma on offree and entrapped nonwe ng–we ng fluid interfacial areas in porous media. Adv. Water Resour. 24:887–898. doi:10.1016/S0309-1708(01)00017-3

Ouyang, Y., D. Shinde, R.S. Mansell, and W. Harris. 1996. Colloid enhanced
transport of chemicals in subsurface environments: A review. Crit. Rev. Environ. Sci. Technol. 26:189–204. doi:10.1080/10643389609388490
Powelson, D.K., and A.L. Mills. 2001. Transport of Escherichia in sand columns
with constant and changing water contents. J. Environ. Qual. 30:238–245.
doi:10.2134/jeq2001.301238x
Powelson, D.K., J.R. Simpson, and C.P. Gerba. 1990. Virus transport and survival
in saturated and unsaturated flow through soil columns. J. Environ. Qual.
19:396–401. doi:10.2134/jeq1990.00472425001900030008x
Pyrak-Nolte, L.J., D.D. Nolte, D. Chen, and N.J. Giordano. 2008. Rela ng
capillary pressure to interfacial areas. Water Resour. Res. 44:W06408.
doi:10.1029/2007WR006434
Raoof, A., and S.M. Hassanizadeh. 2010. A new method for genera ng porenetwork models of porous media. Transp. Porous Media 81:391–407.
doi:10.1007/s11242-009-9412-3
Robinson, D.A., and S.P. Friedman. 2001. The eﬀect of par cle size distribu onon the eﬀec ve dielectric permi vity of saturated granular media. Water
Resour. Res. 37:33–40. doi:10.1029/2000WR900227
Ryan, J.N., T.H. Illangasekare, M.I. Litaor, and R. Shannon. 1998. Par cle andplutonium mobiliza on in macroporous soil during rainfall simula ons. Environ. Sci. Technol. 32:476–482. doi:10.1021/es970339u
Sadeghi, G., J.F. Schijven, T. Behrends, S.M. Hassanizadeh, J. Gerritse, and P.J.
Kleingeld. 2011. Systema c study of eﬀects of pH and ionic strength on
a achment of phage PRD1. Ground Water 49:12–19. doi:10.1111/j.17456584.2010.00767.x
Saiers, J.E., and J.J. Lenhart. 2003. Colloid mobiliza on and transport within
unsaturated porous media under transient-flow condi ons. Water Resour.
Res. 39:10–19. doi:10.1029/2002WR001370
Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. ROSETTA: A computer program for es ma ng soil hydraulic parameters with hierarchicalpedotransfer func ons. J. Hydrol. 251:163–176. doi:10.1016/S00221694(01)00466-8
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. doi:10.1080/10643380091184174
Shang, J., M. Flury, G. Chen, and J. Zhuang. 2008. Impact of flow rate, water content, and capillary forces on in situ colloid mobiliza on during
infiltra on in unsaturated sediments. Water Resour. Res. 44:W06411.
doi:10.1029/2007WR006516
Sharma, P., H.M. Abdou, and M. Flury. 2008b. Eﬀect of the lower boundary condi on and flota on on colloid mobiliza on in unsaturated sandy sediments.
Vadose Zone J. 7:930–940. doi:10.2136/vzj2007.0163
Sharma, P., M. Flury, and J. Zhou. 2008a. Detachment of colloids from a solid
surface by a moving air–water interface. J. Colloid Interface Sci. 326:143–
150. doi:10.1016/j.jcis.2008.07.030
Sharma, P.K., M.J. Gibcus, H.C. van der Mei, and H.J. Busscher. 2005. Influenceof fluid shear and micro bubbles on bacterial detachment from a surface.
Appl. Environ. Microbiol. 71:3668–3673. doi:10.1128/AEM.71.7.36683673.2005
Šimůnek, J., M. Sejna, and M.Th. van Genuchten. 1998. The HYDRUS-1D so ware package for simula ng the one-dimensional movement of water, heat,and mul ple solutes in variably-saturated media, Version 2.0. U.S. Salinity
Lab., Riverside, CA.
Sirivithayapakorn, S., and A.A. Keller. 2003. Transport of colloids in unsaturated porous media: A pore-scale observa on of processes during
the dissolu on of air–water interface. Water Resour. Res. 39(12):1346.
doi:10.1029/2003WR002487
Torkzaban, S., S.M. Hassanizadeh, J.F. Schijven, A.M. de Bruin, and A.M. de
Roda Husman. 2006a. Virus transport in saturated and unsaturated sandcolumns. Vadose Zone J. 5:877–885. doi:10.2136/vzj2005.0086
Torkzaban, S., S.M. Hassanizadeh, J.F. Schijven, and H.H.J.L. van den Berg. 2006b.
Role of air–water interfaces on reten on of viruses under unsaturated condi ons. Water Resour. Res. 42:W12S14. doi:10.1029/2006WR004904
van Genuchten, M.Th. 1980. A closed-form equa on for predic ng the hydraulic conduc vity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.
doi:10.2136/sssaj1980.03615995004400050002x
Wan, J., and T.K. Tokunaga. 1998. Measuring par on coeﬃcients of colloids
at air–water interfaces. Environ. Sci. Technol. 32:3293–3298. doi:10.1021/
es980228a
Yates, M.V., S.R. Yates, J. Wagner, and C.P. Gerba. 1987. Modeling virus survival and transport in the subsurface. J. Contam. Hydrol. 1:329–345.
doi:10.1016/0169-7722(87)90012-X
Zhuang, J., J.F. McCarthy, J.S. Tyner, E. Perfect, and M. Flury. 2007. In situ colloid mobiliza on in Hanford sediments under unsaturated transient flow
condi ons: Eﬀect of irriga on pa ern. Environ. Sci. Technol. 41:3199–3204.
doi:10.1021/es062757h
Zhuang, J., J.S. Tyner, and E. Perfect. 2009. Colloid transport and remobiliza onin porous media during infiltra on and drainage. J. Hydrol. 377:112–119.
doi:10.1016/j.jhydrol.2009.08.011

www.VadoseZoneJournal.org