A random walk solution for modeling solute transport with network reactions and multi-rate mass transfer in heterogeneous systems: Impact of biofilms

A Random Walk Solution for Modeling Solute Transport with Network Reactions and Multi-Rate Mass Transfer in Heterogeneous Systems: Impact of Biofilms Christopher V. Henria,, Daniel Fern`andez-Garciaa, aDepartment of Geotechnical Engineering and Geosciences, Universitat Polit`ecnica de Catalunya, UPC-Barcelona Tech, Jordi Girona 1-3, 08034 Barcelona, Spain. Abstract The interplay between the spatial variability of aquifer properties, mass transfer and chemical reactions often complicates reactive transport simulations. It is well documented that hydro-biochemical properties are ubiquitously heterogeneous and that rate-limited mass transfer typically leads to the conceptualization of an aquifer as a multi-porosity system. Within this context, chemical reactions taking place in mobile/immobile water regions can be substantially different between each other. This paper presents a particle-based method that can efficiently simulate heterogeneity, network reactions and multi-rate mass transfer. The approach is based on the development of transition probabilities that describe the likelihood that particles belonging to a given species and mobile/immobile domain at a given time will be transformed into another species and mobile/immobile domain afterwards. The joint effect of mass transfer and sequential degradation is shown to be non-trivial. A characteristic double peak of daughter products occurs when the degradation capacity in the immobile domain is relatively small. This late rebound of concentrations is not driven by any change in the flow regime (e.g., pumping ceases in the pump-and-treat remediation strategy) but due to the natural interplay between mass transfer and chemical reactions. To illustrate that the method can simultaneously represent mass transfer, spatially varying properties and network reactions without numerical problems, we have simulated the degradation of tetrachloroethylene (PCE) in a three-dimensional fully heterogeneous aquifer subjected to Email address: [email protected] (Daniel Fern`andez-Garcia ) Preprint submitted to Elsevier May 12, 2015 rate-limited mass transfer. Two types of degradation modes were considered to compare the effect of an active biofilm with that of clay pods present in the aquifer. Results of the two scenarios display significantly differences. Biofilms that promote the degradation of compounds in an immobile region are shown to significantly enhance degradation, rapidly producing daughter products and less tailing. Keywords: network-reaction, multirate mass-transfer, random-walk particle-tracking, heterogeneities. 1. Introduction The assessment of groundwater polluted systems requires efficient and trustable predictive models. By efficient, we usually mean a small computational cost, and by trustable the proper representation of all key processes controlling the fate and transport of contaminants. Unfortunately, these two desirable properties contrast with the common adage opposing efficiency and complexity in transport simulations. Model complexity often arises from the need to jointly incorporate, in three dimensions, the spatial variability of aquifer properties and numerous chemical reactions into the same multispecies reactive transport model. Thus, properties such as the hydraulic conductivity and the different degradation rates can vary several orders of magnitude in an aquifer (e.g., Rubin, 2003; Fenell et al., 2001; Sandrin et al., 2004). Even though the description of the spatial variability of all these properties at high resolution is crucial for making contaminant predictions (e.g., Feehley et al., 2000; Salamon et al., 2007; Riva et al., 2008; Llopis- Albert et al., 2009), its implementation in transport models typically leads to numerical problems. Contaminants in the subsurface are also affected by chemical reactions. Among them, network reactions have been used to model a large variety of contaminants. This includes the degradation of chlorinated solvents (e.g., Clement, 1997, 2001), the decay of radioactive species (e.g., Painter et al., 2007), and the transformation of pesticides, organic phosphates and nitrogen in the environment (e.g., van Genuchten, 1985; Mishra and Mishra, 1991; Vishwanathan et al., 1998). When contaminant concentrations are small, i.e., less than the Michaelis half-saturation constant in the Monod or Michaelis- Menten enzyme kinetic model, the microbial biotransformation rates can be described by pseudo-first-order reaction rates (e.g., Bouwer et al., 1981; 2 Vogel et al., 1987; Haston and McCarty, 1999; Burnell et al., 2014). In this context, organic chlorinated solvents are often described by first-order reaction chains schematically described by A→B→C→D, meaning that species A is transformed into species B, B into C and so on. The quantification of the risk posed by these contaminants is not a trivial problem (Benekos et al., 2006; Henri and Fern`andez-Garcia, 2014). The degradation products can constitute new noxious chemical compounds not necessarily less toxic than its parent product. Even though analytical solutions (e.g., Sun et al., 1999; Zhang and Woodbury, 2002; Sun and Buscheck, 2003; Falta et al., 2007) have been routinely employed by decision makers to efficiently manage the risk posed by these contaminants, the spatial variability of aquifer properties, always observed in natural systems, have seriously questioned the application of such simple approaches (Adrian et al., 1994; Miralles-Wilhelm and Gelhar, 1996; Scholl, 2000; Cunningham and Fadel, 2007; Henri and Fern`andez-Garcia, 2014). The presence of low permeability regions where contaminants can be temporarily trapped by diffusion typically reduces the efficiency of in situ cleanup technologies (e.g., Soga et al., 2004; Stroo et al., 2012). To represent this situation in a macroscopic transport model, the porous medium is typically conceptualized as a multi-porosity system that partitions the domain into one region formed by mobile water where advection and dispersion can take place, and another region with relatively immobile water where transport is diffusion limited (van Genuchten and Wierenga, 1976; Neretnieks, 1980; Cunningham et al., 1997; Carrera et al., 1998; Zinn and Harvey, 2003). This conceptual model has received great attention in the last decades for its success in reproducing tailing and anomalous transport (Harvey and Gorelick, 2000; Huang and Hu, 2000; Li and Brusseau, 2000; Haggerty et al., 2000). Albeit mass transfer models were originally formulated with a single mass transfer coefficient, i.e., a double porosity system (van Genuchten and Wierenga, 1976; Neretnieks, 1980; Harvey and Gorelick, 2000), this particular model has shown drastic limitations in characterizing the long-term behavior of solute transport (Young and Ball, 1995; Haggerty and Gorelick, 1995; Haggerty et al., 2000). The main reason is that the diverse mineralogical composition of a real soil and the complex spatial variability of aquifer properties can produce a mixture of mass transfer processes occurring over multiple scales. To overcome this limitation, the multirate mass transfer model was later on developed by Haggerty and Gorelick (1995), who extended the double porosity model to multiple immobile domains. The multirate mass transfer model 3 represents a complete and practical answer to the complexity of a natural groundwater polluted system. For instance, by choosing appropriate parameter values, this model can simulate a large variety of diffusion processes such as diffusion into cylinders, spheres, layers, rock matrices, and others (Haggerty and Gorelick, 1995; Carrera et al., 1998). Moreover, various authors (e.g. Zinn and Harvey, 2003; Liu et al., 2004; Fern`andez-Garcia et al., 2009) have demonstrated, that solute transport can be properly upscaled using a multirate mass transfer model. The incorporation of all theses processes highly complicates numerical simulations. Reactive transport codes based on Eulerian methods (finite differences or finite-elements) typically encounter numerical problems and suffer from computational burden when dealing with strong heterogeneities (Salamon et al., 2006; Boso et al., 2013). These limitations have been largely emphasized by the recent need to conduct integrated risk analysis studies with many uncertain parameters at high resolution (Maxwell et al., 2008; De Barros et al., 2009; Benekos et al., 2006). In this context, Particle Tracking Methods (PTMs) constitute an efficient numerical alternative to simulate reactive transport (Kitanidis, 1994; Salamon et al., 2006; Henri and Fern`andez- Garcia, 2014). Even though a large variety of methods exist to simulate rate-limited mass transfer processes with particle tracking (Benson and Meerschaert, 2009; Delay and Bodin, 2001; Dentz and Berkowitz, 2003; Salamon et al., 2006; Tsang and Tsang, 2001), this method is still limited in the type of chemical reactions available, which include sorption (Tompson, 1993; Valocchi and Quinodoz, 1989; Michalak and Kitanidis, 2000), radioactive decay (Wen and G´omez-Hern´andez, 1996; Painter et al., 2007), first-order network reactions (Burnell et al., 2014; Henri and Fern`andez-Garcia, 2014), and simple bimolecular reactions (Benson and Meerschaert, 2008; Ding et al., 2013; Edery et al., 2009, 2010; Paster et al., 2014) among others. None of the methods available nowadays supports multi-porosity systems with network reactions in three-dimensional randomly heterogeneous porous media. We present a new particle tracking method that efficiently handles hydrobiogeochemical spatial variability, multi-rate mass transfer processes and first-order network reactions. The approach is based on the development of transition probabilities that describe the likelihood that particles belonging to a given species and mobile/immobile domain at a given time will be transformed into another species and mobile/immobile domain afterwards. These probabilities are derived from the solution matrix of the spatial moments governing equations. The organization of the paper is as follows. First, 4 sections 2 and 3 respectively present the governing transport equations and the development of transition probabilities. These probabilities are then used in section 4 to incorporate mass transfer processes and network reactions into a random walk model. Finally, a 3D high-resolution synthetic example is presented to investigate the impact that the different sources of variability have on the fate and transport of tetrachloroethylene (PCE) and its degradation products in a groundwater polluted system. 2. The Multirate Model with Network Reactions 2.1. Conceptual Model The multirate mass transfer model has been extensively presented in the literature (Haggerty and Gorelick, 1995; Haggerty et al., 2000; Donado et al., 2009). Here, we extend the model for use in modeling contaminant transport with network reaction systems. The porous medium is conceptualized as a multi-porosity system consisting of one mobile water region where solute moves by advection and dispersion, and any number of immobile water regions where solute transport is diffusion limited. A series of mass transfer equations describe the mass exchange between the mobile and any number of immobile regions. Considering this conceptualization of the porous medium, Figure 1 shows a schematic representation of a chain reaction (A→B→C) in this system. This figure represents several physical and biochemical processes that occur simultaneously in porous media. At the large scale, contaminants can diffuse into clay layers/lenses or get transferred into low conductivity regions by slow advection. In this context, it is generally observed that biotransformation rates in the mobile and immobile regions can be substantially different. The main reason is that biotransformation rates and bacterial activity largely depends on the clay content, being often smaller in confining beds than in more permeable systems (Chapelle, 2001). Because bacteria generally have diameters that range between 0.1 and 1 microns, the small effective porosity of clays typically restricts the ability of bacteria to move and reproduce effectively. Also, the natural occurrence of preferential flow channels in porous media (typically represented as a mobile region in the multirate model) favors the movement of groundwater and dissolved elements through certain pathways, which typically harbor larger bacterial densities and microbial activities compare to the adjacent porous media (Pivetz and Steenhuis, 1995; Mallawatantri et al., 1996; Rubol et al., 2014). In fact, Vinther et al. (1999) and Bundt et al. (2001) found 5 that both substrate availability and nutrient supply are largest in preferential flow paths, enhancing bacterial biomass and associated microbial processes. Similar processes occur at a smaller scale where contaminants may diffuse into stagnant water (intraaggregate or dead-end pores) or/and inside biofilms attached to the aquifer soil surface. In this case, the rate of contaminant degradation mainly occurs in the active biofilm rather than in the bulk aqueous phase (Rittmann and McCarty, 1980; Baveye and Valocchi, 1989; Cunningham and Mendoza-Sanchez, 2006). Based on these observations, the reactive multirate model proposed here considers that network reactions can simultaneously occur over multiple scales with different biotransformation rates in the mobile and immobile regions. 2.2. Governing Equations We consider a network reaction system formed by Nsspecies moving through a mobile domain of porosity φ0and Nim immobile domains of porosity (φ1, ..., φNim ). Let us denote the biotransformation rate associated with species iand jin the `domain by rij` =yijkj`φ`cj`,(1) where rij` is the production of mass of the ith species per unit of time and aquifer volume [ML−3T−1] due to the degradation of the jth species in the ` domain, ki` [T−1] is the first-order contaminant destruction rate constant of the ith species in the `domain, yij [M M−1] is the effective yield coefficient for any reactant or product pair (i, j). This coefficient is a stoichiometric coefficient defined as the ratio of mass of species igenerated to the amount of mass of species jconsumed. The yield coefficients yii is equal to −1 and represent the first-order decay of the ith species. Similar reaction terms have been presented by many authors (Clement, 1997, 2001; Sun et al., 1999; Falta et al., 2007). The multirate mass transfer model with network reactions can be written as the following system of equations φ0Ri0 ∂ci0 ∂t + Nim X l=1 φ`Ri` ∂ci` ∂t =L(ci0) + Ns X j=1 Nim X `=0 rij`, i = 1, Ns,(2) φ`Ri` ∂ci` ∂t =φ`α0 i` (ci0−ci`) + Ns X j=1 rij`, ` = 1, Nim, i = 1, Ns.(3) 6 Without the chemical reaction term rij`, these equations form the standard multirate mass transfer model (Haggerty and Gorelick, 1995) that describes advective-dispersive transport with rate-limited mass transfer between a mobile domain and any number of immobile domains for each species. In these equations, the variable ci0[M L−3] is the concentration of the ith species in the mobile domain (denoted always by the subscript index `= 0), ci` [M L−3] is the concentration of the ith species in the `th immobile domain (`= 1, ..., Nim), Ri0(dimensionless) is the retardation factor of the ith species in the mobile domain, and Ri` (dimensionless) is the retardation factor of the ith species in the `th immobile domain (`= 1, ..., Nim). Sorption is considered in local equilibrium (linear isotherm), and L(c) is the transport operator of the mobile concentrations defined by L(c) = ∇·(φ0D∇c)−∇·(qc),(4) where q[L T−1] is the groundwater flux, and Dis the dispersion tensor [L2T−1]. The first equation (2) is actually the mass balance associated with any of the species involved in the network reaction system, and equation (3) describes the mass transfer of the ith species between the mobile domain and the `th immobile domain. We have assumed that only aqueous concentrations are susceptible to undergo chemical reactions, i.e., no biodegradation in the sorbed phase occurs. Nevertheless, we note that other situations can be simulated by properly redefining the degradation rates (van Genuchten, 1985). 3. Development of Transition Probabilities Transition probabilities denote the probability that a system that was in a given state at time t= 0 jumps to another state at some later time t > 0 (e.g., Ross, 2003). In subsurface hydrology, this concept has been used in the past to simulate solute transport with sorption/desorption processes (Valocchi and Quinodoz, 1989; Michalak and Kitanidis, 2000), rate-limited mass transfer processes (Salamon et al., 2006), and kinetic network reactions (Henri and Fern`andez-Garcia, 2014). Following Michalak and Kitanidis (2000) and Henri and Fern`andez-Garcia (2014), transition probabilities can be determined from the evolution of the zeroth spatial moments of the solute plume. This is the procedure employed here. Let us express the governing system of equations (2)-(3) in terms of the total densities, ρi`, defined as the total amount of aqueous and sorbed mass 7 of a given species iper unit volume in the `domain, i.e., ρi` =φ`Ri`ci`. From this, the transport equation reads as Nim X `=0 ∂ρi` ∂t =Lρi0 φ0Ri0+ Ns X j=1 Nim X `=0 Kij`ρj`, i = 1, Ns,(5) ∂ρi` ∂t =αi`βi`ρi0−αi`ρi` + Ns X j=1 Kij`ρj`, ` = 1, Nim, i = 1, Ns,(6) where Kij` =yijkj`/Rj`,βi` is the field capacity coefficient of the ith species in the `immobile domain, and αi` is the apparent mass transfer coefficient of the ith species in the `immobile domain, βi` =φ`Ri` φ0Ri0 ,(7) αi` =α0 i` Ri` .(8) Let us also define the zeroth spatial moment of the solute plumes associated to species iin the mobile or any immobile domain `by, µi`(t) = ZV ρi`(x, t)dV, ` = 0, Nim, i = 1, Ns,(9) where Vis the volume of the domain. Knowing that a particle located at position xtat time tcan be seen as an infinitely small plume (Kitanidis, 1994), the total density of the particle in either the mobile or immobile domain can be represented by ρi`(x, t) = µi`(t)δ(x−xt), ` = 0, Nim, i = 1, Ns,(10) The governing equations of the spatial moments can be derived by integrating by parts the transport equations given by (2)-(3). This leads to Nim X `=0 dµi`(t) dt = Ns X j=1 Nim X `=0 Kij`(xt)µj`(t),∀i= 1, Ns,(11) 8 dµi` dt =αi`βi`µi0−αi`µi` + Ns X j=1 Kij`µj`, ` = 1, Nim, i = 1, Ns.(12) Let us now define the vector mtof size n=Ns(1 + Nim) as mt=   m1(t) . . . mNs(t)   ,where mi(t) =      µi0(t) µi1(t) . . . µiNim (t)      .(13) The system of equations (11)-(12) can then be expressed in matrix format as Admt dt =B(xt)mt.(14) Note that the subscript xtexpresses that the matrix is evaluated at the particle position xt. The Amatrix is an n×nblock diagonal matrix defined as A=     A0 11 0··· 0 0 A0 22 .... . . . . .......0 0··· 0 A0 NsNs      ,(15) where 0is the zero matrix of size (1+Nim)×(1+Nim), and A0 ii {i= 1, ..., Ns}, is the following (1 + Nim)×(1 + Nim) constant matrix A0 ii =     1 1 ··· 1 0 1 0 0 . . .......0 0··· 0 1      .(16) The B(xt) matrix is the following n×nblock matrix evaluated at the xt location, B(xt) =    B0 1 1 ··· B0 1Ns . . ..... . . B0 Ns1··· B0 NsNs   xt ,(17) 9 The variable DaII can be seen as the second Damkh¨oler number defined as the ratio of the chemical reaction rate to the mass transfer rate. The variable χis the ratio between the immobile and the mobile decay rate. As shown in Appendix A, in this case, the transition probability matrix Ptcan be explicitly determined. From this, we evaluate the influence of DaII and χ on the temporal evolution of Ptassociated with a particle initially belonging to species 1 and the mobile domain. The effect that the field capacity β has on the transition probabilities is similar to the effect induced by the Damkh¨oler number and is therefore not shown. Figure 2 shows the influence of mass transfer on the temporal evolution of Ptfor χ= 1 and β= 10. When the mass transfer rate is larger than the decay rate in the mobile domain (i.e., DaII <1), the probability to remain in the mobile domain (Figure 2a) drops to an early equilibrium between the mobile and the immobile domain. After this, biochemical reactions start dominating and this probability decays with time to almost zero. In accordance with this result, the probability that species 1 is in the immobile domain after time t increases to a plateau until chemical reactions take place at τ≈1 (Figure 2b). Consequently, the probability that this particle is in the immobile domain increases with the mass transfer rate (inverse of DaII). Interestingly, at early times, mass transfer is still not active and the probability that species 1 turns into 2 in the mobile domain does not depend on DaII (Figure 2c). This probability increases linearly with time due to the degradation of species 1 in the mobile domain. When DaII <1 and mass transfer starts to take place there are less particles of species 1 to degrade in the mobile domain and the probability stabilizes. With time, this effect vanishes and the probability increases again linearly with time until species 2 starts to degrade. A decay difference between the mobile and immobile domain (χ6= 1) can also have relevant consequences on the temporal evolution of Pt. This is shown in Figure 3 for DaII = 1 and β= 10 (rate-limited mass transfer). The probability that the particle still belongs to species 1 in the mobile domain after a time tis shown in Figure 3a. Results show that when degradation in the immobile domain is smaller than that of the mobile domain (χ < 1), for example due to the presence of an aquitard, the early equilibrium between the mobile and the immobile domain discussed previously is also observed. However, when degradation is higher in the immobile domain (χ > 1), for example due to the presence of biofilms at the pore-scale, the mass transferred into the immobile domain will be rapidly consumed (see Figure 3b), preventing the previously observed mass transfer equilibrium 16 (Figure 3a). During pump-and-treat remediation strategies, it is often observed that once pumping ceases, a rebound of concentrations at the well takes place (e.g., De Barros et al., 2013). Figure 3c shows that in a double-porosity system with network reactions this effect can also occur without any change in the pumping regime. We call this effect ”natural rebound” which is explained as it follows. When degradation is active in the immobile domain, the parent species transferred into the immobile domain will be transformed into degradation products. Once the reaction chain in the mobile domain has occurred, these products will be allowed to back diffuse into the mobile domain causing the rebound of concentrations. This explains the double peak observed in Figure 3c for the degradation product. When degradation is not active in the immobile domain, a double peak is also observed at a later time (χ= 0). In this case, once the chain reaction in the mobile domain has occurred, a second chain reaction in the mobile domain can be triggered by the release of the parent species previously stored in the immobile domain due to back-diffusion. 6. An Example of Application: Effect of the parameters spatial variability 6.1. Problem Setup The main advantage of our method is the possibility to simultaneously represent mass transfer, spatially varying properties (heterogeneity) and network reactions without numerical problems. To illustrate this, we consider a three-dimensional high-resolution synthetic aquifer initially contaminated by tetrachloroethylene (PCE) and affected by rate-limited mass transfer and degradation. With time, PCE is sequentially transformed into TCE (trichloroethylene), TCE into DCE (cis-Dichloroethylene), and DCE into VC (vinyl chloride). Groundwater flow takes place in a heterogeneous hydraulic conductivity field obtained from a single realization of a sequential Gaussian simulation. The natural logarithm of the hydraulic conductivity (ln K) is described by a random function of mean 3.55 m/d and an isotropic variogram of range 10.0 m and variance of 2.5. The domain is a rectangular block of 120×100×40 m3, descritized into 300×250×100 cubic grid cells of size 0.4 m (Figure 4). The flow is driven by a hydraulic gradient of 0.01 oriented along the x-axis and solved using the finite difference code Modflow (Harbaugh et al., 2000). The transport parameters adopted are summarized 17 in Table 1. Porosity, local dispersivities and retardation factors are always considered homogeneous. The source of contamination is represented by a PCE instantaneous point injection in the mobile domain (x=20,y=50,z=20) of unit mass equally partitioned into 100,000 particles. Concentration breakthrough curves of all species (BTCs) were recorded at two control planes located at 1 and 5 variogram ranges from the injection location (see Figure 4). Rate-limited mass transfer is simulated using the spherical diffusion model presented by Rao et al. (1980). Haggerty and Gorelick (1995) showed that this model and the multi-rate model are mathematically equivalent provided that the series of porosities and mass transfer coefficients are specifically chosen. For practical purposes, the series has to be truncated to a discrete number. Haggerty (2009) explained that the truncation of the multi-rate series becomes an acceptable approximation when the final term is defined in an appropriate manner. Table 3 shows the series of mass transfer and porosities for spherical diffusion. The number of terms used is 10. The total field capacity (βtot) is fixed to 1.0 and the apparent pore diffusion (Dp/a2) is set to 0.0023. These properties were chosen from the review of mass transfer data given by Haggerty et al. (2004). Mass transfer parameters are summarized in Table 2. The degradation rates of contaminants in aquifers can substantially vary in space (Allen-King et al., 2006) due to, for instance, changes in the bacteria activity responsible for biodegradation (e.g., Fenell et al., 2001; Sandrin et al., 2004). To represent this, we consider that a perfect negative linear correlation between ln Kand the decay rates exist (Miralles-Wilhelm and Gelhar, 1996; Cunningham and Fadel, 2007). We admit this way that small decay rates relate to high ln Kvalues and vice versa. This can happen for instance when water fluxes can effectively deplete electron acceptors/donors in the porous medium (Cunningham and Fadel, 2007). The minimum and maximum values of PCE, TCE, DCE and VC decay rates were defined based on the range of possible values provided by the U.S. Environmental Protection Agency as a support to decision makers (Environmantal Protection Agency, 1999, 2002). Thus, decay rates reflect those obtained in several field sites and laboratory observations. The statistics of the decay rates are given in Table 2. The spatial distribution of the lnK and the decay rates in the mobile domain is shown in Figure 4. Two different conceptual models for the degradation in the immobile domain were considered. The first model assumes that degradation is not taking place in the immobile domain. This repre- 18 sents for instance the presence of small clay layers or pods in the aquifer, preventing bacteria to move and reproduce effectively. The second model considers the enhancement of degradation due to the existence of an active biofilm at the pore scale. In this case, the decay rate in the immobile domain is assumed to be 10 times larger than in the mobile domain. 6.2. Results and Discussion The distribution of mobile and immobile particles simulated with the proposed random walk method is respectively shown in Figure 5 and 6 for two different times. The species compound is denoted by the particle color. The particle size is proportional to the log of the cancer potency factor to visually display not only the density of particles (concentrations) but also the potential threat that these contaminants pose to human health. Important differences between the two conceptual models can be distinguished. The biofilm model clearly displays an enhancement of degradation which now mostly occurs in the immobile domain. Comparing the two models, we see that, in the biofilm model, a larger portion of the PCE has been already transformed into daughter products at time t=160 d. This means that daughter products can now be produced at earlier times and closer to the source area than expected from the mere diffusion of products into clay regions. This has important consequences for risk assessment. In heterogeneous porous media without local mass transfer effects, Henri et al. (submitted) have demonstrated that the hot spot of the risk posed by a chemical mixture (co-existence of the original pollutants and their daughter products) depends on the joint effect of degradation, advection and toxicity. Results here show that the conceptualization of the immobile domain as an active degradation region can largely complicate the quantification and interpretation of human health risk. The concentration breakthrough curves of all the species obtained at two different control planes during the simulations are shown in Figure 7. Results show that the biofilm model and the clay model give also very distinct signals. The clay model displays breakthrough curves with long tails due to ratelimited mass transfer. In this model, particles can be temporally trapped in the immobile domain by diffusion without undergoing degradation in this region. These particles will be capable to back diffuse into the mobile region at a later time. However, in the biofilm model, these trapped particles can be transformed into daughter products. This process prevents tailing and increases the amount of degradation products that can be transferred later on by diffusion into the mobile domain. This explains why the breakthrough 19 curves in the biofilm model display more peaked distributions with more mass and less tailing. 7. Conclusions The interaction between the spatial variability of aquifer properties, mass transfer and chemical reactions often complicates reactive transport simulations. It is well documented that hydro-biochemical properties are ubiquitously heterogeneous and that rate-limited mass transfer typically leads to the conceptualization of an aquifer as a multi-porosity system. In this conceptual model, contaminants can be transferred to a number of immobile water regions so as to represent different phenomena observed at multiple scales, i.e., diffusion into stagnant water at the pore-scale, diffusion into biofilms attached to soil surface, diffusion into low permeability inclusions at the centimeter to meter scale, and diffusion into aquitards or the rock matrix of a fractured system. Importantly, the chemical reactions taking place in these mobile/immobile water regions can be substantially different between each other. Along this line, we have presented a random walk solution that is capable to efficiently simulate multi-rate mass transfer and first-order network reactions in heterogeneous porous media without restrictions in the spatial variability of biochemical and hydrodynamic properties. First-order rate coefficients vary in space and the type of water region involved. The approach is based on the development of transition probabilities that describe the likelihood that particles belonging to a given species and mobile/immobile domain at a given time will be transformed into another species and mobile/immobile domain afterwards. This is important for assessing the risk posed by a large variety of chemical systems that otherwise suffer from numerical problems in dealing with heterogeneities during reactive transport modeling, e.g., the degradation of chlorinated solvents, the decay of radioactive species, and the transformation of pesticides, organic phosphates and nitrogen in the environment. The method is limited to first-order network reactions. The approach is used to investigate the joint effect of network reactions and mass transfer on the spatiotemporal behavior of the sequential degradation of tetrachloroethylene(PCE). Transition probabilities show that a double peak of daughter products can occur when the degradation capacity in the immobile domain is relatively small. This late rebound of concentrations is not driven by any change in the flow regime (e.g., pumping ceases in the pump-and-treat remediation strategy) but due to the natural interplay be- 20 tween mass transfer and chemical reactions. To illustrate that the method can simultaneously represent mass transfer, spatially varying properties and network reactions without numerical problems, we have simulated the degradation of PCE in a three-dimensional fully heterogeneous aquifer subjected to rate-limited mass transfer. Two types of degradation modes were considered to compare the effect of an active biofilm with that of clay pods in the aquifer. Both hydraulic conductivity and biochemical parameters were considered spatially variable and described at high resolution. Results of the two scenarios display significant differences. Biofilms that promote the degradation of compounds in an immobile region are shown to significantly enhance degradation, rapidly producing daughter products and less tailing. Acknowledgements The authors acknowledge the financial support provided by the Spanish Ministry of Science and Innovation through the SCARCE Consolider-Ingenio 2010 program (reference CSD2009-00065) and project FEAR (CGL2012- 38120). 21 Figure 1: Illustration of mass transfer occurring at different scales and rates in case of a sequential three-species reactive system (A →B→C). The thickness of the bicolor arrow representing biochemical reaction is proportional to the expected intensity of the reaction. The bottom right magnifying frame illustrates the heterogeneity in the grain distribution and related contaminant diffusion from/into aggregates. The top right magnifying frame illustrates the potential diffusion-rate limited transfer of contaminant toward an active biofilm and the subsequent enhanced degradation. 22 10-5 10-3 10-1 101103 10-3 10-2 10-1 100 Transition Probability [-] DaII = ¥ DaII = 101 DaII = 100 DaII = 10-1 DaII = 10-2 10-5 10-3 10-1 101103 10-3 10-2 10-1 100 10-5 10-3 10-1 101103 t [-] 10-5 10-4 10-3 10-2 10-1 100 Transition Probability [-] 10-5 10-3 10-1 101103 t [-] 10-5 10-4 10-3 10-2 10-1 100 Final State: Species 1, immobile zone Final State: Species 1, mobile zone Final State: Species 2, mobile zone Final State: Species 2, immobile zone a. b. c. d. Figure 2: Impact of the second Damk¨ohler number on the transition probability of a particle initially being species 1 in the mobile domain: (a) probability to stay in the initial state after a time t; (b) probability to be transferred to species 1 in the immobile domain; (c) probability to turn into species 2 in the mobile domain; and (d) probability to turn into species 2 in the immobile domain. 23 10-4 10-2 100102104 10-3 10-2 10-1 100 Transition Probability [-] = 0.0 = 10-1 = 100 = 101 = 102 10-4 10-2 100102104 10-3 10-2 10-1 100 10-4 10-2 100102104 t [-] 10-4 10-3 10-2 10-1 100 Transition Probability [-] 10-4 10-2 100102104 t [-] 10-4 10-3 10-2 10-1 100 Final State: Species 1, immobile zone Final State: Species 1, mobile zone Final State: Species 2, mobile zone Final State: Species 2, immobile zone a. b. c. d. Figure 3: Impact of the ratio between decay rate in the mobile and immobile domain on the transition probability of a particle initially being species 1 in the mobile domain: (a) probability to stay in the initial state after a time t; (b) probability to be transferred to species 1 in the immobile domain; (c) probability to turn into species 2 in the mobile domain; and (d) probability to turn into species 2 in the immobile domain. 24 Figure 4: Simulation setup displaying the lognormal distribution of the hydraulic conductivity and first-order degradation rates of PCE, TCE, DCE and VC. 25 which give the matrix A−1B A−1B=    −Da−1 II β−1Da−1 II 0 0 Da−1 II β−Da−1 II −χ0 0 1 0 −Da−1 II β−1Da−1 II 0χ Da−1 II β−Da−1 II −χ     . (A.7) The eigensytem of a such matrix can be obtained analytically. The four eigenvalues are given by λ=λ1λ1λ2λ2Twhere λ1=1 2qDa−2 II β2+ 2 Da−1 II Da−1 II −χ+ 1β+Da−1 II +χ−12 +1 2(−β−1) Da−1 II −χ 2−1 2, λ2=−1 2qDa−2 II β2+ 2 Da−1 II Da−1 II −χ+ 1β+Da−1 II +χ−12 +1 2(−β−1) Da−1 II −χ 2−1 2. The matrix of eigenvectors is defined by: S=    S110S120 S210S220 S310S320 1 0 1 0     (A.8) where S1i=γiDa−1 II Da−1 II +λi+χDa−1 II β+λi+ 1Ψi S2i=γi Da−1 II +λi+χΨi S3i=Da−1 II β+λi+ 2χ+Da−1 II +λiDa−1 II Ψi after the definition of Ψiand γias: Ψi=−Da−1 II β−λi−1χ2+Da−2 II β2+ ((β−1) λi+ 2 β−1) Da−1 II +λi+ 1χ +Da−2 II β, 32 γi=Da−1 II β+λi+ 1χ+ (1 + (β+ 1) λi)Da−1 II +λi2+λiβ(χ+ 1) Da−2 II 33 Adrian NR, Robinson JA, Suflita JM. Spatial Variability in Biodegradation Rates as Evidenced by Methane Production from an Aquifer Spatial Variability in Biodegradation Rates as Evidenced by Methane Production from an Aquifer. Applied and environmental microbiology 1994; 60(10): 3632- 3639. Allen-King RM, Divine DP, Robin MJL, Alldredge JR, Gaylord DR. Spatial distributions of perchloroethylene reactive transport parameters in the Borden Aquifer, Water Resour. Res. 2006; 42(1): p. W01413. Anderson, E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J., Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A. and Sorensen, D. (1999), LAPACK Users’ Guide, Third Ed., Society for Industrial and Applied Mathematics. Baveye, P., and A. Valocchi (1989), An evaluation of mathematical models of the transport of biologically reacting solutes in saturated soils and aquifers, Water Resour. Res., 25(6), 1413 1421. Benekos, I. D., Shoemaker, C. A., and Stedinger, J. R. (2006). Probabilistic risk and uncertainty analysis for bioremediation of four chlorinated ethenes in groundwater, Stochastic Environmental Research and Risk Assessment, 21(4), 375-390. doi:10.1007/s00477-006-0071-4. Benson, D. A. and Meerschaert, M. M. (2008). Simulation of chemical reaction via particle tracking: Diffusion-limited versus thermodynamic rate-limited regimes. Water Resour. Res., 44(12), 1-7. doi:10.1029/2008WR007111. Benson, D. a., and Meerschaert, M. M. (2009). A simple and efficient random walk solution of multi-rate mobile/immobile mass transport equations, Advances in Water Resources, 32(4), 532-539. doi:10.1016/j.advwatres.2009.01.002. Boso, F., A. Bellin, M. Dumbser (2013), Numerical simulations of solute transport in highly heterogeneous formations: A comparison of alternative numerical schemes, Advances in Water Resources, 52, 178-189. Bouwer, E. J., B. E. Rittmann, and P. L. McCarty (1981), Anaerobic degradation of halogenated 1- and 2-carbon organic compounds, Environ. Sci. Technol., 15(5), 596599. 34 Bundt, M., Widmer, F., Pesaro, M., Zeyer, J., Blaser, P. (2001), Preferential flow paths: biological hot spots in soils, Soil Biol. Biochem., 33, 729-738. Burnell, D. K., J. W. Mercer, and C. R. Faust (2014), Stochastic modeling analysis of sequential first-order degradation reactions and non-Fickian transport in steady state plumes, Water Resour. Res., 50, doi:10.1002/ 2013WR013814. Carrera, J., S´anchez-vila, X., Benet, I., Medina, A., Galarza, G., and Guimer , J. (1998). On matrix diffusion : formulations, solution methods and qualitative effects. Hydrogeology Journal, 178–190. Chapelle, F. (2001), Groundwater microbiology and geochemistry, Wiley, 477 pp. Clement, T. P. (1997), Technical Report PNNL-SA-11720: RT3D - A modular computer code for simulating reactive multispecies transport in 3- Dimensional groundwater systems. Version 1.0. Pacific Northwest National Laboratory, Richland. Clement, T. P. (2001), Generalized solution to multispecies transport equations coupled with a first-order reaction network, Water Resour. Res., 37, 157-163. Cunningham, J. A., C. J. Werth, M. Reinhard, and P. V. Roberts (1997), Effects of grain-scale mass transfer on the transport of volatile organics through sediments: 1. Model development, Water Resour. Res., 33(12), 27132726, doi:10.1029/97WR02425. Cunningham, J. A., and I. Mendoza-Sanchez (2006), Equivalence of two models for biodegradation during contaminant transport in groundwater, Water Resour. Res., 42, W02416, doi:10.1029/2005WR004205. Cunningham J. A. and Z. J. Fadel (2007), Contaminant degradation in physically and chemically heterogeneous aquifers,J. Cont. Hydrology, 94 (3-4), 293-304. De Barros, F. P. J., Rubin, Y., and Maxwell, R. M. (2009). The concept of comparative information yield curves and its application to risk-based site characterization. Water Resources Research, 45(6), 1–16. doi:10.1029/2008WR007324. 35 De Barros, F. P. J., Fern`andez-Garcia, D., Bolster, D., and Sanchez-Vila, X. (2013). A risk-based probabilistic framework to estimate the endpoint of remediation: Concentration rebound by rate-limited mass transfer. Water Resour. Res., 49(4), 1929–1942. doi:10.1002/wrcr.20171. Delay F. and J. Bodin, Time domain random walk method to simulate transport by advection-dispersion and matrix diffusion in fracture networks, Geophys. Res. Lett. 28 (21) (2001) 4051. Dentz M. and B. Berkowitz, Transport behavior of a passive solute in continuous time random walks and multirate mass transfer, Water Resour. Res. 39 (5) (2003) 1111–1131. Ding, D., Benson, D. a., Paster, A. and Bolster, D. (2013). Modeling bimolecular reactions and transport in porous media via particle tracking. Advances in Water Resources, 53, 56–65. doi:10.1016/j.advwatres.2012.11.001. Donado, L., X. Sanchez-Vila, M. Dentz, J. Carrera, and D. Bolster (2009), Multicomponent reactive transport in multicontinuum media, Water Resour. Res., 45,W11402, doi:10.1029/2008WR006823. Edery, Y., H. Scher, and B. Berkowitz (2009), Modeling bimolecular reactions and transport in porous media, Geophys. Res. Lett.,36, L02407, doi:10.1029/2008GL036381. Edery, Y., H. Scher, and B. Berkowitz (2010), Particle tracking model of bimolecular reactive transport in porous media, Water Resour. Res., 46, W07524, doi:10.1029/2009WR009017. Environmantal Protection Agency (EPA) (1999). Anaerobic biodegradation rates of organic chemicals on groundwater: A summary of field and laboratory studies. Rep. EPA S0535, Washington D.C. Environmantal Protection Agency (EPA) (2002). Calculation and use of firstorder rate constants for monitored natural attenuation studies - Ground water issue. Rep. EPA/540/S-02/500, Washington D.C. Falta R. W., M. B. Stacy, A. N. M. Ahsanuzzaman, M. Wang, and R. C. Earle (2007), REMChlor user’s manual, Remediation evaluation model for chlorinated solvents, Center for Subsurface Modeling Support Ground Water 36 and Ecosystems Restoration Division, US EPA R.S. Kerr Environmental Research Center, Ada, Oklahoma. Feehley, C. E., C. Zheng, and F. J. Molz (2000), A dual-domain mass transfer approach for modeling solute transport in heterogeneous aquifers: Application to the Macrodisperison Experiment (MADE) site, Water Resources Research,36 (9), 2501–2515. Fennell, D.E., A.B. Carroll, J.M. Gossett, S.H. Zinder (2001), Assessment of indigenous reductive dechlorinating potential at a TCE-contaminated site using microcosms, polymerase chain reaction analysis, and site data, Environ. Sci. Technol., 35 (9), pp. 1830-1839. Fern`andez-Garcia, G. Llerar-Meza, Gmez-Hernndez, J. J. (2009), Upscaling transport with mass transfer models: Mean behavior and propagation of uncertainty, Water Resour. Res., 45, W10411, doi:10.1029/2009WR007764. Fern`andez-Garcia D., T. H. Illangasekare and H. Rajaram (2005), Differences in the scale-dependence of dispersivity estimated from temporal and spatial moments in chemically and physically heterogeneous porous media, Adv. Water Resour. Res, 28(7), 745-759. Fern`andez-Garcia D. and X. Sanchez-Vila (2011), Optimal reconstruction of concentrations, gradients and reaction rates from particle distributions, J. of Cont. Hydro., 120-121, 99-114. Haggerty, R., and Gorelick, S. M. (1995). Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity, Water Resour. Res., 31(10), 2383-2400. Haggerty, R., Mckenna, S. A., and Meigs, L. C. (2000). On the late-time behavior of tracer test breakthrough curves. Water Resour. Res., 36(12), 3467–3479. Haggerty, R. (2009). STAMMT-L: Solute Transport and Multirate Mass Transfer, version 3.0, user’s manual, ERMS549160. Haggerty, R., C. F. Harvey, C. Freiherr von Schwerin, and L. C. Meigs (2004), What controls the apparent timescale of solute mass transfer in aquifers 37 and soils? A comparison of experimental results, Water Resources Research,40, W01510, doi:10.1029/2002WR001716. Harbaugh A., E. Benta, M. Hill and M. McDonald (2000),MODFLOW 2000 the US Geological Survey Modular groundwater model-user guide to modularization concepts and the ground-water flow process, Open File Rep. U.S. Geol. Surv., 00–92, 121pp. Harvey, C., and Gorelick, S. M. (2000). Rate-limited mass transfer or macrodispersion : Which dominates plume evolution at the Macrodispersion Experiment (MADE) site? Water Resour. Res., 36(3), 637–650. Haston, Z. C., and P. L. McCarty (1999), Chlorinated ethene half-velocity coefficients (Ks) for Reductive dehalogenation, Environ. Sci. Technol., 33, 223226. Henri, C. V., and D. Fern`andez-Garcia (2014), Toward efficiency in heterogeneous multispecies reactive transport modeling: A particle-tracking solution for first-order network reactions, Water Resour. Res., 50, 72067230, doi:10.1002/2013WR014956. Henri, C. V., Fernandez-Garcia, D., de Barros, F. P. J., submitted. Probabilistic human health risk assessment of chemical mixtures in heterogeneous aquifers: Risk statistics, hot spots and preferential flow channels. Water Resour. Res. Huang, H., and B. X. Hu (2000), Nonlocal nonreactive transport in heterogeneous porous media with interregional mass diffusion, Water Resour. Res., 36(7), 1665–1675. Kitanidis, P. K. (1994), Particle-tracking equations for the solution of the advection-dispersion equation with variable coefficients, Water Resour. Res., 30(11), 3225–3227. Lichtner, P. C.,S. Kelkar, and B. Robinson (2002), New form of dispersion tensor for axisymmetric porous media with implementation in particle tracking, Water Resources Research,38, 1146, doi:10.1029/2000WR000100. Li, Z., and M. L. Brusseau (2000), Nonideal transport of reactive solutes in heterogeneous porous media: 6. Microscopic and macroscopic approaches 38 for incorporating heterogeneous rate-limited mass transfer, Water Resour. Res., 36(10), 2853–2867. Liu G., C. Zheng and S. Gorelick, Limits of applicability of the advectiondispersion model in aquifers containing connected high-conductivity channels., Water Resources Research 40 (2004) W08308. Llopis-Albert, C. , and Capilla, J. E. (2009). Gradual conditioning of non- Gaussian transmissivity fields to flow and mass transport data. 3. Application to the macrodispersion experiment (MADE-2) site, on Columbus Air Force Base in Mississippi (USA). J. Hydrol., 371 (14 ), 7584. Mallawatantri, A.P., McConkey, B.G., Mulla, D.J. (1996), Characterization of pesticide sorption and degradation in macropore linings and soil horizons of Thatuna silt loam. J. Environ. Qual., 25, 227-235. Maxwell, R. M., Carle, S. F., and Tompson, A. F. B. (2008). Contamination, risk, and heterogeneity: on the effectiveness of aquifer remediation. Environmental Geology, 54(8), 1771–1786. doi:10.1007/s00254-007-0955-8 Michalak, A. M., and P. K. Kitanidis (2000), Macroscopic behaviour and random-walk particle tracking of kinetically sorbing solutes, Water Resour. Res., 36(8), 2133–2146. Miralles-Wilhelm, F. and L.W. Gelhar (1996), Stochastic analysis of transport and decay of a solute in heterogeneous aquifers, Water Resour. Res., 32 (12), pp. 3451?3459 Mishra, B. K., and C. Mishra (1991), Kinetics of nitrification and nitrate reduction during leaching of ammonium nitrate through a limed Ultisol profile, J. Indian Soc. Soil Sci., 39, 221228. Moler, C., and Van Loan, C. (2003). Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later. SIAM Review, 45(1), 3–49. doi:10.1137/S00361445024180. Neretnieks, I. (1980), Diffusion in rock matrix: An important factor in radionuclide retardation?, J. Geophys. Res., 85(B8), 4379–4397. Painter, S., V. Cvetkovic, J. Mancillas, and O. Pensado (2008), Time domain particle tracking methods for simulating transport with retention and firstorder transformation, Water Resour. Res.,44, W01406. 39 Paster, A., D. Bolster, and D. Benson (2014), Connecting the dots: Semianalytical and random walk numerical solutions of the diffusion-reaction equation with stochastic initial conditions, Journal of Computational Physics, 263, 91-112. Pivetz, B.E., Steenhuis, T.S. (1995), Soil matrix and macropore biodegradation of 2,4-D. J. Environ. Qual., 24 (4), 564-570. Rao, P. S. C., D. E. Rolston, R. E. Jessup, and J. M. Davidson, Solute transport in aggregated porous media: Theoretical and experimental evaluation, Soil Sci. Soc. Am., 1, 44, 1139–1146, 1980. Rittmann, B. E., and P. L. McCarty (1980), Model of steady state-biofilm kinetics, Biotechnol. Bioeng., 22(11), 2343 2357. Riva, M., A. Guadagnini, D. Fernandez-Garcia, X. Sanchez-Vila and T. Ptak (2008), Relative importance of geostatistical and transport models in describing heavily tailed breakthrough curves at the Lauswiesen site, Journal of Contaminant Hydrology,101, 1–13. Ross, S. M. (2003), Introduction to probability models, 8th ed., 755 pp., Academic Press, Oxford. Rubin, Y. (2003), Applied Stochastic Hydrogeology, Oxford Univ. Press, Oxford. Rubol, S., A. Freixa, A. Carles-Brangari, D. Fern`andez-Garcia c, A.M. Romani b, X. Sanchez-Vila (2014), Connecting bacterial colonization to physical and biochemical changes in a sand box infiltration experiment, Journal of Hydrology, 517, 317–327. Salamon, P., D. Fern`andez-Garcia, J. Gomez-Hernandez, A review and numerical assessment of the randomwalk particle tracking method, Journal of Contaminant Hydrology 87 (1-3) (2006) 277–305. Salamon, P., Fern`andez-Garcia, D., and G´omez-Hern`andez, J. J. (2006). Modeling mass transfer processes using random walk particle tracking, Water Resour. Res., 42(11), 1-14. doi:10.1029/2006WR004927. Salamon, P., D. Fern`andez-Garcia, and J. J Gomez-Hernandez (2007), Modeling tracer transport at the made site: The importance of heterogeneity, Water Resour. Res.,43, W08404. 40 Sandrin, S. K., M.L. Brusseau, J.J. Piatt, A.A. Boudour, W.J. Blanford, N.T. Nelson (2004), Spatial variability of in situ microbial activity: biotracer tests, Ground Water, 42 (3), 374-383. Scholl, M. A. (2000). Effects of heterogeneity in aquifer permeability and biomass on biodegradation rate calculations-Results from numerical simulations. Ground Water, 38, 5, 702–712. Soga, K., J. Page, and T. Illangasekare (2004), A review of NAPL source zone remediation efficiency and the mass flux approach, J. Hazard. Mater., 110(13), 13–27. Stroo, A. Leeson, J. A. Marqusee, P. C. Johnson, C. H. Ward, M. C. Kavanaugh, T. C. Sale, Ch. J. Newell, K. D. Pennell, C. A. Lebr´on, and M. Unger (2012), Chlorinated Ethene Source Remediation: Lessons Learned, Environ. Sci. Technol. 2012, 46, 6438-6447, dx.doi.org/10.1021/es204714w. Sun, Y., J. N. Petersen, T. P. Clement, and R. S. Skeen (1999), Development of analytical solutions for multispecies transport with serial and parallel reactions, Water Resources Research, 35(1), 185–190. Sun, Y., and Buscheck, T. a. (2003). Analytical solutions for reactive transport of N-member radionuclide chains in a single fracture. Journal of Contaminant Hydrology, 62–63, 695–712. doi:10.1016/S0169-7722(02)00181-X Tompson, F. B. (1993). Numerical simulation of chemical migration in physically and chemically heterogeneous porous media, Water Resour. Res., 29(11), 3709–3726. Tsang Y. W. and C. F. Tsang, A particle-tracking method for advective transport in fractures with diffusion into finite matrix blocks, Water Resources Research 37 (3) (2001) 831–835. Valocchi, A., and H. A. M. Quinodoz (1989), Application of the random walk method to simulate the transport of kinetically adsorbing solutes, Groundwater Contamination,185, 35–42. van Genuchten, M. T., and P. J. Wierenga, Mass transfer studies in sorbing porous media, I, Analytical solutions, Soil Sci. Soc. Am. J., 40, 473-480, 1976. 41