Journal of Contaminant Hydrology 55 (2002) 113 – 135
www.elsevier.com/locate/jconhyd

Kinetic modeling of virus transport at the field scale
Jack F. Schijven a,b,*, Jiřı́ Šimůnek c,1aNational Institute of Public Health and the Environment, Microbiological Laboratory for Health Protection,
P.O. Box 1, 3720 BA Bilthoven, The NetherlandsbDelft University of Technology, Faculty of Civil Engineering and Geosciences, P.O. Box 5048,
2600 GA Delft, The NetherlandscGeorge E. Brown Jr. Salinity Laboratory, USDA, ARS, 450 West Big Springs Road, Riverside,CA, 92 507, USA
Received 18 December 2000; received in revised form 8 August 2001; accepted 5 September 2001

Abstract
Bacteriophage removal by soil passage in two field studies was re-analyzed with the goal toinvestigate differences between one- and two-dimensional modeling approaches, differences
between one- and two-site kinetic sorption models, and the role of heterogeneities in the soil
properties. The first study involved removal of bacteriophages MS2 and PRD1 by dune recharge,
while the second study represented removal of MS2 by deep well injection. In both studies, removal
was higher during the first meters of soil passage than thereafter. The software packages HYDRUS1D and HYDRUS-2D, which simulate water flow and solute transport in one- and two-dimensionalvariably saturated porous media, respectively, were used. The two codes were modified by incorporating reversible adsorption to two types of kinetic sites.
Tracer concentrations were used first to calibrate flow and transport parameters of both models
before analyzing transport of bacteriophages. The one-dimensional one-site model did not fully
describe the tails of the measured breakthrough curves of MS2 and PRD1 from the dune recharge
study. While the one-dimensional one-site model predicted a sudden decrease in virus concentrations
immediately after the peaks, measured data displayed much smoother decline and tailing. The onedimensional two-site model simulated the overall behavior of the breakthrough curves very well. The
two-dimensional one-site model predicted a more gradual decrease in virus concentrations after the
peaks than the one-dimensional one-site model, but not as good as the one-dimensional two-site
model. The dimensionality of the problem hence can partly explain the smooth decrease inconcentration after peak breakthrough. The two-dimensional two-site model provided the best
results. Values for katt2 and kdet2 could not be determined at the last two of four monitoring wells,
thus suggesting that either a second type of kinetic sites is present in the first few meters of dunepassage and not beyond the second monitoring well, or that effects of soil heterogeneity and*Corresponding author. Tel.: +31-30-274-2994; fax: +31-30-274-4434.
E-mail addresses: Jack.Schijven@rivm.nl (J.F. Schijven), Jsimunek@ussl.ars.usda.gov (J. Šimůnek).
1
Tel.: +1-909-369-4865; fax: +1-909-342-4964.
0169-7722/02/$ - see front matter D 2002 Elsevier Science B.V. All rights reserved.
PII: S 0 1 6 9 - 7 7 2 2 ( 0 1 ) 0 0 1 8 8 - 7

114

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

dimensionality of the problem overshadowed this process. Variations between single collector
efficiencies were relatively small, whereas collision efficiencies varied greatly. This implies that thenonlinear removal of MS2 and PRD1 is mainly caused by variations in interactions between grain
and virus surfaces rather than by physical heterogeneity of the porous medium.
Similarly, a two-site model performed better than the one-site model in describing MS2 concentrations for the deep well injection study. However, the concentration data were too sparse in this study
to have much confidence in the fitted parameters. D 2002 Elsevier Science B.V. All rights reserved.
Keywords: Virus removal; Virus transport; MS2; PRD1; Water treatment; Artificial recharge

1. Introduction
Much interest exists, worldwide, in removing viruses by soil passage, either forprotection of groundwater, or as treatment of surface water that is subsequently used for
drinking water. Both groundwater as well as surface water may be contaminated withpathogenic viruses from various fecal sources (e.g., Yates, 1985; Havelaar et al., 1993). A
considerable amount of research has been carried out on the processes that determine virus
removal by soil passage, both at the laboratory and field scale (see Schijven and
Hassanizadeh, 2000, for a review). Virus removal during subsurface transport is due toa complex interplay of processes, of which inactivation and adsorption are of major
importance (Yates et al., 1987). In sandy aquifers, adsorption is reversible and kineticallylimited, while likely several types of kinetic sites exist (Hassanizadeh and Schijven, 2000).
In addition, processes of advection, dispersion and dilution cause spreading of viruses andthus affect virus concentrations.
Transport and removal of viruses even in field studies is often still modeled as a onedimensional problem, assuming constant and steady-state water flow and considering only
longitudinal dispersion (e.g. Bales et al., 1997; DeBorde et al., 1999; Pieper et al., 1997;
Ryan et al., 1999; Schijven et al., 1999, 2000). In reality, field situations are generally much
more complex because of temporal and spatial variations in water flow velocities, diverging
or converging flow lines, the effects of transverse dispersion and dilution. For these reasonsa two-dimensional, or even a three-dimensional, transport model may be needed.
Another complication in understanding field-scale virus transport is that removal rates
of viruses are often higher initially than further along the transport pathway (Bales et al.,
1995; DeBorde et al., 1998; Pieper et al., 1997; Ryan et al., 1999; Schijven et al., 1999,
2000). This may be due to heterogeneity within the soil environment at various scales(Schijven and Hassanizadeh, 2000), such as spatial variations of preferential attachment
sites formed by ferric oxyhydroxides (Schijven, 2001; Schijven et al., 2000) or in the soilorganic matter content (Schijven, 2001).
In the present study, data on virus removal from two field studies (Schijven et al., 1999,
2000) were re-analyzed with the goal to investigate (1) differences between one- and twodimensional modeling approaches, (2) differences between one- and two-site kinetic
sorption models, and (3) the role of heterogeneities in subsurface hydraulic properties onvirus removal. The first study involved removal of bacteriophages MS2 and PRD1 by
dune recharge at Castricum, The Netherlands (Schijven et al., 1999), while the second
study represented removal of MS2 by deep well injection at Someren, The Netherlands

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

115(Schijven et al., 1999). Details about both studies will be given below. Bacteriophages
MS2 and PRD1, used in both field studies, represent the low adsorptive behavior ofnegatively charged viruses and therefore may represent a relatively worst case scenario interms of rapid virus transport (Schijven and Hassanizadeh, 2000). This was confirmed byresults from column experiments with dune sand from the Castricum field site (Schijven,2001), which showed that a negatively charged pathogenic virus, Coxsackievirus B4, wasremoved almost as little as MS2, while a less negatively charged virus, Poliovirus 1, wasremoved much more effectively.
For our analyses, we used the software packages HYDRUS-1D (Šimůnek et al., 1998)
and HYDRUS-2D (Šimůnek et al., 1999), which simulate water flow and solute transport
in one- and two-dimensional variably saturated porous media, respectively. The two codes
were adapted by incorporating reversible adsorption to two types of kinetic sites. Data
from the two field studies were previously analyzed using one-dimensional analyticalmodels, assuming constant flow velocities, a homogeneous soil environment, and constantreaction and transport parameters between the source of bacteriophages and several
collecting wells. In this study, we will consider the two-dimensionality of both fieldproblems, two-site kinetic sorption, and the heterogeneity of subsurface.

2. Brief description of field experiments
2.1. Removal of MS2 and PRD1 by dune recharge
Removal of bacteriophages MS2 and PRD1 during dune recharge has been studied at
Castricum, The Netherlands. Details of this field experiment and an extensive evaluationof the results are found in Schijven et al. (1999). The site consists of fine aeolian dunesands that form a phreatic aquifer around 10 m thick, with a hydraulic conductivity of
about 12 m day  1 (at a water temperature of 5 C) and an effective porosity of
approximately 0.35. Surface water infiltrates the sand through a number of parallel canals
that are about 120 m apart. Water is produced through rows of wells that are placed along
midlines of the canals. Within one of the canals, a 150-m2 (10  15 m) compartment wasconstructed using PVC sheeting along wooden poles. The distance from this compartment
to extraction well PCQ414 was 63 m, with a water travel time of approximately 40 days as
determined by salt tracing. Almost perpendicular to the bank of the canal, six monitoring
wells, W1 to W6, were installed up to a distance of 30 m along a flow line to extraction
well PCQ414 (Fig. 1). Well screens W1 to W4 were 0.25 m long and well screens W5 andW6 0.50 m long. The well screens were positioned at different depths to locate them along
a calculated flow line. Table 1 lists the exact locations of these screens.
Prior to the dosage of the water with bacteriophages, sodium chloride was added to thecompartment as a conservative salt tracer to estimate the interstitial flow velocity and the
medium dispersivity. The chloride concentration of the water in the compartment wasincreased in one step from 150 to 750 mg/l and kept at this level for exactly 7 days. This
represented an increase of electric conductivity (EC) from 890 to about 2700 mS/cm. For aperiod of 35 days, the EC was measured at 10 min intervals in the six monitoring wells bymeans of fixed EC-sensors connected to a data logger.

116

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

Fig. 1. Schematic cross-section of the dune recharge field site with a source compartment and monitoring wells
W1 to W6. The curved arrow indicates the flow path. Depth is given in meters above sea level (msl). M1 to M8
denote the different materials.

Bacteriophages were seeded only after salt concentrations were back to their background levels in order to avoid any possible density effects due to dosage of the salt tracer
on estimates of attachment/detachment of bacteriophages. Dosage of the water withphages started by emptying bottles with a phage concentration of 2.5  1013 plaqueforming units (pfu) of MS2 and 2.5  1012 pfu of PRD1 at the center of the compartment. Thiswas done in order to raise the concentration in the compartment immediately to about 108
and 107 pfu/l, respectively. Although we tried to keep the concentrations constant bypumping phages from 1000-l containers into the compartment for a period of exactly 11days, they varied by a factor of about two due to variations in power supply to the dosage
pump. Water samples from the compartment and monitoring wells were taken regularly fora period of 4 months. The field study was performed during winter, when watertemperatures were low and inactivation hence was minimal. The inactivation rate ofphages in the water phase was measured both in the field and by means of laboratory
experiments.

Table 1
Location of monitoring screens at the field site for the dune recharge experiment (Vogelaar et al., 1997)
Monitoring
screenHorizontal distance tobank of canal [m]Top of screen[m + msl]

Bottom of screen[m + msl]

Mean hydraulichead [m + msl]

W1
W2
W3
W4
W5
W6

1.40
2.65
5.30
9.05
15.95
29.00

0.89
0.3
 0.03
0.03
0.17
0.01

0.64
0.05
 0.28
 0.22
 0.33
 0.49

2.37
2.26
2.07
1.83
1.44
0.77

msl is the mean sea level.

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

117

Concentrations of MS2 and PRD1 were reduced about 3 log10 within the first 2.4 m andanother 5 log10 in a linear fashion within the following 27 m. The higher initial removal
rate was investigated in more detail in column experiments (Schijven, 2001). In that study,
soil samples were taken at the field site along the same flow line where the monitoring
screens were located. These samples were analyzed physicochemically. Grain size
distribution showed little variation amongst the soil samples. But the fraction of organiccarbon and ferric oxalate, as a measure of ferric oxyhydroxides, decreased exponentiallywith distance. Following continuous seeding of bacteriophages in 1.4 m columns filledwith these different soil materials, concentrations of free bacteriophages were measured
every 10 cm along the columns to determine their removal rate with travel time or
distance. Bacteriophage removal rates were found to be highly and positively correlatedwith soil organic carbon content and to a lesser extent with the presence of ferric
oxyhydroxides (Schijven, 2001). Hence, we believe that the high initial removal within
the first meters of dune passage is due to attachment to preferable sites formed by ferric
oxyhydroxides and/or organic carbon.
2.2. Removal of MS2 by deep well injection
Removal of MS2 by deep well injection was studied at a field site near Someren, The
Netherlands. Details of this study are given in Schijven et al. (2000). Pre-treated surface
water was injected at a depth of 280 –310 m below the surface (Fig. 2). A 50-m thick layerof clay overlies the highly permeable aquifer. The aquifer itself lies on top of another 2– 3
m thick clay layer. The aquifer consists of several sandy layers (fluvial sediments) thatdiffer in permeability. Water was abstracted at a distance of 98 m from the injection point,
from a depth interval of 278 – 298 m below the surface. Monitoring wells were located at
distances of 8 m (WP3), 12 m (WP2), 22 m (WP4) and 38 m (WP1) from the injection
well along the line of flow symmetry connecting IP2 and PP1. Samples were taken from 2m long polyvinyl chloride (PVC) screens at a depth of 310 m below the surface. The
screens were all situated within the same layer of sand that had the highest permeability(25 m day 1).
Sodium chloride was injected as a conservative tracer to estimate the interstitial flow
velocity of the injected water and the dispersivity of the porous medium. For a period ofexactly 5 days, 265 kg NaCl in 880 l was injected each day into the main injection pipeline
through which pretreated surface water was injected at a rate of 960 m3 day 1. This
resulted in a salt concentration of 0.275 kg m 3 in the injection water, with the electricalconductivity (EC) of the injection water increasing in one step from 450 to about 1000 mS/
cm. During the first 2 weeks of the tracer study, the EC was measured automatically at the
number-2 screens of IP2, WP3 (8 m) and WP2 (12 m) at 4-h intervals. EC-measurements
in samples from the number-2 screens of WP1 (38 m) and WP4 (22 m) were carried out
manually at daily and weekly intervals. Similarly as for the dune sand experiment,
bacteriophages were seeded only after sodium chloride concentrations were back to theirbackground levels. A phage concentration of 7.0  107 pfu/l of MS2 was applied to theinjection well for a period of 5 days. The inactivation rate of MS2 in the water phase wasmeasured in a sample from the injection water, and from the water collected after 8 m ofsoil passage.

118

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

Fig. 2. Schematic cross-section of a field site with injection well IP2, monitoring wells WP1 to WP4 andproduction well PP1. Samples were taken from the number-2 screens. Also shown are the redox-zones in theaquifer (Stuyfzand, 1999). * O2 = dissolved oxygen concentration declining with distance from 10 to 0 mg l  1.
NO3 = nitrate concentration declining with distance from 20 to 0 mg l  1. SO4 = sulfate concentration higher than
pretreated surface water due to pyrite oxidation (Schijven et al., 2000).

Within the first 8 m of soil passage, concentrations of MS2 were reduced by 6 log10, but
only by about 2 log10 in the following 30 m. From geochemical mass balances, it could be
deduced that within the first 8 m distance from the injection well ferric oxyhydroxides
precipitated as a consequence of pyrite oxidation, but not at larger distances. Ferric
oxyhydroxides provide positively charged patches on to which fast attachment ofnegatively charged bacteriophages may take place. The non-linear removal could therefore
be ascribed to preferable attachment of MS2 to patches of ferric oxyhydroxides that arepresent within the first 8 m from the injection point, but not thereafter (Schijven et al.,2000).

3. Theory
3.1. Conceptual model
3.1.1. Water flow
Because of the unsaturated nature of part of the transport domain of the dune recharge
experiment at Castricum, subsurface flow was simulated with the variably-saturated(Richards’) equation as follows



@h
@
A @h
A
¼
KðhÞ Kij
þ Kiz
@t
@xi
@xj

ð1Þ

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

119

where xi represents the spatial coordinate [m], with x2 = z taken here to be positive upward
for the vertical cross section (i.e., for dune recharge experiment) and x2 = y for flow in ahorizontal plane (i.e., the deep well injection experiment), t is time [day], h is the watercontent [m3 m  3], h is the pressure head [m], K is the hydraulic conductivity [m day  1],
and KijA are components of a dimensionless anisotropy tensor. Knowledge of the functional
relationships relating the water contents and the hydraulic conductivities to the pressure
heads, in addition to the initial and boundary conditions, is needed to numerically solve
Eq. (1). An advantage of the Richards’ equation is that it can be used to describe
simultaneously flow in both the saturated and unsaturated zones. Eq. (1) hence can beapplied to variably saturated water flow in the dune recharge experiment, and to fully
saturated flow in the deep well injection experiment.
3.1.2. Solute and virus transport
Transport of the conservative (NaCl) tracer is described using the convection –dispersion equation. In two-dimensional form, this equation is as follows
@hc
@
@c
@c
¼
hDij
 qi
@t
@xi
@xj
@xi

ð2Þ

where c is the tracer solution concentration [g m  3], qi is the ith component of thevolumetric flux density [m day  1], and Dij is the dispersion coefficient tensor [m2 day  1]
defined as (Bear, 1972)
hDij ¼ kL AqAdij þ ðkL  kT Þqi qj
þ hDd sdij
AqA

ð3Þ

where Dd is the molecular diffusion coefficient in free water [m2 day  1], s is a tortuosity
factor [  ], AqA is the absolute value of the Darcian fluid flux density [m day  1], dij is the
Kronecker delta function, and kL and kT are the longitudinal and transverse dispersivities[m], respectively.
Schijven et al. (1999) used a one-dimensional (1D) one-site kinetic sorption model todescribe the transport of viruses at both field sites. They assumed that both attachment anddetachment of viruses to the solid phase could be described using a first-order kinetic rate
process. The governing equations describing 1D virus transport are then as follows
h

@C
@S
@2C
@C
þ qB
¼ kL hv 2  hv
 ll hC  ls qB S
@t
@t
@x
@x

ð4Þ

@S
¼ katt hC  kdet qB S  ls qB S
@t

ð5Þ

qB

where C is the concentration of free phages [pfu m  3], S is the concentration of attachedphages [pfu kg  1], h is the water content [m3 m  3], assumed here to be equal to porosityn [  ] for the fully saturated conditions, v ( = q/h) is the average interstitial water velocity[m day  1], qB is the dry bulk density [kg m  3], katt and kdet are attachment and120

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

detachment rate coefficients, respectively [day 1]; and ll and ls are inactivation rate
coefficients of the free and attached phages, respectively [day 1]. It was shown that thisone-site sorption model could not accurately duplicate the tailing of virus breakthroughcurves measured at various monitoring wells (Schijven et al., 1999), nor the breakthroughcurves measured on the laboratory samples (Hassanizadeh and Schijven, 2000).
Hassanizadeh and Schijven (2000) showed that the breakthrough curves of bacteriophages MS2 and PRD1 in the column experiments could be better described by means ofa convection –dispersion model that incorporates adsorption to two types of kinetic sites.
Their model assumes that the sorption sites on the solid phase can be divided into two
fractions with different properties and various attachment and detachment rate coefficients.
The governing equations of the two-site sorption model are as follows (Hassanizadeh and
Schijven, 2000)
h

@C
@S1
@S2
@2C
@C
þ qB
 ll hC  ls1 qB S1  ls2 qB S2
þ qB¼ kL hv 2  hv
@t
@x
@x
@t
@t
ð6Þ

qB

@S1
¼ katt1 hC  kdet1 qB S1  ls1 qB S1
@t

ð7Þ

qB

@S2
¼ katt2 hC  kdet2 qB S2  ls2 qB S2@tð8Þwhere all variables are the same as for Eqs. (4) and (5), with the subscripts 1 and 2
referring to the two different kinetic sites.
3.2. Numerical model
The governing water flow (Eq. (1)) and solute transport (Eqs. (2) – (8)) equations weresolved numerically subject to appropriate initial and boundary conditions using the
HYDRUS-1D (Šimůnek et al., 1998) and HYDRUS-2D (Šimůnek et al., 1999) models
for the one- and two-dimensional problems, respectively. Galerkin-type linear finite
element methods were used for the spatial discretization, and a finite difference method
for temporal discretization. A fully implicit mass conservative modified Picard iteration
scheme (Celia et al., 1990) was used to solve the Richards’ equation.
HYDRUS-1D and HYDRUS-2D simulate water flow and solute transport in variably
saturated porous media. The two codes considered the concept of two-site sorption (Selim
et al., 1987; van Genuchten and Wagenet, 1989) to permit consideration of nonequilibrium
adsorption –desorption reactions. The two-site sorption concept assumes that the sorption
sites can be divided into two fractions, with sorption on one fraction of the sites (the type-1
sites) being instantaneous, and with sorption on the remaining (type-2) sites being timedependent. The HYDRUS codes hence permitted only one type of sorption sites to be
kinetic. To make the codes applicable to our study, HYDRUS-1D and HYDRUS-2D wereboth modified by implementing reversible kinetic adsorption to two sets of kinetic sites asmodeled with Eqs. (6) –(8).

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

121

3.3. Inverse modeling
Solute transport and reaction parameters, as well as saturated hydraulic conductivities
and porosities were obtained by calibrating the HYDRUS-1D and HYDRUS-2D model
solutions against the measured data. This was done by iteratively changing model
parameters and thus improving model fits to the measured data until a desired degree ofprecision is obtained. The objective function, SSQ, that was minimized during theparameter estimation procedure contains weighted residuals between the measured data
and the model predictionsSSQðb; qÞ ¼

m
X
j¼1

vj

nj
X


2
wij qj ðx; ti Þ  qj ðx; ti ; bÞ

!
ð9Þ

i¼1

where the right-hand side represents weighted deviations between the measured andcalculated space-time variables (e.g., concentrations and/or pressure heads at different
times in the flow domain). In this term, m is the number of different sets ofmeasurements, nj is the number of measurements in a particular measurement set,
qj* (x,ti) represents specific measurements at time ti for the jth measurement set at
location x,qj (x,tib) are the corresponding model predictions for the vector of optimized
parameters b (e.g., Ks, n, kL, katt and kdet), and vj and wi,j are weights associated with aparticular measurement set or point, respectively. Minimization of SSQ was accomplished using the Levenberg– Marquardt nonlinear minimization algorithm (Marquardt,
1963).

4. Analysis of the dune recharge experiment4.1. Water flow and tracer transport
Two types of analyses were carried out with data collected during the dune recharge
experiment. First, HYDRUS-1D was used to optimize the hydraulic (i.e., the saturated
hydraulic conductivity and porosity) and solute transport (i.e., longitudinal dispersivity
and attachment and detachment constants) parameters along a one-dimensional flow path,similarly as in previous studies (Schijven et al., 1999, 2000). In the initial analyses weassumed that the different parameters were constant between the source of tracer or phages
and particular wells. However, because model parameters calibrated against data measured
at one well often could not reproduce concentrations at the other wells, in later analysis weallowed various variables to vary between individual wells. For comparison purposes welist in Table 2 results obtained previously (Schijven et al., 1999) using the one-dimensional
one-site sorption model as modeled with the CXTFIT2 code (Toride et al., 1995). The
optimized dispersivities and average pore-water velocities were fitted to the EC and virus
data assuming stream tubes from the source to a particular well and with the aquifer having
homogeneous properties and constant porosity n of 0.35.

122

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

Table 2Estimated values of model parameters determined using the CXTFIT2 analytical model (one-dimensional, onesite kinetic sorption) from individual breakthrough curves of NaCl, MS2 and PRD1 for the dune recharge
experiment (Schijven et al., 1999)
Solute

NaCl
MS2

PRD1

Variable

Units

X
v
k*L
ll
ls
katt
kdet
ll
ls
katt
kdet

m
m day 1mday 1day 1day 1day 1day 1day 1day 1day 1

Wells
W1

W2

W3

W4

W5

W6

2.4
1.41
0.008
0.03
0.085
4.1
0.00087
0.12
0.071
4.0
0.00077

3.8
1.56
0.012
0.03
0.092
3.2
0.0016
0.12
0.067
3.1
0.0011

6.4
1.59
0.017
0.03
0.092
2.8
0.0026
0.12
0.067
2.2
0.0018

10.2
1.57
0.017
0.03
0.092
2.0
0.0018
0.12
0.067
1.5
0.002517.1
1.52
0.0096
0.03
0.092
1.3
0.00052
0.12
0.067
1.3
0.0021

30.1
1.19
0.08
0.03
0.092
0.8
0.0030
0.12
0.067
0.7
0.0034

X is the distance along a streamline.

Table 3 presents stream-tube dispersivities and average pore-water velocities as fittedwith the HYDRUS-1D code to the measured EC data. As expected, the optimization
results are almost identical to those in Table 2, thus verifying the numerical solution in the
HYDRUS-1D code. Optimized pore-water velocities were smaller for the first well (W1),almost constant for the next four wells, and then significantly lower for the last well (W6).
The optimized longitudinal dispersivities were very small (on the order or 1 cm) with the
exception of the last well, for which we obtained a dispersivity of 8 cm. Although
excellent agreement was obtained between measured and fitted concentrations for aparticular well (not shown), the optimized parameters could not always be used successfully to predict concentrations at other locations.
To improve the description of concentrations at all wells simultaneously, in the second
step of our one-dimensional analyses we assumed that several variables could change
between different wells. We used this simplified conceptual model of heterogeneity, i.e.,based on the location of wells, because our physical and geochemical analyses of soil

Table 3One-dimensional optimizations using EC data from the dune recharge experiment and assuming stream tubes
from a source to a particular well having constant properties (n = 0.35). Each breakthrough curve was optimized
independently using HYDRUS-1D
Well

q [m day  1]

kL [m]v [m day  1]

SSQ

W1
W2
W3
W4
W5
W6

0.49 F 0.001
0.54 F 0.001
0.55 F 0.0015
0.56 F 0.0008
0.53 F 0.0008
0.41 F 0.0008

0.0069 F 0.0005
0.010 F 0.0007
0.016 F 0.0016
0.016 F 0.001
0.0074 F 0.0004
0.083 F 0.0054

1.4
1.54
1.57
1.6
1.51
1.17

0.00200.00430.013
0.00690.030
0.033

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

123

samples from different locations did not show significant differences between samples andit appeared that the sand was physically relatively homogeneous. Thus, we did not haveevidence for a more complex physically based model of soil heterogeneity. We defined an
objective function using EC breakthrough curves from all six wells, assumed a constantwater flux in the subsurface equal to the arithmetic average of the fluxes in Table 3, usedlongitudinal dispersivities between individual wells equal to those optimized for particular
breakthrough curves, and optimized the porosities between the source and the first well,and between all subsequent two consecutive wells. Although the pore-water velocities
could vary along a stream tube depending upon the porosity, the volumetric flux q must be
constant unless sources or sinks are present. This optimization assumes that instead ofusing different stream tubes with particular pore velocities and dispersivities for each well,
the solute arriving at the first well is moved further to the second well (with different
porosities and dispersivities between the first two wells), and then to the third well, etc.
Parameters derived in this way represent effective parameters, characteristic of transportdomains between particular wells. Results of the optimization are given in Table 4 and in
Fig. 3. Notice an excellent correspondence between the measured and fitted EC data using
this simplified model for describing heterogeneity of the experimental site. We obtained alarge value of porosity between wells W5 and W6 (equal to 0.6), which causes a relatively
small value of the optimized flux for the breakthrough curve of well W6, including its latearrival time.
For the two-dimensional analyses, the transport domain was first divided into four
major material layers: the surface layer M1 (from 3.0 to 0.4 m relative to the mean see
level), two intermediate layers M2 (0.4 to  3.5 m) and M3 (  3.5 to  7.5 m), and theless permeable layer M4 (  7.5 to  12 m) at the bottom of the transport profile (Fig.
1). These four layers were determined from bore profiles (Vogelaar et al., 1997). As
boundary conditions, we used a constant head in the source pool, corresponding to theheight of the ponded water, and a constant flux along the active part of the pumping well,corresponding to the pumping rate. Next we calibrated HYDRUS-2D by optimizing thesaturated hydraulic conductivities for these boundary conditions against the mean
hydraulic heads recorded for all six wells. Defining the objective function using only
the mean hydraulic heads in the observation wells produced conductivities that predictedunrealistic arrival times for the tracer fronts at these monitoring wells. Thus, in the next
step, we defined the objective function using also values of the mean flux along thestreamline as obtained from the one-dimensional analysis of the tracer fronts. This

Table 4One-dimensional optimization using EC data from the dune recharge experiment and assuming the presence of
one stream tube from the source through all wells
Well

W1

W2

W3

W4

W5

W6

SSQ

Porosity

0.384 F0.00200.282 F0.0044

0.338 F0.00260.352 F0.00220.371 F0.00110.596 F0.00130.00115All EC data were analyzed simultaneously using HYDRUS-1D; porosities and dispersivities varied between
wells. For dispersivities, see Table 3.

124

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

Fig. 3. Measured and optimized EC data assuming a stream tube from the source through all wells of the dune
recharge experiment. All EC data were analyzed simultaneously using HYDRUS-1D; porosities and dispersivities
varied between wells.

calibration led to the saturated hydraulic conductivities for particular soil layers as given
in Table 5 for materials M1 through M4.
In order to obtain better description of the EC breakthrough curves at particularmonitoring wells, material layer M2 was further subdivided into another five submaterials(M2, M5, M6, M7 and M8), corresponding to areas between individual wells. We did thisadditional subdivision only for material layer M2 since this was the only layer for which
we had information on NaCl and phages concentrations. Material M2 was situatedbetween the left side of the transport domain and well W2, material M5 was between
wells W2 and W3, M6 between W3 and W4, M7 between W4 and W5, and material M8 to
the right of well W5 (Fig. 1). We then calibrated the porosities of each of these materials

Table 5Two-dimensional simulations of EC data from the dune recharge experiment. All EC data were analyzed
simultaneously using HYDRUS-2D
Material

M1

M2

M5

M6

M7

M8

M3

M4

Ks
[m day  1]
Porosity

9.55 F
1.84
0.388 F
0.007811.2 F
1.71
0.308 F
0.060

11.2 F
1.71
0.419 F
0.021

11.2 F
1.71
0.403 F
0.018711.2 F
1.71
0.326 F0.010011.2 F
1.71
0.502 F0.008126.0 F
4.40
0.350

2.0
0.350

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

125using all measured breakthrough curves simultaneously. The results of the two-dimensional analyses of water flow and the EC data are given in Table 5. The optimizations
produced similar agreement between measured and optimized NaCl breakthrough curves
as those given in Fig. 3 for the one-dimensional numerical analysis.
Calibrated values of the saturated hydraulic conductivity, porosity and longitudinal
dispersivity were assumed to correctly describe the solute transport system and were usedin our subsequent analyses of virus transport.
4.2. Analysis of virus transport
In the case of the dune recharge study, the measured inactivation coefficients of freebacteriophages, as well as the inactivation coefficients of attached bacteriophages,
determined from the slope of the tails of the breakthrough curves, were reported by
Schijven et al. (1999). The inactivation rate coefficients of virus attached to site 2 wereTable 6
Parameter values for MS2 breakthrough curves for the dune recharge experiment
katt1

kdet1

katt2

1D One-siteW1
4.48 F 0.170
W2 0.548 F 0.241
W32.50 F 0.105
W4 0.822 F 0.0887

0.00195 F 0.00050
3.6  10  8 F 0.0061
0.0121 F 0.00270
0.0252 F 0.0120

1D Two-siteW1
4.28 F 0.0969
W2 0.585 F 0.173
W3
2.17 F 0.0708
W4 0.827 F 0.0666

0.000862 F 0.00020
0.838 F 0.231
1  10  6c F 0.00006
0.963 F 1.55
0.003110 F 0.00096 0.0001c F 0.406
0.0150 F 0.0069
0.0001c F 1.56

2D One-siteW1
5.87 F 0.161
W2 0.415 F 0.143
W3
1.74 F 0.0756
W4 0.0001 F 1.59

0.00132 F 0.00028
0.181 F 0.0885
0.00209 F 0.00186
NA

2D Two-site
W1
4.51 F 0.278
W2 0.447 F 0.321
W3
1.57 F 3.50
W4 0.0001 F 1.05

0.00132 F 0.00042
0.00202 F 0.00538
0.00804 F 0.0209
NA

1.45 F 0.363
0.0001 F 0.486
0.291 F 6.84
0.0001 F 0.550

kdet2

0.0413 F 0.135
2.10 F 6.58
NA
NA

SSQ

ga

ab

0.03910.03860.01400.0135

0.53
0.66
0.58
0.60

1.4  10 3
8.1  10 5
5.6  10 4
2.1  10 40.01880.02690.00580.0112

0.53
0.66
0.58
0.60

1.3  10 3
8.6  10 5
4.8  10 4
2.1  10 4

0.52
0.0295d 0.37
0.42
0.0324d 0.49

1.8  10 3
1.4  10 4
6.1  10 4
3.3  10 80.215 F 0.0759
NA
0.0175d52.1 F 1086.9
NA0.0183d

0.52
0.37
0.42
0.49

1.4  10 3
1.5  10 4
5.5  10 4
3.3  10 8

One stream tube from source through all wells was considered in one-dimensional optimizations. Optimizations
were carried out sequentially from well to well.
aSingle collector efficiency.
bCollision efficiency associated with katt1.
c
Parameter fixed or constraint reached.
d
In the 2D case, breakthrough curves from wells W1 and W2 (and W3 and W4) were optimized
simultaneously.

126

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

assumed to be the same as those for site 1. All inactivation coefficients were assumed to be
known and were not further included in the inverse analysis. The observations wereassumed to be log-normally distributed and hence logarithmically transformed prior to thefitting procedure.
Application of the one-dimensional one-site kinetic sorption stream-tube transportmodel produced a reasonable fit of the breakthrough curves except for considerable
discrepancy at the end of the rising and the start of the declining limbs of the
breakthrough curves (Schijven et al., 1999). Table 2 lists the corresponding parametervalues.
The measured MS2 and PRD1 breakthrough data were subsequently analyzed withthe one- and two-dimensional HYDRUS-1D and HYDRUS-2D models using both the
one-site and two-site kinetic sorption models. Attachment and detachment coefficientswere again optimized sequentially starting with the breakthrough curve measured at wellW1, then moving forward to well W2, and so on. In several optimization runs model

Table 7
Parameter values for PRD1 breakthrough curves for the dune recharge experiment
katt1

kdet1

1D One-site
W1 4.56 F 0.153
W2 0.280 F 0.308
W3 1.57 F 0.136
W4 120 F 617

0.00217 F 0.000052
1  10  6 F 0.0129
0.0353 F 0.000892
167 F 898

1D Two-site
W1 4.45 F 0.0927
W2 0.517 F 0.247
W3 0.975 F 0.0741
W4 36.8 F 682

0.00116 F 0.00022
0.00583 F 0.00638
0.00756 F 0.00302
56.8 F 1071

2D One-site
W1 5.59 F 0.298
W2 0.471 F 0.222
W3 0.215 F 0.304
W4 0.001d0.00175 F 0.00056
0.0302 F 0.0166
0.536 F 1.47
8.1  108

2D Two-site
W1 4.90 F 0.393
W2 0.437 F 1.804
W3 0.181 F 0.296
W4 0.818 F 410.5

0.00175 F 0.00078
0.00840 F 0.0437
0.162 F 0.631
88.2 too large

katt2

0.876 F 0.243
1.42 F 1.33
0.442 F 0.460
0.101 F 0.101

kdet2

0.524 F 0.170
1.86 F 3.36
1.230 F 1.76
0.00942 F 0.0481

SSQ

ga

ab

0.051
0.062
0.027
0.062

0.30
0.37
0.33
0.33

2.5  10  3
7.41  10  5
6.2  10  4
5.4  10  2

0.026
0.04350.00153
0.0624

0.30
0.37
0.33
0.33

2.4  10  3
1.4  10  4
3.9  10  4
1.7  10  2

0.06220.0622c

0.29
0.21
0.24
0.27

3.1  10  3
2.7  10  4
1.3  10  4
6.0  10  8

0.29
0.21
0.24
0.27

2.7  10  3
2.5  10  4
1.1  10  4
4.9  10  4

0.0910c

0.961 F 0.475
0.001d F 3.55
0.001d F 0.687
0.001d F 2.05

0.282 F 0.2403
NA
NA
NA

0.01390.0517c0.103cOne stream tube from source through all wells was considered in one-dimensional optimizations. Optimizations
were carried out sequentially from well to well.
aSingle collector efficiency.
bCollision efficiency associated with katt1.
c
In the 2D case, breakthrough curves from wells W1 and W2 (and W3 and W4) were optimized
simultaneously.
dParameter fixed or constraint reached.

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

127

parameters were determined simultaneously from two corresponding breakthroughcurves. Optimized parameter values determined from the MS2 and PRD1 breakthroughcurves, using both one- and two-dimensional analyses and one- and two-site kinetic
sorption models, are given in Tables 6 and 7, respectively. Measured and optimized MS2
and PRD1 breakthrough curves for the different models and viruses are presented in
Figs. 4 – 7.
The description of all breakthrough curves using models of various dimensionality andcomplexity is reasonably good. However, as can be seen from Figs. 4 and 6, the onedimensional one-site model does not completely describe the tails of the measured
breakthrough curves. While the one-dimensional one-site model predicts a sudden
decrease of virus concentrations immediately after their peak, measured data display arelatively smooth transition of the declining limb to the more straight later part of the tail.
The one-dimensional two-site model did fit the overall shape of the breakthrough curves
very well (Figs. 4 and 6).
Figs. 5 and 7 show fits obtained using the two-dimensional model and assuming onesite sorption. The two-dimensional one-site model predicted a more gradual decrease of
the virus concentrations after the peaks than the one-dimensional one-site model, but not
as good as the one-dimensional two-site model. While the two-dimensional one-site

Fig. 4. Measured and optimized MS2 breakthrough curves for the dune recharge experiment. One-dimensional
one- and two-site models were used in sequential optimizations. Dots are observations; lines are the model fits
(dashed lines for a one-site model, full lines for a two-site model).

128

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

Fig. 5. Measured and optimized MS2 breakthrough curves for the dune recharge experiment. Two-dimensional
one- and two-site models were used in sequential optimizations. Dots are observations; lines are the model fits
(dashed lines for a one-site model, full lines for a two-site model).

model seems to follow the course of the MS2 tail reasonably well, this is less so for
PRD1. Apparently, dimensionality of the problem can partly explain the smootherappearance of the concentration after peak breakthrough. Still, the two-dimensional twosite model fitted (visually as well as based on SSQ values) the breakthrough curve
better.
With the two-site model, both one- and two-dimensional, estimates of katt2 and kdet2
for MS2 and PRD1 at W2 were highly uncertain and irregular. The value of katt1 at W2
for both MS2 and PRD1 was relatively low, thus suggesting that little adsorption
occurred in the vicinity of W2. This may explain the large uncertainty found in theestimates of kdet2 at W2. Notice also that the models fitted the maximum breakthroughconcentration at W2 very well, but generally underestimated the maximum breakthroughat W1. This is because of the skewness of the climbing limb of the breakthrough curveat W1. At W3 and W4, values for katt2 and kdet2 could not be determined with a twodimensional model. When the attachment coefficients for the second type of sorption
sites were small (often estimated with large confidence intervals) or could not be
determined altogether, then obviously also the detachment coefficients could not be
determined. This suggests that a second type of kinetic sites is either present only in thefirst few meters of dune passage (due to the presence of organic carbon), but not beyond
W2, or that the effects of soil heterogeneity and dimensionality of the problem

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

129

Fig. 6. Measured and optimized PRD1 breakthrough curves for the dune recharge experiment. One-dimensional
one- and two-site models were used in sequential optimizations. Dots are observations; lines are the model fits
(dashed lines for a one-site model, full lines for a two-site model).

overshadowed this process. The latter explanation seems more plausible since the twosite kinetic sorption process was determined independently on laboratory soil columns(Hassanizadeh and Schijven, 2000). Most of the interaction with the type-1 sites took
place before W2, and thus had the greatest effects on the shapes of the breakthroughcurves at the subsequent wells. Beyond W2 little attachment took place and the minorchanges in the shape of the breakthrough curves are then mainly due to effects ofdispersion.
Estimates of katt1 between the four different models were quite similar. For both MS2
and PRD1, the highest values were found at W1, reflecting the high initial removal ratethat occurred within the first few meters of dune recharge. The other estimates of katt1 werelower, varying between 0.0 and 2.5. These variations reflect differences in adsorption
properties of the soil and determine variations in removal rates. No attachment could bedetermined between wells W3 and W4 using a two-dimensional model for bothbacteriophages. When modeling breakthrough from the source to a particular well, thevariations in rates were less explicit (see Table 2). This is because the higher initialattachment rate was redistributed along the entire stream tube, which leads to gradually
decreasing attachment coefficients for wells further along the streamline.
Tables 6 and 7 also list corresponding single collector and collision efficiencies
associated with katt1. The single collector efficiency is the rate at which particles, in our130

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

Fig. 7. Measured and optimized PRD1 breakthrough curves for the dune recharge experiment. Two-dimensional
one- and two-site models were used in sequential optimizations. Dots are observations; lines are the model fits
(dashed lines for a one-site model, full lines for a two-site model).

case viruses, strike a soil grain divided by the rate at which particles move towards thegrain and represents the physical factors determining particle collision. The collision
efficiency is the efficiency with which collided particles remain attached to the soil grains
due to electrostatic interaction. Collision efficiency a was calculated using the following
equation (Yao et al., 1971)a¼

2dc
katt1 1
3 ð1  hÞ v g

ð10Þ

where dc is the size of the sand grains. The single collector efficiency g was calculated as
described by Martin et al. (1992). For each well, velocity at that well was used for
calculation of the single collector efficiency.
Variations between single collector efficiencies were relatively small, whereas collision
efficiencies varied greatly. This implies that the nonlinear removal of MS2 and PRD1 ismainly caused by variations in the interactions between grain and virus surfaces, rather
than by physical heterogeneities of the porous medium. This is in agreement with findings
from the column experiments (Schijven, 2001) where removal was found to be strongly
and positively correlated with such surface characteristics as soil organic carbon content
and ferric oxyhydroxide content.

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

131

5. Deep well injection experiment5.1. Analyses of water flow and tracer transport
Because of the dimensionality of this problem, only two-dimensional modeling wasconsidered. We simulated the deep well injection experiment on a horizontal quadrilateraltransport domain of 200  100 m that was perpendicular to the monitoring, pumping andinjecting wells. The transport domain was discretized into two finite element meshes (3508
and 13 136 finite elements) with a very fine spatial discretization around the injecting andpumping wells and in the area between the two wells. Coarser discretizations were used
towards the external boundaries. Wells IP2 and PP1 were simulated as circles with adiameter of 80 cm and were discretized using 15 or 25 nodes for the coarser and finer finite
element mesh, respectively. A finer finite element mesh was used for simulation of NaCltransport because the small dispersivities and the absence of sorption and inactivation led
to relatively large Peclet numbers over the transport domain, a situation that usually leads
to significant numerical oscillations. Coarser meshes could be used for the virus transport
calculations since the sorption and inactivation produced considerably smoother (dispersed) curves.
Injection and pumping rates (13 m day  1) were used as boundary conditions for water
flow, and relative concentrations of injected water for solute transport. NaCl concentrations measured at wells WP2, WP3, and WP4 were used to calibrate the longitudinal
dispersivity of the aquifer. The transverse dispersivity was assumed to be one tenth of thelongitudinal dispersivity. Results of the calibration for the optimized value of thelongitudinal dispersivity of 0.163 cm are shown in Fig. 8. Note that without any
calibration of the water flow field, the fit of the first two breakthrough curves was
excellent, while only third breakthrough curve arrived slightly later than calculated. Also
notice that the tracer not only moved in the direction of the pumping well (monitoring well
WP3), but also in the opposite direction (monitoring wells WP2 and WP4). Good
calibration of the tracer breakthrough curves at wells WP2 and WP3 should be adequate
for our subsequent virus transport analysis. The breakthrough curves for MS2 at WP1 and
WP4 were not analyzed because of the limited number of measurements at these wells.
5.2. Analysis of virus transport
Calibrated value of the longitudinal dispersivity and the independently determinedsaturated hydraulic conductivity and porosity were used as starting points of the virustransport analyses. Again, both one-site and two-site kinetic sorption models were used in
the analyses of the MS2 breakthrough curves at wells WP2 and WP3. Optimization results
for both models are given in Table 8 and in Fig. 9. In one set of optimization runs, we
assumed that the transport parameters are constant throughout the subsurface. In thesecond set of optimization runs, we allowed the immobilization, attachment, and detachment coefficients to be different in the direction of each well. The second set of
optimizations we ran only for the one-site kinetic sorption model.
Similarly as for the dune recharge study, the two-site model performs better than the
one-site model in predicting the measured MS2 concentrations, both in terms of the peak

132

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

Fig. 8. Measured and optimized NaCl breakthrough curves at three wells for the deep well injection experiment.
The distance in the legend corresponds with the location of the monitoring well with respect to the injection
and pumping well, with positive number corresponding with monitoring well located between the two majorwells.

concentration and subsequent declining concentrations (Fig. 9). Also, similarly as for thedune recharge study, the sum of the attachment coefficients katt1 and katt2 for the two-site
model is approximately equal to the attachment coefficient katt1 for the one-site model.
Attachment to the type-2 sites is again significantly slower than attachment to the type-1
sites, with detachment from the type-2 sites faster than that from the type-1 sites. This

Table 8Parameter values for MS2 breakthrough curves for the deep well injection experiment
Model

ll
a

One-site (2 wells
simultaneously)
Two-site (2 wells
simultaneously)
One-site (1st well)0.039aOne-site (2nd well)0.039aaParameter fixed.

0.039

0.039a

ls

katt1

kdet1

0.163 F0.08120.0551 F
0.438
0.112 F
0.04870.121 F
0.065

6.276 F
0.597
4.91 F
5.30
5.85 F
0.187
2.49 F
0.060

0.0165 F
0.02060.00170 F
0.01430.00278 F
0.002100.000366 F
0.00050

katt2

kdet2

SSQ
0.672

1.68 F
6.12

0.479 F
2.01

0.599
0.07890.0690J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

133

Fig. 9. Measured and optimized MS2 breakthrough curves at two observation wells for the deep well injection
experiment assuming homogeneous and heterogeneous subsurface. Two-dimensional one-site and two-site kinetic
sorption models were used in the optimization.

leads to a much larger residence time for MS2 on the type-1 sites and more MS2
becoming attached to these sites compared to sites-2. Thus, again, attachment to type-1
sites mainly determines virus removal. Because of the sparseness of MS2 data, theconfidence of optimized parameters is relatively small, especially for the two-site model.
The estimates of ls from the one-site model are very high and likely an overestimation,
whereas that from the two-site model approximates that of ll. Possibly the value of ls isvery similar to that of ll, because the tails of the breakthrough curves may be evenflatter. This is indicated by the three to four observations that make up the later part of
the tails.
Excellent agreement between measured and optimized MS2 concentrations was
obtained when we allowed the immobilization, attachment, and detachment coefficientsfor the one-site kinetic sorption model to vary in the direction of each well (Fig. 9).
Obviously, this can be explained by the fact that WP3 (8 m) lies within the oxygen zone,
whereas WP2 (12 m) lies outside this zone (Fig. 2). Within the oxygen zone, favorablesites of ferric oxyhydroxides are present, but not beyond that zone (Schijven et al., 2000).
The MS2 concentration data were, however, too sparse to fit parameters for the two-sitekinetic sorption model with varying parameters in the two directions. Possibly, some
improvement may be expected when assuming ls is equal ll.

6. Concluding remarks
We implemented a reversible two-site kinetic sorption model into the HYDRUS-1Dand HYDRUS-2D numerical codes and used the resulting models to analyze the transport

134

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

of MS2 and PRD phages in two field studies. Differences between the one- and twodimensional modeling approaches, differences between the one- and two-site kinetic
models, and the role of heterogeneities in the soil properties were investigated.
For the dune recharge experiment, the one-dimensional one-site model did notcompletely describe the tails of the measured breakthrough curves. While the onedimensional one-site model predicted a sudden decrease in the virus concentrations
immediately after their peaks, measured data displayed a relatively smooth transition of
the declining limb to the later straight part of the tail. The one-dimensional two-site model
fitted the course of the breakthrough curves very well. The two-dimensional one-site
model predicted a more gradual decrease in the virus concentrations after their peaks thanthe one-dimensional one-site model, but not as good as the one-dimensional two-site
model. Because the two-dimensional one-site sorption model could describe the measured
breakthrough curves already reasonably well, the breakthrough curves did not provide
enough information for successful estimation of two-site sorption parameters. Transverse
dispersion seemed to partially explain features of the breakthrough curve that were thought
to be the effects of a second type of kinetic sites when the one-dimensional model wasused. Apparently, dimensionality of the problem can partly predict the smooth decrease inconcentration after peak breakthrough. Still, the two-dimensional two-site model did
provide better predictions of the breakthrough curves, partly due to the higher order of
freedom (more fitted parameters) of the inverse solution. Values for katt2 and kdet2 could not
be determined at the larger distances from the source (third and fourth monitoring wells)
using the two-dimensional model, thus suggesting that either a second type of kinetic sites
is present within the first few meters of dune passage, but not beyond the second
monitoring well, or that effects of soil heterogeneity and dimensionality of the problem
overshadowed this process. Variations between single collector efficiencies were relatively
small, whereas collision efficiencies varied greatly. This implies that the nonlinear removalof MS2 and PRD1 is mainly caused by variations in interactions between grain and virus
surfaces rather than by physical heterogeneity of the porous medium.
Similarly, a two-site model performed better than the one-site model in describing MS2
concentration for the deep well injection study. However, the concentration data were too
sparse in this study to have much confidence in the fitted parameters.

Acknowledgements
S.M. Hassanizadeh and M.Th. van Genuchten are greatly acknowledged for their expert
comments.

References
Bales, R.C., Li, S., Maguire, K.M., Yahya, M.T., Gerba, C.P., Harvey, R.W., 1995. Virus and bacteria transport ina sandy aquifer, Cape Cod, MA. Ground Water 33, 653 – 661.
Bales, R.C., Li, S., Yeh, T.C.J., Lenczewski, M.E., Gerba, C.P., 1997. Bacteriophage and microsphere transport insaturated porous media: forced-gradient experiment at Borden, Ontario. Water Resour. Res. 33, 639 – 648.
Bear, J., 1972. Dynamics of Fluid in Porous Media. Elsevier, New York.

J.F. Schijven, J. Šimůnek / Journal of Contaminant Hydrology 55 (2002) 113–135

135

Celia, M.A., Bouloutas, E.T., Zarba, R.L., 1990. A general mass-conservative numerical solution for the unsaturated flow equation. Water Resour. Res. 26, 1483 – 1496.
DeBorde, D.C., Woessner, W.W., Lauerman, B., Ball, P.N., 1998. Virus occurrence in a school septic system andunconfined aquifer. Ground Water 36, 825 – 834.
DeBorde, D.C., Woessner, W.W., Kiley, Q.T., Ball, P.N., 1999. Rapid transport of viruses in a floodplain aquifer.
Water Res. 33, 2229 – 2238.
Hassanizadeh, S.M., Schijven, J.F., 2000. Use of bacteriophages as tracers for the study of removal of viruses. In:
Dassargues, A. (Ed.), Tracers and Modeling in Hydrogeology. Proceedings of TRAM’2000 held in Liege,
Belgium, 23 – 26 May 2000, pp. 167 – 174.
Havelaar, A.H., van Olphen, M., Drost, Y.C., 1993. F-specific RNA bacteriophages are adequate model organisms for enteric viruses in fresh water. Appl. Environ. Microbiol. 59, 2956 – 2962.
Marquardt, D.W., 1963. An algorithm for least-squares estimation of nonlinear parameters. J. Soc. Ind. Appl.
Math. 11, 431 – 441.
Martin, R.E., Bouwer, E.J., Hanna, L.M., 1992. Application of clean bed filtration theory to bacterial depositionin porous media. Environ. Sci. Technol. 26, 1053 – 1058.
Pieper, A.P., Ryan, J.N., Harvey, R.W., Amy, G.L., Illangasekare, T.H., Metge, D.W., 1997. Transport andrecovery of bacteriophage PRD1 in a sand and gravel aquifer: effect of sewage-derived organic matter.
Environ. Sci. Technol. 31, 1163 – 1170.
Ryan, J.N., Elimelech, M., Ard, R.A., Harvey, R.W., Johnson, P.R., 1999. Bacteriophage PRD1 and silica colloidtransport and recovery in an iron oxide-coated sand aquifer. Environ. Sci. Technol. 33, 63 – 73.
Schijven, J.F., 2001. Virus Removal from Groundwater by Soil Passage. Modeling, Field and Laboratory Experiments. Ph.D.-thesis, Delft University of Technology, ISBN 90-646-4046-7.
Schijven, J.F., Hassanizadeh, S.M., 2000. Removal of viruses by soil passage: overview of modeling, processes
and parameters. Crit. Rev. Environ. Sci.Technol. 30, 49 – 127.
Schijven, J.F., Hoogenboezem, W., Hassanizadeh, S.M., Peters, J.H., 1999. Modelling removal of bacteriophages
MS2 and PRD1 by dune infiltration at Castricum, the Netherlands. Water Resourc. Res. 35, 1101 – 1111.
Schijven, J.F., Medema, G., Vogelaar, A.J., Hassanizadeh, S.M., 2000. Removal of microorganisms by deep wellinjection. J. Contam. Hydrol. 44, 301 – 327.
Selim, H.M., Schulin, R., Flühler, H., 1987. Transport and ion exchange of calcium and magnesium in an
aggregated soil. Soil Sci. Soc. Am. J. 51 (4), 876 – 884.
Šimůnek, J., Šejna, M., van Genuchten, M.Th., 1998. The HYDRUS-1D software package for simulating the onedimensional movement of water, heat, and multiple solutes in variably-saturated media. Version 2.0, IGWMCTPS-70. International Ground Water Modeling Center, Colorado School of Mines, Golden, Colorado, 202 pp.
Šimůnek, J., Šejna, M., van Genuchten, M.Th., 1999. The HYDRUS-2D software package for simulating twodimensional movement of water, heat, and multiple solutes in variably saturated media. Version 2.0, IGWMCTPS-53. International Ground Water Modeling Center, Colorado School of Mines, Golden, Colorado, 251 pp.
Stuyfzand, P.J., 1999. Deep well injection in Zuid-Oost Nederland (DIZON). Final report on quality changes
during soil passage. Kiwa, Nieuwegein, The Netherlands, report KOA 99.054 (In Dutch).
Toride, N., Leij, F.J., van Genuchten, M.Th., 1995. The CXTFIT code for estimating transport parameters fromlaboratory or field tracer experiments. Version 2.0. US Salinity Laboratory, Agricultural Research Service, US
Department of Agriculture, Riverside, CA, Report No. 137.
van Genuchten, M.Th., Wagenet, R.J., 1989. Two-site/two-region models for pesticide transport and degradation:
theoretical development and analytical solutions. Soil Sci. Soc. Am. J. 53, 1303 – 1310.
Vogelaar, A.J., van Baar, M.J.C., Peters, J.H., Bergsma J., Report of NaCl-tracer experiments in December 1996
at Castricum, 1997 Kiwa, Nieuwegein, The Netherlands, report SWI 97.115 (in Dutch).
Yao, K.M., Habibian, M.T., O’Melia, C.R., 1971. Water and wastewater filtration: concepts and applications.
Environ. Sci. Technol. 5, 1105 – 1112.
Yates, M.V., 1985. Septic tank density and groundwater contamination. Ground Water 23, 586 – 591.
Yates, M.V., Yates, S.R., Wagner, J., Gerba, C.P., 1987. Modelling virus survival and transport in the subsurface.
J. Contam. Hydrol. 1, 329 – 345.