Journal of Hydrology 578 (2019) 124082

Contents lists available at ScienceDirect

Journal of Hydrology
journal homepage: www.elsevier.com/locate/jhydrol

Research papers

Modeling Virus Transport and Removal during Storage and Recovery in
Heterogeneous Aquifers

T

Saeed Torkzabana,b,⁎, Mark Hockinga, Scott A. Bradfordc, Sara S. Tazehkanda,b,Salini Sasidharanc,d, Jiří Šimůnekd
aGeological Survey of Victoria, PO Box 500, Melbourne, Victoria 3002, Australia
National Centre for Groundwater Research and Training, Flinders University, Adelaide, SA 5001, Australiac
US Salinity Laboratory, USDA, ARS, Riverside, CA, USAd
Department of Environmental Sciences, University of California Riverside, Riverside, CA, USAb

ARTICLE INFOABSTRACTThis manuscript was handled by Peter K.
Kitanidis, Editor-in-Chief, with the assistance ofJian Luo, Associate Editor

A quantitative understanding of virus removal during aquifer storage and recovery (ASR) in physically andgeochemically heterogeneous aquifers is needed to accurately assess human health risks from viral infections. A
two-dimensional axisymmetric numerical model incorporating processes of virus attachment, detachment, andinactivation in aqueous and solid phases was developed to systematically evaluate the virus removal performance of ASR schemes. Physical heterogeneity was considered as either layered or randomly distributed hydraulic conductivities (with selected variance and horizontal correlation length). Geochemical heterogeneity inthe aquifer was accounted for using Colloid Filtration Theory to predict the spatial distribution of attachment
rate coefficient. Simulation results demonstrate that the combined effects of aquifer physical heterogeneity andspatial variability of attachment rate resulted in higher virus concentrations in the recovered water at the ASR
well (i.e. reduced virus removal). While the sticking efficiency of viruses to aquifer sediments was found to
significantly influence virus concentration in the recovered water, the solid phase inactivation under realistic
field conditions combined with the duration of storage phase had a predominant influence on the overall virus
removal. The relative importance of physical heterogeneity increased under physicochemical conditions that
reduced virus removal (e.g. lower value of sticking efficiency or solid phase inactivation rate). This study provides valuable insight on site selection of ASR projects and an approach to optimize ASR operational parameters(e.g. storage time) for virus removal and to minimize costs associated with post-recovery treatment.

Keywords:
Aquifer storage and recovery (ASR)
Virus attachment, detachment, and inactivation
Heterogeneous aquifers

1. Introduction
Managed aquifer recharge (MAR) is increasingly used worldwide to
replenish aquifers to meet growing demands for a secure and sustainable supply of potable and non-potable water. MAR is a broad term for avariety of methods of recharging and recovering water from aquifers.
For example, one of the most common types of MAR is the Aquifer
Storage and Recovery (ASR) technique (Page et al., 2016) which injects
and recovers water via the same well. A variety of water types such as
stormwater runoff and secondary treated wastewater can be used assource water in ASR schemes. The source water (e.g., urban stormwater) is usually contaminated with pathogenic microbes such as
viruses, bacteria, and protozoa. Viruses are generally designated as the
microbial pathogen of greatest risk because of their low infectious doseand potential to travel long distances in aquifers (Reynolds et al., 2008).

⁎

Several studies have reported that viruses can travel several hundred
meters or may be viable for several months (or even years) under certain environmental conditions (Sidhu et al., 2015; Torkzaban et al.,2006; Schijven et al., 2016). The potential for large transport distanceand long-term survival for pathogenic viruses in aquifers is a significant
public health concern. Indeed, one of the main impediments to theuptake of ASR has been attributed to the presence of pathogenic virusesin the recovered water (NRMMC-EPHC-NHMRC, 2009). Recovered ASR
water is therefore subjected to expensive post-treatments to remove/inactivate pathogens (Dillon et al., 2010).
Schijven and Hassanizadeh (2000) presented a comprehensive review on factors and processes influencing virus removal in laboratory
and field studies. Processes influencing virus transport and removal inaquifers include advective and dispersive transport, attachment to anddetachment from sediment surfaces, and inactivation in aqueous and on

Corresponding author at: Geological Survey of Victoria, PO Box 500, Melbourne, Victoria 3002, Australia.
E-mail address: saeed.torkzaban@ecodev.vic.gov.au (S. Torkzaban).

https://doi.org/10.1016/j.jhydrol.2019.124082
Received 19 June 2019; Received in revised form 21 August 2019; Accepted 26 August 2019
Available online 27 August 20190022-1694/ © 2019 Elsevier B.V. All rights reserved.

Journal of Hydrology 578 (2019) 124082

S. Torkzaban, et al.

solid phases (Sasidharan et al., 2018; Anders and Chrysikopoulos,
2005). In addition, previous research suggests that physical and geochemical heterogeneities of the aquifer play an important role in theprocess of microbial transport and removal (Bhattacharjee et al., 2002;
Maxwell et al., 2007; Rehmann et al., 1999; Katzourakis andChrysikopoulos, 2018). For example, preferential flow paths, such as
fractures or regions with high hydraulic conductivity, in a physicallyheterogeneous aquifer are recognized to significantly increase thetransport of contaminants and microbial pathogens in aquifers(Masciopinto et al., 2008; Bradford et al., 2017). Accurate predictionsof virus removal during ASR, therefore, require consideration of removal factors and aquifer heterogeneity. These predictions can then be
used to evaluate the suitability of the recovered water for its intended
use and to determine the required degree of post-recovery treatment.
Aquifers are known to be physically and geochemically heterogeneous (Harvey et al., 1993, Kitanidis, 1997; Dagan et al., 2003). Poor
characterization of aquifer heterogeneity is an important factor that
hampers accurate predictions of both flow and contaminant transport.
In contrast to homogeneous conditions, aquifer heterogeneity is expected to produce irregular concentration fronts with enhanced virustransport in certain parts of the aquifer. The last two decades have seen
the development of a relatively large number of models that predictflow and solute transport in heterogeneous porous media (Jarvis andLarsbo, 2012; Šimůnek et al., 2003; Wang et al., 2014). However, mostmicrobial modeling studies reported in the literature have only dealt
with homogeneous porous media (Anders and Chrysikopoulos, 2005;
Jin and Flury, 2002; Rockhold et al., 2004; Unc and Goss, 2004). Relatively speaking, there are fewer studies that have considered the effectof physical and chemical aquifer heterogeneity on virus transport andremoval (Hornstra et al., 2018; Bhattacharjee et al., 2002; Pang, 2009;
Sinreich and Flynn, 2011). Recently, Katzourakis and Chrysikopoulos(2018) extensively investigated the effect of spatially variable attachment coefficient on biocolloid transport in geochemically heterogeneous porous formations and suggested that neglecting to account foraquifer chemical heterogeneity may lead to erroneous predictions ofvirus transport in porous media. To date, the effect of aquifer heterogeneity on the virus removal performance during ASR, in which injection and recovery occur from the same well, has not yet been studied.
Few ASR field experiments have been conducted to investigate virus
removal because they are expensive, time-consuming, and experimentally challenging (Sidhu et al., 2015). Nonetheless, several laboratory
and field studies using bacterial viruses as surrogates for human viruses
during aquifer recharge have yielded significant insight into the processes of virus transport and removal during ASR. For example, virus
attachment rates appear to be faster than detachment rates in fieldstudies (Bales et al., 1997; Schijven et al., 1999, Schijven and
Hassanizadeh, 2000, Anders and Chrysikopoulos, 2005, Hornstra et al.,2018). Several lab studies have shown that virus inactivation rate ismuch greater when they are attached on mineral surfaces than in freewater (Hsu et al., 2011; Sasidharan et al., 2018). In contrast, several
other studies reported that mineral surfaces such as quartz sand andclay-sized particles offered protection against virus inactivation undercertain conditions (Anders and Chrysikopoulos, 2006; Chrysikopoulos
and Aravantinou, 2012). Consequently, the relative importance of solid
phase inactivation remains a controversial (or perhaps case-specific)
issue. Recently, Sasidharan et al. (2018) reported that virus removal byattachment and solid phase inactivation may be enhanced by adjusting
the storage time in the ASR operation. However, no studies have systematically investigated the influence of storage time on the removal ofviruses due to attachment and solid phase inactivation of viruses during
ASR.
In this paper, we have developed a two-dimensional axisymmetricvirus transport model to assess the removal of viruses in an ASR scheme
under various operational conditions. In this model, aquifer heterogeneity is considered either as layered or random distributions ofhydraulic conductivity with various standard deviations and correlation
lengths. The virus transport equation is coupled with the fluid flow
equation. The model is used to perform a systematic investigation on
the influence of virus removal parameters (attachment, detachment,
and solid and aqueous phase inactivation), physical and chemical
aquifer heterogeneity, and storage time on virus removal during ASR.
We show that the incorporation of aquifer heterogeneity into thetransport model results in predictions which are significantly different
than those obtained from models developed for homogeneous aquifers.
2. Materials and method
2.1. Governing equations
We assumed an ASR system with a fully penetrated pumping wellinstalled in a confined aquifer. The three-dimensional domain wasmodeled using a two-dimensional axisymmetric cross-section. The
vertical axis of rotation was located at the pumping well (herein referred to as the ASR well). An ASR cycle consisted of an injection,
storage, and recovery phase. The groundwater velocity was assumed to
be at steady-state during each phase and the transition period between
two phases was neglected. Furthermore, the effect of flow reversal
during the recovery phase on mobilization of attached viruses wasneglected based on available results in the literature (Sasidharan et al.,2018). During the injection phase, water with a constant virus concentration was injected into the aquifer and a concentration front
moved away from the ASR well. During the recovery phase, stored
water was extracted via the same ASR well and the concentration front
moved toward the well. For simplicity, we assume that the injection andrecovery flow rates were the same and the effect of regional groundwater flow was not considered during these phases. Our axisymmetricmodel does not allow the simulation of non-symmetric conditions suchas regional groundwater flow. Modelling virus transport in an ASR
scheme with a regional hydraulic gradient requires a 3-dimensionalmodel which introduces extra variables and complicates the investigation. In this paper, we have sought to understand the fundamental removal processes of viruses in ASR systems in which theaverage water velocity during injection and recovery phases is significantly greater than the native groundwater velocity. However, tostudy virus removal during storage phase, the water velocity was determined based on a regional hydraulic gradient of 0.001 to estimate
the attachment rate coefficient (to be discussed later).
The governing equation for the two-dimensional, axisymmetric, andtransient water flow in a heterogeneous confined aquifer is as follows(Schijven et al., 2010):

Ss

h
h
1
h
=
kz
+
kr r
t
z
z
r r
r

(1)

where h [L, where L denotes units of length] is the hydraulic head, t [T,where T denotes units of time] is the time, Ss [L−1] is the specificstorage, and kz and kr [L T−1] are the hydraulic conductivity in the z(vertical) and r (radial) [L] directions. The spatially distributed hydraulic heads were determined and then used to calculate the watervelocity by applying Darcy’s law as vr = ( kr / n)( h/ r ) andvz = ( k z /n)( h/ z ) where vr [L T−1] and vz [L T−1] are the vertical andradial pore water velocities, respectively, with porosity denoted by n[–].
Virus transport in the aquifer was described using the advection–dispersion equation along with appropriate terms for virus attachment, detachment, and inactivation:

C
S
+
= Dzz
t
t
Q = µl C + µs S
2

2C

z2

+

1
C
Drr r
r r
r

vz

C
z

vr

C
r

Q
(2)

Journal of Hydrology 578 (2019) 124082

S. Torkzaban, et al.

S= K att CtK det S

size particles. The sticking efficiency parameter (α) is defined as theratio of the number of collisions that result in attachment to the total
number of collisions (Tufenkji and Elimelech, 2004). Essentially, thevalue of α represents the probability that a collision will succeed inattachment. The significance of the sticking efficiency is that, in contrast to Katt, it is considered to be independent of flow velocity.
In this study, the representative sediment grain size was taken to bed10 [L] (Bradford et al., 2017) i.e., where 10% of the sediment mass isfiner than d10. The value of d10 was related to the hydraulic conductivity (K) distribution using the Kozeny–Carman equation (Bear,
1972) as:µs S

where C* is the normalized virus concentration defined as C/Cin (where
Cin [N L−3] is the number [N] of virus per unit volume in water during
the injection phase), S* is the normalized adsorbed virus concentration
S
defined as S = nCb (S [N M−1] is the number of attached viruses perinunit mass [M] of sediment and ρb [M L−3] is the bulk density of thesediment), µl and µs [T−1] are the inactivation rate coefficients for thefree viruses and attached viruses, respectively, and Katt and Kdet [T−1]are rate constants for virus attachment and detachment, respectively,that will be described later. Finally, Dzz and Drr [L2 T−1] are the dispersion coefficients in z and r directions, given by:

Dzz =
Drr =

z |v| + Ddiff + ( rz |v| + Ddiff + ( r

z)

z)

vz vz
|v|

d10 =

(3)

Spatial variations in the hydraulic conductivity result in a heterogeneous flow field that influences virus transport and the concentration
distribution in the aquifer. Two types of physical heterogeneity were
investigated: layered and random heterogeneities. In a layered system,
the aquifer was divided into three layers with separate hydraulic conductivities. Aquifers with fractures, large blocks of macro-pores, orvarious sedimentary deposits may be described as a layered aquifer(Ibaraki and Sudicky, 1995).
Evidence from field-scale hydraulic conductivity measurements indicates that the spatial distribution of hydraulic conductivity is lognormally distributed (Sudicky and MacQuarrie, 1989). Indeed, Freeze(1975) pointed out that the lognormal distribution of K typically has astandard deviation (in log base 10 units) ranging from 0.2 to 2.0. In thisstudy, we employed the HYDRUS (2D/3D) computer software (Simunek
et al., 2016) to generate stochastic (random) distributions of hydraulicconductivity scaling factor. HYDRUS uses the spectral approach of
Mejia and Rodriguez-Iturbe (1974) which is based on an auto-correlated log-normal distribution. Using this method, a specific K value canbe assigned to the flow domain that is then multiplied in each node by ascaling factor (λ). The stochastic distribution requires three input
parameters: the standard deviation (σ) of log10((λ)), and its correlation
length in the lateral (the X-correlation length) and vertical (the Z-correlation length) directions. The value of σ determines the extent ofvariations in the scaling factors, with higher values leading to higher
variations in scaling factors. On the other hand, the correlation length isa measure of the distance in a specific direction that scaling factors are
related. A high value of the X-correlation length means that the scaling
factor maintains similar values for a greater horizontal distance. In thisstudy, various combinations of σ and lateral correlation length wereconsidered.

(4)

where d50 [L] is the median grain size, α [−] is the sticking efficiency
and η [−] is the single-collector contact efficiency representing theratio of the number of particles approaching the collector to the numberof particles striking a collector. A correlation equation developed by
Messina et al. (2015) was used to calculate η as a function of parameterssuch as Peclet number, grain size, virus size, and virus density(1.05 g cm−3). This recently developed equation provides η values
lower than one over a wide range of parameters and it is valid for virusTable 1
Hydrologic and transport parameters used in the model simulations.
Parameter

Basic value

Well radiusAquifer thickness
Anisotropy ration (r:z)
Injection and recovery rates
Hydraulic conductivity, K
Specific storage, Ss
Dispersity (longitudinal), αr
Dispersity (transverse), αz
Porosity, n

0.1 m
20 m
1
300 m3 day−1
5 m day−1
1 × 10−5 m−1
0.1 m
0.01 m
0.4

Transport Parameters
Median grain diameter, d50
Virus diameter, dp
Sticking efficiency, αDetachment rate coefficient, Kdet
Aqueous inactivation rate coefficient, µl
Solid phase inactivation rate coefficient, µs

2.5 × 10−4 m
5.0 × 10−8 m
1 × 10−3
0.001 day−1
0.01 day−1
0.05 day−1−1

2.3. Physical heterogeneity of aquifer

Table 1 shows the baseline flow and transport parameters used in
the simulations. The ranges of these parameters, typically estimated inlaboratory and field studies, are also listed in the table. All simulations
were conducted using parameter values that lie within the specified
ranges. Even in cases of randomly heterogeneous conditions, the hydraulic conductivity and geochemical distributions were determinedsuch that the specified ranges of these heterogeneities were never exceeded.
For our modeling purposes, attachment and detachment weremodeled as first-order processes. According to the Colloid Filtration
Theory (CFT), the attachment coefficient is assumed to be related to theaverage flow velocity as (Yao et al., 1971):

v

(5)
−1

where w [M L ] is the density of water, Xw [M L T ] is the dynamic viscosity of water, and g [L T−2] is the acceleration due togravity. Note that d50 = d10Ui, where Ui is the uniformity coefficientthat was taken to be 4. The use of Eq. (5) in this manner allowed d50 in
CFT calculations to be related to K.

vr vr
|v|

2.2. Model parameters

3 (1 n)
2 d50

0.5

−3

where z and z [L] are the vertical and radial dispersivities,|v| [L T−1]
is the velocity magnitude and Ddiff [L2 T−1] is the molecular diffusion
coefficient.

K att =

KXw 180(1 n) 2
n3
wg

Comments

2.4. Correlation of the sticking efficiency with the hydraulic conductivity

100–1000
0.1–100

A randomly distributed K leads to spatial variations of the porewater velocity and d10 as described by Eqs. (1) and (5). Moreover, theparameter α in Eq. (4) was assumed to be correlated with K. Consequently, spatial variability in K produced random values of Katt. Unlike
the grain size parameter (d50), there are no established empirical expressions relating the value of α to K. The empirical approach of
Garabedian et al. (1988) was therefore employed to relate α to K as:

0–1 × 10−2
0–0.1 day−1
0–0.1 day−1
0–0.1 day−1ln( ) = ln( ) + ln(K )where3is the correlation coefficient and(6)
represents the portion ofJournal of Hydrology 578 (2019) 124082

S. Torkzaban, et al.

the sticking efficiency that depends on solution (e.g., pH and ionicstrength), sediment (e.g., metal oxide content), and virus properties(e.g., surface properties) that cannot be accounted for by variability in
the K field. Therefore, the effects of variability of solution, sediment,
and virus properties were not directly considered and were beyond thescope of this paper. Using this approach, it is possible for the value of α
to be positively correlated to K ( > 0) , negatively correlated to K
( < 0) , or uncorrelated to K ( = 0) . Note that in case of a constantvalue of α the variability of Katt only occurs through spatial variabilityof d50 and v using Eq. (4). In an aquifer with randomly distributed
heterogeneities, we assumed ( < 0) because previous research foundthat the value of α was higher in a low-conductivity than in a highconductivity porous media (Harvey et al., 1993; Morley et al., 1998).
This finding was attributed to the presence of iron-rich coatings, minerals which were most abundant in the finer-grained sediments,
leading to greater colloid sticking efficiencies. The value of
waschosen as −0.1. The value of was chosen to range between 10-5 forhighly unfavorable attachment conditions (e.g. anoxic conditions) to
10-2 for modestly favorable conditions reported for limestone aquifer(Sasidharan et al., 2017). For the simulations performed in this paper,
detachment was assumed to be uncorrelated with K as there is no established empirical expression in the literature relating the value of Kdet
to K.
Several virus transport models have been developed that incorporate constant first-order inactivation of free (μl) and attached (μs)
viruses (Chrysikopoulos and Sim, 1996; Sim and Chrysikopoulos,
1996). Based on the existence two or more subpopulations of viruses ina suspension with different inactivation rate coefficients, Sim andChrysikopoulos (1996) also developed a transport model, incorporating
kinetic reversible adsorption and different time-dependent inactivation
rate coefficients for suspended and attached viruses. This more sophisticated way of modeling virus inactivation during subsurface
transport under saturated conditions may also be useful when considering heterogeneous populations of different viruses that exhibit
multiphasic inactivation. However, in this study, we modelled virus
inactivation as a first-order decay process with different rate coefficients for free and attached viruses. Literature data show that inactivation rate coefficients of viruses in the solution phase range between 0.005 and 0.1 day−1, depending on the temperature and virus
type (Sidhu et al., 2015, Hornstra et al., 2018). In our simulations, thesolution-phase inactivation coefficient was set to 0.01 day−1. For surface inactivation, the values of µs were chosen to range from 0 to
0.1 day−1. Recent research has shown that values of µs for bacteriophages PRD1 and ΦX174 were several orders of magnitude greaterthan µl (Sasidharan et al., 2018). Bacteriophages MS2, PRD1 andΦX174 are commonly used as surrogates to human viruses (Schijven
and Hassanizadeh, 2000). However, recent experimental studies withreal pathogenic enteric viruses (adenovirus) have reported that theattachment behavior of surrogate phages MS2 and ΦX174 is substantially different to that of adenovirus (Bellou et al., 2015; Schijven
et al., 2003). Therefore, results based on bacteriophages may not always be representative to true pathogenic enteric viruses.

Fig. 1. Cross-section of the axisymmetrical modeling domain and boundary
conditions BC0 (no flux of water or virus), BC1 (water and virus enter and leave
the domain), and BC2 (a constant head for water and zero virus concentration).
Corresponding boundary condition equations for water flow and virus transport
are listed in Table 2.
Table 2
Boundary conditions (see Fig. 1 for a schematic of the modeling domain).
Boundary

Description

Water flow

Virus transport

BC0

No flux of water and virus
across boundary
Screen of the ASR well with aconstant flux of water into (−)
and out (+) of the domain

h
or
r

h
=0
z

C
r

h
=
r

Q
Aw kr

C = 1 (during
injection)

BC1

BC2

Right outer boundary, aconstant head for water andzero virus concentration

C
r

h=B

or

C
z

=0

=0

(during

recovery)C =0

head is constant (equal to the aquifer thickness B in these simulations),
and the virus concentration was zero. BC0 designates no flux of waterand viruses across the boundaries. The coupled groundwater flow andvirus transport equations accounting for all the processes of attachment,
detachment, inactivation, as well as physical and geochemical heterogeneities, were implemented and solved using COMSOL Multiphysics
software (COMSOL, Inc., Palo Alto, CA 94301). We made use of theInfinite Element Domain feature of COMSOL. This feature performs acoordinate scaling in the radial direction to the selected domain suchthat boundary conditions on the outside of the infinite layer are effectively applied at a very large distance. Therefore, unwanted effects ofartificial boundary conditions on the region of interest are suppressed.
This allows modeling details in a transport domain, which is actually
very large or infinite. The simulations were performed using a finite
element mesh comprising 76,822 triangular elements. The optimal
mesh resolution was determined by starting from a smaller number of
elements and doubling the number until the concentration profiles fromtwo consecutive meshes were found to be within a few percent (1–3%).
2.6. Numerical experiments
2.6.1. Virus transport in heterogeneous aquifer
Numerical experiments were initially conducted to highlight thecapabilities of the model and provide some insight on the processes ofvirus transport during ASR in aquifers with different types of heterogeneity. In these simulations, aquifer heterogeneity was consideredeither as layered or random distributions of hydraulic conductivity withvarious standard deviations and correlation lengths. The basic values ofthe model parameters and their range of variation for the numericalsimulations are listed in Table 1. The range of parameter values covered
possible scenarios for virus transport in sandy aquifers, which has been
determined in many field or laboratory studies (Schijven and
Hassanizadeh, 2000).

2.5. Numerical model approach and model setup
The axisymmetric modeling domain is shown in Fig. 1. The initialconditions were defined as a constant hydraulic head and zero virus
concentration everywhere in the domain. The boundary conditions forthe flow and virus transport equations are listed in Table 2. At boundary
BC1, the screen of the ASR well is situated where there is a constant fluxof water into (during the injection phase) and out (during the recovery
phase) of the domain. The value of dispersivity for nodes located in thewell was set to zero allowing to implement Dirichlet boundary conditions (C*=1) for virus injection at BC1 during the injection phase.
During the recovery phase, a zero dispersive flux was assigned forviruses at BC1. At the boundary BC2, it was assumed that the hydraulic4

Journal of Hydrology 578 (2019) 124082

S. Torkzaban, et al.

2.6.2. Virus concentration in the recovered water
Additional simulations were conducted to examine the effect ofremoval parameters (attachment, detachment, and solid and aqueous
phase inactivation), physical and chemical aquifer heterogeneity, andstorage time on the extent of virus removal during ASR. In this case,normalized virus concentrations (C*) in the recovered water at the ASR
well were calculated for a range of aquifer heterogeneities (σ for ln(K),radial correlation length, and correlations of α with K), sticking efficiencies, detachment rates, solid phase inactivation rates, and storage
durations that would be expected in the field. Note that the values of C*
were flux-averaged concentrations during the recovery phase. They
were obtained after an injection phase of 60 days and plotted on alogarithmic scale as a function of recovery time. A constant mean K
value of 5 m day−1 was assumed in all subsequent simulations.
Moreover, all simulations were performed for a log-normally distributed K field with a given parameter set of σ, X, and Z. Unlessotherwise mentioned the storage phase was neglected, and the recovery
phase started immediately after the completion of 60 days of injection.

A
K=50, =1×10-5

C*

K=5, =1×10-3
K=25, =1×10-4

B
K=50, =1×10-3
K=5, =1×10-5
K=25, =1×10-4

Fig. 3. Relative importance of combined effects of layered physical and geochemical heterogeneity on the transport of viruses in a model aquifer. The
normalized virus concentration profiles are shown at t = 30 days for an aquifer
containing three layers with different properties. The values of K and α for eachlayer are shown in each figure. The values of other parameters used in thesimulations are listed in Table 1. The contour labels indicate the normalized
virus concentration C*.

3. Results and discussion
3.1. Virus transport in an aquifer with distinct layered heterogeneity
Simulations in this section considered three distinct formation
layers with different hydraulic conductivities. The value of K equaled
50, 5, and 25 m day−1 in the top, middle, and bottom layers, respectively. Fig. 2 presents distribution of normalized virus concentrations in
the model domain after 30 days of injection. A constant value of α = 0(completely unfavorable for attachment) and α = 1 × 10-4 was used in
the top and bottom figures, respectively. The K value of the top layer
was an order of magnitude larger than the K value in the layer below.
Therefore, water in the top layer flowed faster than the other layers,
and therefore most of the viruses migrated with the water flow throughthis more permeable layer. This example points out the paramountimportance of preferential flow paths on virus transport. It is also notedthat increasing α from 0 to 1 × 10-4 considerably lowered virus
movement in the aquifer.

Fig. 3 shows the combined effects of layered geochemical andphysical heterogeneity on distribution of normalized virus concentrations in the aquifer after 30 days. The aquifer shown here consisted ofthree layers with the same value of K as in Fig. 2, but with different
combinations of α. In Fig. 3A, an inverse relationship between K and α
was assumed; e.g., values of α were assigned to be 1 × 10-5, 1 × 10-3,
and 1 × 10-4 in the top, middle, and bottom layers, respectively. Conversely, a positive correlation between K and α was considered inFig. 3B by assigning α = 1 × 10-3, 1 × 10-5, and 1 × 10-4 in the top,
central, and bottom layer, respectively. All other conditions were kept
the same in the two simulations. The results clearly show that a decreased attachment of viruses in the layer with the lowest value of α
(1 × 10-5) and highest value of K resulted in preferential migration ofviruses through that layer (Fig. 3A). We observe that the combination ofphysical and geochemical heterogeneity significantly affected the
movement patterns of viruses in the aquifer and that an inverse relationship between K and α was a worst-case (maximum) scenario forvirus transport.

AC*

3.2. Virus transport in an aquifer with randomly distributed heterogeneity
The above simulations considered highly idealized scenarios,
whereas natural aquifers can be highly heterogeneous with a lognormal distribution of K. In Fig. 4, we show the influence of a lognormal distributed K (with a mean value of 5 m day−1) on virustransport in the aquifer. The simulations were performed for a valueof 1 × 10-4 and the localized value of α was correlated to the value of K
through Eq. [6] with = −0.1. Under such conditions, variations of K
led to local variations of the pore water velocity, d50, and α in theaquifer. Therefore, spatial variability of these parameters resulted inlocal variations of Katt (Eqs. [4–6]). Fig. 4 shows the simulated distribution of normalized virus concentrations in the aquifer after 30 daysof injection. In this case, the value of σ for ln(K) (e.g., the physicalheterogeneity) increased from 0.7, 1.4, and 2.1 in the top, middle, andbottom figures, respectively. The correlation lengths of X and Z weretaken as 0.5 in these simulations. The virus concentration variability
and location of the front tended to increase with an increase in σ.
Maximum and minimum values of K in the model domain increased anddecreased, respectively, with increasing σ. For example, values of K

B

Fig. 2. Virus concentration profiles during the injection phase at t = 30 days fordifferent layered hydraulic conductivity field as shown; (A) a uniform field of
α = 0 (zero attachment); (B) a uniform field of α = 1 × 10-4. The values ofother parameters are shown in Table 1. The contour labels indicate the normalized virus concentration C*.
5

Journal of Hydrology 578 (2019) 124082

S. Torkzaban, et al.

hydraulic conductivity, acting as preferential pathways, were moreconnected in which the virus advection velocity was significantly
higher. In addition, the rate of virus attachment through the high K
zones decreases with increasing flow velocity, resulting in less attenuated virus concentration profile. The low velocity regions with depleted virus transport are also more prevalent for larger value of X.
Consequently, greater non-uniformity of the virus concentration front
occurred for larger radial correlation length which yielded greater deviations from the homogeneous aquifer with uniform fronts. The effectof heterogeneity was also evident in the distribution of conservative
tracer concentrations which showed more dispersion than the homogeneous case (data not shown).
3.3. Virus concentration in the recovered water
3.3.1. Sensitivity to aquifer heterogeneity
To explore the impact of spatial variability of K on virus concentration in the recovered water, a series of simulations was performed
with different values of σ of ln(K). Values of σ = 0.7, 1.44, and 2.1 werechosen to represent low, intermediate, and high levels of heterogeneityin K, respectively. In these simulations, detachment and inactivation
coefficients were set to zero and thus virus removal occurred only due
to irreversible virus attachment to aquifer sediments. In addition, thewas set to 1 × 10-3, which was a relatively conservativevalue ofvalue based on many laboratory and field studies (Sasidharan et al.,2018; Sasidharan et al., 2017). Note that the localized values of α andd50 were correlated to the value of K through equations (5) and (6). A
similar simulation was also performed for the homogeneous aquifer(σ = 0) with a hydraulic conductivity of 5 m day−1 for comparison
purposes.
Fig. 6 shows values of C* as a function of recovery time for aquifers
with various degrees of heterogeneity. Note that at the start of the recovery phase, the virus concentration was the same as the input concentration during the injection phase (i.e. C*=1) because the inactivation of viruses was not considered. As expected, virus
concentration in the recovered water decreased over the recovery time.
For example, a more than 2 log10 reduction in virus concentration

Fig. 4. Influence of randomly distributed hydraulic conductivity on virustransport. The simulations were performed for a continuous injection of viruses
at the ASR well for 30 days. The value of α was set to 1 × 10−4 and α wascorrelated with K through Eq. (6). The three plots represent the virus concentration distribution (C*) after 30 days. In all cases, the K field was generatedusing a mean value of K = 5 m day−1 and correlation lengths of 0.5. The
standard deviations (σ) of ln(K) were (A) 0.7, (B) 1.4, and (C) 2.1. The values ofother parameters used in the simulations are shown in Table 1.

exhibited two orders of magnitude variations for σ = 2.1. Locations
with high K provided preferential pathways for virus transport, whereas
low K regions behaved as natural barriers to flow that decreased thetransport of viruses. Virus concentrations were therefore very low inregions with low K (see lighter shades in Fig. 4C).
Fig. 5 presents the virus concentration distribution when σ of ln(K)was 2.1, the vertical correlation length (Z) was 0.5 m, and the radial
correlation length (X) was 5 m (Fig. 5A) and 50 m (Fig. 5B). The results
illustrate that the virus concentration fronts moved faster in regionswith higher value of X. When the value of X = 50 m, zones with higher

Fig. 6. Normalized virus concentrations at the ASR well as a function of recovery time obtained from heterogeneous and homogeneous simulations following 60 days of injection. Values of σ of ln(K) were 0.7, 1.44, and 2.1 torepresent low, intermediate, and high levels of heterogeneity in K field withmean K of 5 m day−1. Detachment and inactivation coefficients were set to zeroand thus virus removal occurred only due to virus attachment to aquifer sediments. The value of was set to 1 × 10−3 and the localized values of α and d50
were correlated to K according to equations (5) and (6). A similar simulation
was also performed for the homogeneous aquifer (σ = 0) with a hydraulicconductivity of 5 m day−1 for comparison purposes. The values of other parameters used in the simulations are shown in Table 1.

Fig. 5. Influence of randomly distributed hydraulic conductivity on virustransport. The value of was set to 1 × 10−4 and α was correlated with K using
Eq. (6). The two plots represent the virus concentration distribution (C*) after
30 days of injection. In all cases, the K field was generated using a mean valueof K = 5 m day−1. The standard deviation (σ) of ln(K) was 2.1 with (A) X = 5 m
and Z = 0.5 m and (B) X = 50 m and Z = 0.5 m. The values of other parameters
used in the simulations are listed in Table 1.
6

Journal of Hydrology 578 (2019) 124082

S. Torkzaban, et al.

occurred after approximately 10 days of water recovery. Differences in
C* for the heterogeneous and homogeneous aquifers were initiallysmall because the recovered viruses had only traveled short distances in
the aquifer and experienced little attachment. However, differences in
C* for the heterogeneous and homogeneous aquifers increased at laterstages of recovery, especially for higher values of σ. For example, thevalues of C* for the heterogeneous aquifer with σ = 2.1 was around 4
orders of magnitude higher than those for the homogeneous aquifer by
the end of the recovery phase (day 60). This reflects a decrease in virus
removal with an increase in aquifer physical heterogeneity. For example, a 3 log10 (0.001) reduction in the virus concentration occurred
after 21, 13, 8, and 7 days when σ of ln(K) equaled 2.1, 1.4, 0.7, and 0,
respectively.
The difference between virus removal for heterogeneous andhomogeneous cases can be attributed to spatially variable v, α, and d50
which also results in heterogeneity in Katt. In higher permeable zones adecrease in virus attachment occurs due to lower Katt. Conversely, an
increased virus removal is expected in lower permeable zones due tohigher Katt. The slight difference between the results from homogeneousand heterogeneous aquifers with σ = 0.7 is attributed to the fact thatthe volume of high permeability zones was small and they were mainly
isolated. Therefore, viruses were likely to encounter low permeablezones during their transport in the aquifer which produced higher virus
attachment. Consequently, spatial variability of K in an aquifer withlow amounts of heterogeneity may not have a large influence on virus
removal during ASR. Conversely, more volume of the aquifer wascharacterized by the tails of the K distribution when σ increased to 1.4
or 2.1. The “high K” tail may be thought of as highly conductive regionswith small values Katt. Note that the relationship between Katt and K was
highly nonlinear, and resulted in significant decreases in Katt with increasing K. Therefore, when the volume of highly conductive zones wasincreased with increasing σ, these zones acted as preferential pathwaysfor virus transport with less chance of encountering low permeability
zones. The model simulations suggest that reduced virus removal is
expected in aquifers with high σ compared to homogeneous aquifers.
An additional suite of simulations was conducted to explore theimpact of the radial correlation length (X) on virus removal during the
ASR operation. Three different values for the radial correlation length(X = 0.5, 5, and 50) were examined for a fixed value of σ = 1.4. Recall
that the radial correlation length is a measure of the radial distance thatsimilar values of K are maintained. Fig. 7 shows the values of C* in therecovered water when X = 0.5, 5, and 50 m. Less virus removal

occurred as the value of X was increased. For example, the virus concentration after 30 days of recovery only decreased by 2.5 log10 when
X = 50 m, whereas it was reduced by 5 log10 when X = 0.5 (Fig. 7).
These results indicate that virus removal might be significantly lower inaquifers when high permeability zones with low virus attachment rates(higher d50 and lower α) are laterally connected to greater lengths(higher X). In contrast, an aquifer with a smaller value of X , comprising
many isolated high permeability zones, provides significantly higher
levels of virus removal. The significant difference between the homogeneous and heterogeneous aquifers when σ and X were high, highlights
the importance of physicochemical heterogeneity on the virus removal
performance during ASR.
3.3.2. Influence of sticking efficiency
For the previous simulations, we assumed that the mean sticking
efficiency ( ) was 1 × 10−3. However, this model parameter can varyby orders of magnitude (Hornstra et al., 2018; Bradford and Torkzaban,
2013; Bradford and Torkzaban, 2015; Chu et al., 2003). The sticking
efficiency is a function of the water chemistry as well as virus and sediment surface properties (Knappett et al., 2008; Pham et al., 2009;
Sadeghi et al., 2013). An additional series of simulations was therefore
carried out to explore the impact of the sticking efficiency on virus
removal during ASR in heterogeneous aquifers. In these simulations,the value of varied across the heterogeneous aquifer because it wascorrelated with K using Eq. [6], and the values of detachment and inactivation rate coefficients were set to zero. Fig. 8 presents plots of thepredicted influence of on recovered C* in an aquifer with σ = 1.4 anda horizontal correlation length of 50 m. This value of σ = 1.4 waschosen to represent an intermediate level of heterogeneity in K. Other
model parameters are given in Table 1. It is noted that increasingfrom 1 × 10−5 to 1 × 10−2 significantly increased the removal ofviruses. For example, whenwas 1 × 10−2, the values of C* were
found to reduce by almost 5 log10 in less than 5 days of water recovery.
Conversely, decreasing to 1 × 10−5 resulted in little virus removal;e.g., C* was reduced by less than 0.5 log10 after 60 days. In this case, thebehavior of virus transport in the aquifer was somewhat similar to thatof a conservative tracer. These simulation results indicate that characterization of for a given aquifer appears to be an important consideration for accurately quantifying field-scale virus removal during
the ASR operation.
3.3.3. Influence of virus detachment
Consistent with the colloid filtration theory (Eq. (4)) the previous
simulations assumed irreversible virus attachment onto the sedimentFig. 7. Normalized virus concentrations at the ASR well as a function of recovery time obtained from a heterogeneous aquifer when σ of ln(K) was 1.44
and the correlation length in radial direction (X) was 0.5, 5, and 50 m. Results
obtained after 60 days of injection and for the mean K value of 5 m day−1.
Detachment and inactivation coefficients were set to zero and the value of
was set to 1 × 10−3. The values of other parameters used in the simulations areshown in Table 1.

Fig. 8. Influence of sticking efficiency (α) on virus removal during ASR in aheterogeneous aquifer. Normalized virus concentrations were obtained at the
ASR well during the recovery phase after 60 days of injection. The standard
deviation (σ) of ln(K) was 1.4 with X = 50 m and Z = 0.5 m. The localized valueof α was correlated with K using Eq. (6). The values of other parameters used in
the simulations are shown in Table 1.
7

Journal of Hydrology 578 (2019) 124082

S. Torkzaban, et al.

Fig. 9. Influence of virus detachment rate constants (Kdet) on virus removal
during ASR in a heterogeneous aquifer. Normalized virus concentrations wereobtained at the ASR well during the recovery phase after 60 days of injection.
The standard deviation (σ) of ln(K) was 1.4 with X = 5 m and Z = 0.5 m. The
value of was set to 1 × 10-3 and the localized α was correlated with K using
Eq. [6]. The values of other parameters used in the simulations are shown in
Table 1.

Fig. 10. Influence of inactivation rate constants (µs) of attached viruses on virus
removal during ASR in a heterogeneous aquifer. Normalized virus concentrations were obtained at the ASR well during the recovery phase after 60 days of
injection. The transport parameters were taken as follows: σ = 2.1; X = 5 m;
Z = 0.5 m; =0.001; Kdet = 0.001 day−1; µl = 0.01 day−1. The values of otherparameters used in the simulations are shown in Table 1.

of nucleic acid (Gerba, 1983). The most important factors affecting thevirus inactivation rate include temperature, groundwater microbialactivity, pH, salt species and concentration, and virus type (McCarthy
and McKay, 2004; Schijven and Hassanizadeh, 2000). The inactivation
rate coefficient for human enteric viruses in groundwater, determinedat several MAR sites in Australia, was found to be very small(~0.01 day−1) (Sidhu et al., 2015). However, attachment of viruses tomineral surfaces appears to play an important role in virus inactivation.
Many researchers have reported that virus inactivation is generally
accelerated when they are attached to solid surfaces compared to thosein the liquid phase (Ryan et al., 2002; Sasidharan et al., 2018; Schijven
et al., 1999). This enhanced inactivation has been attributed to strong
adhesive forces between viruses and mineral surfaces (Harvey and
Ryan, 2004). In addition, some research has indicated that the strengthof virus attachment may also increase with the residence time (Mondon
et al., 2003; Torkzaban et al., 2013; Vadillo-Rodriguez et al., 2004; Xu
and Logan, 2006), and this may further enhance solid phase inactivation.

surfaces. This assumption has been confirmed in studies involving oppositely charged viruses and collector surfaces (Foppen et al., 2006;
Torkzaban et al., 2006). On the other hand, transport studies involving
similarly charged viruses and collector surfaces have shown a finite rateof virus detachment (Sasidharan et al., 2017; Schijven et al., 2003).
Hence, a series of simulations were performed to assess the importanceof virus detachment on the removal of viruses in a heterogeneous
aquifer with σ = 1.44 and X = 5 m. Fig. 9 presents illustrative examplesof the predicted influence of Kdet (0, 0.001, 0.005, 0.01 day−1) on C* in
was set to
the recovered water. In these simulations, the value of1 × 10−3 and no virus inactivation was considered. Increasing thevalue of Kdet from zero to 0.01 day−1 increased the concentration ofviruses in the recovered water, especially at later stages of recovery. For
example, the concentration of viruses after 30 days of recovery was
reduced by 3 log10 when the Kdet was zero compared to less than 1 log10
reduction when Kdet was 0.01 day−1. It should be mentioned that themagnitude of detachment rate constants of viruses from natural surfaces can vary widely and depends on the relative strengths of the adhesive and the diffusive forces (Schijven and Hassanizadeh, 2000;
Torkzaban et al., 2015). However, the rate coefficient for virus attachment has typically been reported to be higher than that for virus
detachment in field (Bales et al., 1997; Schijven et al., 1999) and laboratory (Sasidharan et al., 2018; Shen et al., 2012) studies.

3.3.5. Effects of storage phase
Up to this point, the storage phase of the ASR operation has been
neglected to focus the discussion on virus removal processes in heterogeneous aquifers during the injection and recovery phases.
However, the storage phase is an important component of an ASR
scheme during which viruses are transported in the aquifer with thenative groundwater flow. In this section, the effect of the duration ofstorage phase on virus removal is investigated. The groundwater velocity during the storage phase was determined based on a regional hydraulic gradient of 0.001. We also assumed that Katt during the storage
phase was described with Eq. (4). In these simulations, the value of µs
was set to 0.05 day−1 and values of other parameters were the same as
those for the simulations shown in Fig. 10. For the conditions investigated in Fig. 11, increasing the duration of storage phase resulted
in a greater amount of virus removal. The marked effect of storage
phase on the extent of virus removal is attributed to the higher value ofµs employed in these simulations. As expected, larger values of µs resultin a greater decrease of the concentration of infective attached viruses
during the storage phase. Therefore, the rate of virus detachment isdecreased and less viruses can be detected in the recovered water.

3.3.4. Influence of solid phase inactivation and storage time
To explore the impact of solid phase inactivation (µs) on virus removal during ASR, a series of simulations was conducted with different
values of µs. Fig. 10 shows values of C* during the recovery phase froma heterogeneous aquifer when µs ranged from 0 to 0.1 day−1. The othertransport parameters were taken as follows: σ = 2.1; X = 5 m;
Z = 0.5 m; =0.001; Kdet = 0.001 day−1; and µl = 0.01 day−1. Notethat a heterogeneous domain with a very large σ provides a worst-casescenario for determining virus transport and risk assessment. As expected, increasing µs produced greater amounts of virus reduction in therecovered water due to enhanced inactivation of attached viruses on thesediment surfaces. In the initial periods of recovery phase, the effect ofµs on enhancing virus removal was small because of the short residence
time for virus attachment and inactivation. Conversely, µs had a marked
effect on reducing C* in later stages of recovery. Virus transport models
which does not incorporate the effect of solid phase inactivation of
attached viruses may, therefore, result in erroneous results for the extent of virus removal in ASR schemes.
Inactivation occurs when viruses lose their ability to infect host cells
and replicate because of disruption of proteins and/or the degradation

3.4. Virus removal for an ASR case study
An additional suite of simulations was conducted to evaluate virus
removal in an ASR scheme by assuming the transport parameters as
8

Journal of Hydrology 578 (2019) 124082

S. Torkzaban, et al.

Fig. 11. Normalized virus concentrations at the ASR well as a function of recovery time obtained from a heterogeneous aquifer for four different durationsof storage phase. Results obtained after 60 days of injection and the mean K
value of 5 m day−1. The model parameters were taken as follows: σ = 2.1;
X = 5 m; Z = 0.5 m;= 0.001; Kdet = 0.001 day−1; µl = 0.01 day−1,µs = 0.05 day−1. The groundwater velocity was determined based on a regional
hydraulic gradient of 0.001. The values of other parameters used in the simulations are shown in Table 1.

Fig. 12. Normalized virus concentrations at the ASR well as a function of recovery time obtained from a heterogeneous and homogeneous aquifer. Results
obtained after 60 days of injection and the mean K value was 5 m day−1. The
model parameters were taken as follows:= 0.145; Kdet = 0.23 day−1;µs = 2.0 day−1; µl = 0.07 day−1. These values are reported in Sasidharan et al.
(2017) for three different bacteriophages obtained from column experiments
packed with sediment and stormwater collected from a MAR site in Parafield,
Australia. The stochastic parameters for the heterogeneous aquifer were taken
as σ = 2.07; X = 50 m in order to represent a highly heterogeneous aquifer.
Moreover, the storage phase was ignored, and water recovery started immediately after 60 days of injection. The values of other parameters used in thesimulations are shown in Table 1.

follows:= 0.145;
Kdet = 0.23 day−1;µs = 2.0 day−1;µl = 0.07 day−1. These are averaged values obtained for three different
bacteriophages (MS2, PRD1, and X174) for a limestone aquifer reported in (Sasidharan et al. (2017). In that study, the transport ofbacteriophages was studied in column experiments packed with sediment and stormwater collected from an ASR site in Parafield, SA,
Australia. Virus breakthrough concentrations (BTCs) were successfully
simulated using an advective–dispersive model that accounted for rates
of attachment (Katt), detachment (Kdet), and solid phase inactivation(µs). We employed the reported fitted values of Kdet, µs, measured µl,
and calculated in our model to evaluate the removal of viruses in bothhomogeneous and heterogeneous aquifers. The stochastic parameters
for the heterogeneous aquifer were taken as σ = 2.07, Z = 5 m, andX = 50 m to represent a highly heterogeneous aquifer with K values
exhibiting three orders of magnitude variations and many connectedhigh permeable lenses in the radial direction. Other model parameters
are given in Table 1. Moreover, the storage phase was ignored, and therecovery phase was immediately started after 60 days of injection.
The combined effects of , Kdet, µs, and µl yielded interesting results(Fig. 12). Under the examined conditions, model results demonstratethat virus removal was significantly higher than those obtained in ourprevious simulations because of greater rate constants of virus attachment and inactivation reported for the aquifer sediment. Consequently,
the vast majority of viruses were attached and subsequently inactivatedon the solid phase within just a few meters from the ASR well. For
example, virus concentration decreased by 6 log10 after about 5 days ofrecovery. This time corresponds to a travel distance of 7.6 m from thewell. The results show that differences between the magnitudes of virus
removal for heterogeneous and homogeneous cases were negligible.
Because values of α and µs reported for this case study were significantly
high, the spatial distribution of water velocities only had a slight influence on the recovered concentrations of viruses. These findings
imply that if stormwater needs a 6 log10 (1 million-fold) reduction invirus concentration before it can serve as drinking water, the practice of
ASR will lead to 92% of the recovered water will meet the requirement.
The remaining 8% of the initial recovered water that exceeds thisthreshold should be further treated or discarded.

and removal during ASR. Both layered and log-normally distributed
heterogeneity can result in substantial changes in virus transport patterns. Our modeling results indicate that incorporation of aquifer heterogeneity resulted in faster virus transport and smaller virus removal
during ASR compared to those of a homogeneous aquifer. In order topredict the virus treatment capacity of the aquifer in an ASR scheme,
some information regarding the degree of aquifer heterogeneity would,
therefore, be required.
The model also provides a valuable insight on the impact of various
removal processes such as attachment, detachment, and solid phasevirus inactivation. The virus sticking efficiency and surface inactivation
terms were found to have a great influence on virus concentrations
during the recovery phase. Large detachment constants may result in an
increased virus concentration in the recovered water during the recovery phase unless the surface inactivation rate is very large.
Furthermore, the model simulations showed that the solid phase inactivation had a greater influence on the overall virus removal compared to the liquid phase inactivation under realistic field conditions.
Results of this study suggest that simpler models that account forvirus transport through a homogeneous aquifer and ignore the influence of aquifer heterogeneity and surface inactivation may yield amisleading picture of virus removal during ASR. It is important to accurately quantify the virus sticking efficiency, detachment, and inactivation of viruses by at least conducting a series of laboratory
column experiments employing aquifer sediments. One way to potentially minimize the risk of viruses during ASR operations is to provide
an adequate residence time for solid phase inactivation. This model,
coupled with reliable estimates of virus removal parameters, may leadto adequate predictions of natural treatment performance of aquifers
targeted for ASR schemes.
Quantitative microbial risk assessment calculations for MAR systems have considered that liquid phase virus inactivation was the only
reliable mechanism for virus removal in the aquifer, and neglected
processes of virus attachment, detachment, and solid phase inactivation(Dillon et al., 2008; Page et al., 2016). Results from this study provide
valuable insight on the relative importance of these removal processes

4. Concluding remarks
The results of this study demonstrate that incorporation of aquifer
heterogeneity has a marked effect on the prediction of virus transport
9

Journal of Hydrology 578 (2019) 124082

S. Torkzaban, et al.

for viruses during ASR and indicate that microbial risk assessments that
only consider liquid phase inactivation may be overly conservative inmany instances. Recognizing the removal of viruses via irreversible
attachment and/or solid phase inactivation during aquifer storage
would help to eliminate some of the expensive post-treatment for recovered MAR water. Stormwater harvesting coupled to ASR provides agreat opportunity to reuse this valuable resource that is often wasted.
Aquifers can act as a treatment barrier within a multiple-barrier approach to harvest and reuse urban stormwater.

Hsu, B.B., Wong, S.Y., Hammond, P.T., Chen, J., Klibanov, A.M., 2011. Mechanism ofinactivation of influenza viruses by immobilized hydrophobic polycations. Proc. Natl.
Acad. Sci. U.S.A. 108 (1), 61–66.
Ibaraki, M., Sudicky, E., 1995. Colloid-facilitated contaminant transport in discretely
fractured porous media: 1. Numerical formulation and sensitivity analysis. Water
Resour. Res. 31 (12), 2945–2960.
Jarvis, N., Larsbo, M., 2012. MACRO (v5. 2): Model use, calibration, and validation.
Transactions of the ASABE, 55(4): 1413–1423.
Jin, Y., Flury, M., 2002. Fate and transport of viruses in porous media. Adv. Agronomy 77,
39–102.
Katzourakis, V.E., Chrysikopoulos, C.V., 2018. Impact of spatially variable collision efficiency on the transport of biocolloids in geochemically heterogeneous porous
media. Water Resour. Res. 54 (6), 3841–3862.
Kitanidis, P.K., 1997. Groundwater Flow in Heterogeneous Formations in Subsurface
Flow and Transport: The Stochastic Approach, edited by G. Dagan, Ê S.P. Neuman,pp. 83–91, Cambridge Univ. Press, 1997.
Knappett, P.S., Emelko, M.B., Zhuang, J., McKay, L.D., 2008. Transport and retention of abacteriophage and microspheres in saturated, angular porous media: effects of ionicstrength and grain size. Water Res. 42 (16), 4368–4378. https://doi.org/10.1016/j.
watres.2008.07.041.
Masciopinto, C., La Mantia, R., Chrysikopoulos, C.V., 2008. Fate and transport of pathogens in a fractured aquifer in the Salento area, Italy. Water Resour. Res. 44 (1),
W01404. https://doi.org/10.1029/2006WR005643.
Maxwell, R.M., Welty, C., Harvey, R.W., 2007. Revisiting the Cape Cod bacteria injection
experiment using a stochastic modeling approach. Environ. Sci. Technol. 41 (15),5548–5558.
McCarthy, J.F., McKay, L.D., 2004. Colloid transport in the subsurface: past, present, andfuture challenges. Vadose Zone J. 3 (2), 326–337.
Mejia, J.M., Rodriguez-Iturbe, I., 1974. On the synthesis of random field sampling fromthe spectrum: an application to the generation of hydrologic spatial processes. Water
Kesour. Hes. 10 (4), 705–711.
Messina, F., Marchisio, D.L., Sethi, R., 2015. An extended and total flux normalized
correlation equation for predicting single-collector efficiency. J. Colloid Interface Sci.
446, 185–193. https://doi.org/10.1016/j.jcis.2015.01.024.
Mondon, M., Berger, S., Ziegler, C., 2003. Scanning-force techniques to monitor timedependent changes in topography and adhesion force of proteins on surfaces. Anal.
Bioanal. Chem. 375 (7), 849–855. https://doi.org/10.1007/s00216-003-1751-2.
Morley, L.M., Hornberger, G.M., Mills, A.L., Herman, J.S., 1998. Effects of transverse
mixing on transport of bacteria through heterogeneous porous media. Water Resour.
Res. 34 (8), 1901–1908.
NRMMC-EPHC-NHMRC, 2009. National Water Quality Management Strategy, Australianguidelines for water recycling: Managing health and environmental risks (phase 2)
Managed Aquifer Recharge. In: Council, N.R.M.M., Council, E.P.a.H., Council, N.H.a.
M.R. (Eds.), National Water Quality Management Strategy, Canberra.
Page, D., et al., 2016. Assessment of treatment options of recycling urban stormwater
recycling via aquifers to produce drinking water quality. Urban Water J. 1–6. https://doi.org/10.1080/1573062X.2015.1024691.
Pang, L., 2009. Microbial removal rates in subsurface media estimated from publishedstudies of field experiments and large intact soil cores. J. Environ. Qual. 38 (4),1531–1559. https://doi.org/10.2134/jeq2008.0379.
Pham, M., Mintz, E.A., Nguyen, T.H., 2009. Deposition kinetics of bacteriophage MS2 tonatural organic matter: role of divalent cations. J. Colloid Interface Sci. 338 (1), 1–9.
https://doi.org/10.1016/J.Jcis. 2009.06.025.
Rehmann, L.L.C., Welty, C., Harvey, R.W., 1999. Stochastic analysis of virus transport inaquifers. Water Resour. Res. 35 (7), 1987–2006.
Reynolds, K.A., Mena, K.D., Gerba, C.P., 2008. Risk of waterborne illness via drinking
water in the United States, Reviews of environmental contamination and toxicology.
Springer, pp. 117–158.
Rockhold, M.L., Yarwood, R., Selker, J.S., 2004. Coupled microbial and transport processes in soils. Vadose Zone J. 3 (2), 368–383.
Ryan, J.N., et al., 2002. Field and laboratory investigations of inactivation of viruses(PRD1 and MS2) attached to iron oxide-coated quartz sand. Environ. Sci. Technol. 36(11), 2403–2413.
Sadeghi, G., Schijven, J.F., Behrends, T., Hassanizadeh, S.M., Van Genuchten, M.T., 2013.
Bacteriophage PRD1 batch experiments to study attachment, detachment and inactivation processes. J. Contaminant Hydrol. 152, 12–17.
Sasidharan, S., Bradford, S.A., Šimůnek, J., Torkzaban, S., Vanderzalm, J., 2017.
Transport and fate of viruses in sediment and stormwater from a Managed Aquifer
Recharge site. J. Hydrol. 555, 724–735.
Sasidharan, S., Bradford, S.A., Šimůnek, J., Torkzaban, S., 2018. Minimizing virustransport in porous media by optimizing solid phase inactivation. J. Environ. Qual.
47, 1058–1067. https://doi.org/10.2134/jeq2018.01.0027.
Schijven, J.F., Hoogenboezem, W., Hassanizadeh, M., Peters, J.H., 1999. Modeling removal of bacteriophages MS2 and PRD1 by dune recharge at Castricum, Netherlands.
Water Resour. Res. 35 (4), 1101–1111.
Schijven, J.F., De Bruin, H.A.M., Hassanizadeh, S.M., de Roda Husman, A.M., 2003b.
Bacteriophages and clostridium spores as indicator organisms for removal of pathogens by passage through saturated dune sand. Water Res. 37 (9), 2186–2194.
Schijven, J., Philip Berger, Ilkka Miettinen, 2003. Removal of pathogens, surrogates, indicators, and toxins using riverbank filtration. Riverbank Filtration. Springer
Netherlands: 73–116.
Schijven J.F, Hassnizadeh S.M, Husman, A.M, 2010. Vulnerability of unconfined aquifers
to virus contamination water research. 44, 4:1170–1181.
Schijven, J.F., Hassanizadeh, S.M., 2000. Removal of viruses by soil passage: overview ofmodeling, processes, and parameters. Crit. Rev. Environ. Sci. Technol. 30 (1),49–127.

Declaration of Competing Interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence the work reported in this paper.
References
Anders, R., Chrysikopoulos, C.V., 2005. Virus fate and transport during artificial recharge
with recycled water. Water Resour. Res. 41 (10), W10415.
Anders, R., Chrysikopoulos, C.V., 2006. Evaluation of the factors controlling the timedependent inactivation rate coefficients of bacteriophage MS2 and PRD1. Environ.
Sci. Technol. 40 (10), 3237–3242.
Bales, R.C., Li, S., Yeh, T.-C.J., Lenczewski, M.E., Gerba, C.P., 1997. Bacteriophage andmicrosphere transport in saturated porous media: forced-gradient experiment at
Borden, Ontario. Water Resour. Res. 33 (4), 639–648.
Bear, J., 1972. Dynamics of Fluids in Porous Materials. Society of Petroleum Engineers,
Dallas, TX, USA.
Bellou, M.I., Syngouna, V.I., Tselepi, M.A., Kokkinos, P., Paparrodopoulos, S.C.,
Vantarakis, A., Chrysikopoulos, C.V., 2015. Interaction of human adenoviruses andcoliphages with kaolinite and bentonite. Sci. Total Environ. 517, 86–95.
Bhattacharjee, S., Ryan, J.N., Elimelech, M., 2002. Virus transport in physically andgeochemically heterogeneous subsurface porous media. J. Contaminant Hydrol. 57(3–4), 161–187. https://doi.org/10.1016/S0169-7722(02)00007-4.
Bradford, S.A., Torkzaban, S., 2013. Colloid interaction energies for physically and chemically heterogeneous porous media. Langmuir 29 (11), 3668–3676. https://doi.org/
10.1021/la400229f.
Bradford, S.A., Torkzaban, S., 2015. Determining parameters and mechanisms of colloid
retention and release in porous media. Langmuir 31 (44), 12096–12105. https://doi.
org/10.1021/acs.langmuir.5b03080.
Bradford, S.A., Leij, F.J., Schijven, J., Torkzaban, S., 2017. Critical role of preferential
flow in field-scale pathogen transport and retention. Vadose Zone J. 16. https://doi.
org/10.2136/vzj2016.12.0127.
Chrysikopoulos, C.V., Aravantinou, A.F., 2012. Virus inactivation in the presence ofquartz sand under static and dynamic batch conditions at different temperatures. J.
Hazard. Mater. 233–234, 148–157. https://doi.org/10.1016/j.jhazmat.2012.07.002.
Chrysikopoulos, C.V., Sim, Y., 1996. One-dimensional virus transport homogeneous
porous media with time dependent distribution coefficient. J. Hydrol. 185, 199–219.
https://doi.org/10.1016/0022-1694(95)02990-7.
Chu, Y., Jin, Y., Baumann, T., Yates, M.V., 2003. Effect of soil properties on saturated andunsaturated virus transport through columns. J. Environ. Quality 32 (6), 2017–2025.
Dagan, G., Fiori, A., Jankovic, I., 2003. Flow and transport in highly heterogeneous
formations: 1. Conceptual framework and validity of first-order approximations. art
No. 1268. Water. Resour. Res. 39 (9), 12. https://doi.org/10.1029/2002WR001717.
Dillon, P., et al., 2008. A critical evaluation of combined engineered and aquifer treatment systems in water recycling. Water Sci. Technol. 753–762. https://doi.org/10.
2166/wst.2008.168.
Dillon, P., et al., 2010. Managed aquifer recharge: rediscovering nature as a leading edge
technology. Water Sci. Technol. 62 (10), 2338–2345. https://doi.org/10.2166/wst.
2010.444.
Foppen, J.W., Okletey, S., Schijven, J.F., 2006. Effect of goethite coating and humic acidon the transport of bacteriophage PRD1 in columns of saturated sand. J. ContaminantHydrol. 85 (3–4), 287–301. https://doi.org/10.1016/j.jconhyd.2006.02.004.
Freeze, R.A., 1975. A stochastic-conceptual analysis of one-dimensional groundwaterflow in nonuniform homogeneous media. Water Resour. Res. 11 (5), 725–741.
Garabedian, S., Gelhar, L., Celia, M., 1988. Large-scale dispersivae transport in aquifers:
Field experiments and reactive transport theory. Parsons Laboratory Report 315.
MIT, Cambridge.
Gerba, C.P., 1983. Virus survival and transport in groundwater. Developments in
Industrial Microbiology (USA).
Harvey, R.W., Kinner, N.E., MacDonald, D., Metge, D.W., Bunn, A., 1993. Role of physicalheterogeneity in the interpretation of small-scale laboratory and field observations ofbacteria, microbial-sized microsphere, and bromide transport through aquifer sediments. Water Resour. Res. 29 (8), 2713–2721.
Harvey, R.W., Ryan, J.N., 2004. Use of PRD1 bacteriophage in groundwater viral transport, inactivation, and attachment studies. FEMS Microbiol. Ecol. 49 (1), 3–16.
https://doi.org/10.1016/j.femsec.2003.09.015.
Hornstra, L.M., Schijven, J.F., Waade, A., Prat, G.S., Smits, F.J.C., Cirkel, G., Stuyfzand,P.J., Medema, G.J., 2018. Transport of bacteriophage ms2 and prd1 in saturated dunesand under suboxic conditions. Water Res. 139, 158–167.

10

Journal of Hydrology 578 (2019) 124082

S. Torkzaban, et al.
Schijven, J.F., Sadeghi, G., Hassanizadeh, S.M., 2016. Long-term inactivation of bacteriophage PRD1 as a function of temperature, pH, sodium and calcium concentration.
Water Res. 103, 66–73.
Shen, C., et al., 2012. Theoretical and experimental investigation of detachment of colloids from rough collector surfaces. Colloids and Surfaces A: Physicochem. Eng.
Aspects 410, 98–110. https://doi.org/10.1016/j.colsurfa.2012.06.025.
Sidhu, J.P.S., et al., 2015. Pathogen decay during managed aquifer recharge at four sites
with different geochemical characteristics and recharge water sources. J. Environ.
Quality 44 (5), 1402–1412. https://doi.org/10.2134/jeq2015.03.0118.
Sim, Y., Chrysikopoulos, C.V., 1996. One-dimensional virus transport in porous media
with time dependent inactivation rate coefficients. Water Resour. Res. 32 (8),2607–2611. https://doi.org/10.1029/96WR01496.
Šimůnek, J., Jarvis, N.J., Van Genuchten, M.T., Gärdenäs, A., 2003. Review and comparison of models for describing non-equilibrium and preferential flow and transport
in the vadose zone. J. Hydrol. 272 (1), 14–35.
Simunek, J., van Genuchten, M.T., Sejna, M., 2016. Recent developments and applications of the HYDRUS computer software packages. Vadose Zone J. 15 (7), 25. https://doi.org/10.2136/vzj2016.04.0033.
Sinreich, M., Flynn, R., 2011. Comparative tracing experiments to investigate epikarst
structural and compositional heterogeneity. Speleogenesis Evol. Karst Aquifers 10,
60–67.
Sudicky, E., MacQuarrie, K., 1989. Behaviour of biodegradable organic contaminants in
random stationary hydraulic conductivity fields, Kobus, HE, and Kinzelbach, W.,
Contaminant Transport in Groundwater, International Symposium on Contaminant
Transport in Groundwater, Stuttgart, pp. 307–315.
Torkzaban, S., Hassanizadeh, S.M., Schijven, J.F., De Bruin, H.A.M., De Roda Husman,
A.M., 2006. Virus transport in saturated and unsaturated sand columns. Vadose Zone

J. 5 (3), 877–885. https://doi.org/10.2136/vzj2005.0086.
Torkzaban, S., Bradford, S.A., Wan, J., Tokunaga, T., Masoudih, A., 2013. Release ofquantum dot nanoparticles in porous media: role of cation exchange and aging time.
Environ. Sci. Technol. 47 (20), 11528–11536. https://doi.org/10.1021/es402075f.
Torkzaban, S., Bradford, S.A., Vanderzalm, J.L., Patterson, B.M., Harris, B., Prommer, H.,
2015. Colloid release and clogging in porous media: effects of solution ionic strengthand flow velocity. J. Contam. Hydrol. 181, 161–171.
Tufenkji, N., Elimelech, M., 2004. Correlation equation for predicting single-collector
efficiency in physicochemical filtration in saturated porous media. Environ. Sci.
Technol. 38 (2), 529–536.
Unc, A., Goss, M.J., 2004. Transport of bacteria from manure and protection of water
resources. Appl. Soil Ecol. 25 (1), 1–18. https://doi.org/10.1016/j.apsoil.2003.08.
007.
Vadillo-Rodriguez, V., Busscher, H.J., Norde, W., de Vries, J., van der Mei, H.C., 2004.
Atomic force microscopic corroboration of bond aging for adhesion of Streptococcus
thermophilus to solid substrata. J. Colloid Interface Sci. 278 (1), 251–254. https://doi.org/10.1016/j.jcis.2004.05.045.
Wang, Y., Bradford, S.A., Šimůnek, J., 2014. Physicochemical factors influencing thepreferential transport of in soils. Vadose Zone J. 13 (1). https://doi.org/10.2136/
vzj2013.07.0120.
Xu, L.-C., Logan, B.E., 2006. Adhesion forces between functionalized latex microspheres
and protein-coated surfaces evaluated using colloid probe atomic force microscopy.
Colloids and Surfaces B: Biointerfaces 48 (1), 84–94.
Yao, K.-M., Habibian, M.T., O'Melia, C.R., 1971. Water and waste water filtration.
Concepts and applications. Environ. Sci. Technol. 5 (11), 1105–1112. https://doi.
org/10.1021/es60058a005.

11