Journal of Contaminant Hydrology 157 (2014) 106–116

Contents lists available at ScienceDirect

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

Stochastic analytical modeling of the biodegradation of steady plumes A. Zarlenga, A. Fiori ⁎ Dipartimento di Ingegneria, Universita di Roma Tre, Rome, Italy

a r t i c l e

i n f o

Article history: Received 11 June 2013 Received in revised form 18 November 2013 Accepted 19 November 2013 Available online 2 December 2013 Keywords: Biodegradation Transport in porous media Heterogeneous porous media Bioremediation Stochastic model

a b s t r a c t We present a stochastic analytical framework to assess the contaminant concentration of a steady plume undergoing biodegradation. The method is focused on heterogeneous formations, and it embeds both fringe and core degradation. The Lagrangian concentration approach of Fiori (2001) was employed, which is suited for describing the interplay between the large scale advection caused by heterogeneity and the local dispersion processes. The principal scope of the model is to provide a relatively simple tool for a quick assessment of the contamination level in aquifers, as function of a few relevant, physically based dimensionless parameters. The solution of the analytical model is relatively simple and generalizes previous approaches developed for homogeneous formations. It is found that heterogeneity generally enhances mixing and degradation; in fact, the plume shear and distortion operated by the complex, heterogeneous velocity field facilitates local dispersion in diluting the contaminant and mixing it with the electron acceptor. The decay of the electron donor concentration, and so the plume length, is proportional to the transverse pore-scale dispersivity, which is indeed the parameter ruling mixing and hence degradation. While the theoretical plume length is controlled by the fringe processes, the core degradation may determine a significant decay of concentration along the mean flow direction, thus affecting the length of the plume. The method is applied to the crude oil contamination event at the Bemijdi site, Minnesota (USA). © 2013 Elsevier B.V. All rights reserved.

1. Introduction Many contaminants present in the subsurface can be degraded by biological processes where microorganisms, naturally present in the groundwater, metabolize organic contaminant producing inorganic material (Farhadian et al., 2008). Natural degradation is thus an attractive technique for its effects and costs, although it may involve long times and large areas. In the degradation processes the contaminant acts as the electron donor and carbon source, and degradation occurs when the contaminant mixes with an electron acceptor, such as oxygen or nitrate. Hence, local mixing is a crucial component which triggers the biodegradation of organic compounds in groundwater. In turn, mixing is enhanced by the spatial heterogeneity of the aquifer hydraulic properties. In ⁎ Corresponding author. E-mail address: aldo.fi[email protected] (A. Fiori). 0169-7722/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jconhyd.2013.11.003

fact, solute plumes do stretch and distort because of heterogeneity, developing solute “fingers” with significant lateral concentration gradients, such that transverse local dispersion/ diffusion is more effective in diluting the contaminant and mixing the electron donor and acceptor. Hence, the interplay between local dispersion and large scale advection due to heterogeneity is crucial in determining local mixing and thus the degradation processes. Plumes originating from old contamination events may easily reach steady state conditions; previous analyses dealing with field investigation (e.g. Essaid et al., 1995; Maier and Grathwohl, 2005; Van Breukelen et al., 2004) suggest adopting the steady state transport to model the key features of the transport processes in such environments. Most of the previous studies on the subject have analyzed steady transport in homogeneous formations (e.g. Cirpka and Valocchi, 2007; Domenico, 1987; Gutierrez-Neri et al., 2009; Liedl et al., 2005), and only a limited number of papers have considered reactive

A. Zarlenga, A. Fiori / Journal of Contaminant Hydrology 157 (2014) 106–116

transport in heterogeneous media, and mainly for transient plumes (e.g. Bellin et al., 2011; Dagan and Cvetkovic, 1996; Severino et al., 2012). In the present work we develop analytical solutions for steady transport of a contaminant undergoing biodegradation. We are interested in relatively simple solutions, with a limited number of parameters, to be used as a screening tool for a rapid analysis of the area subject to contamination and for risk assessment studies. Following a customary approach (e.g. Dagan, 1989; Rubin, 2003) we model the hydraulic conductivity K as space random function with given statistical properties. Under steady conditions, transport is ruled by the lateral and vertical spreading, and longitudinal transport can be assumed as deterministic (Cirpka and Valocchi, 2007; Zarlenga and Fiori, 2013a). The microbial reactions are separated in two components, pertaining to fringe and core biodegradation, respectively (e.g. Gutierrez-Neri et al., 2009). The former refers to the aerobic processes (e.g. denitrification and sulfate reduction) occurring at the plume fringe (Davis et al., 1999), while the latter concerns with anaerobic processes (e.g. iron and manganese reduction and methanogenesis; see Essaid et al., 1995) and mainly occurs at the core of the plume. For the fringe processes, the microbial reaction is assumed as instantaneous, similar to the approach of Borden and Bedient (1986) and Ham et al. (2004); the assumption of instantaneous reaction is supported by the numerical investigations presented by Chu et al. (2005), Cirpka and Valocchi (2007), Maier and Grathwohl (2005), which show that the limiting factor for natural attenuation is the vertical and transversal mixing rather than the reaction kinetics. However, the adopted approximations can be justified considering the uncertainty on the site characterization and on the biological processes involved. In turn, a simple first-order kinetics is adopted for the core degradation processes (Domenico, 1987; Gutierrez-Neri et al., 2009). Flow and transport are solved by the First order approximation approach (Dagan, 1989; Fiori, 1996; Rubin, 2003), which is formally valid for low/moderate heterogeneity. Along the stochastic approach, all the flow and transport variables are considered as space random functions. The Lagrangian concentration definition of Fiori (2001) is employed here; its main scope is to filter out the advective components which do not directly contribute to mixing, and thus focus on the local mixing processes which are the main triggers of degradation. The assumptions adopted allow developing relatively simple analytical or semianalytical solutions for steady transport undergoing biodegradation in heterogeneous porous formations. 2. Mathematical framework Steady flow of mean uniform velocity U takes place in a natural, heterogeneous porous formation, with spatially variable hydraulic conductivity K(x), where x(x1,x2,x3) is the reference system aligned with the mean flow direction. The logconductivity Y = lnK is modeled as a stationary spatial random variable, normally distributed, with mean and variance equal to 〈Y〉 and σ2Y, respectively. The spatial correlation is described by the two-point covariance CY(r), which depends on the vector distance r between two points. As suggested by field evidence (see e.g. Table 2.3 in Rubin, 2003), CY is assumed as axisymmetric, with integral scales equal to (IY,h,IY,h,IY,v) along

107

the directions (x1,x2,x3); thus, the mean velocity U is aligned with the principal anisotropy direction IY,h. With e = IY,v/IY,h the anisotropy ratio, it is usually e = O(10−1) (Rubin, 2003). The velocity field v in the medium is random and spatially variable, with mean equal to 〈v〉 = (U,0,0). A contaminant (e.g. an organic compound which serves as an electron donor) is continuously released over the area A0 = L2L3 placed at x1 = 0, with L2L3 the transverse and vertical dimensions of the plume, respectively. The contaminant concentration CD over the injection area is constant and equal to C0D. Solute undergoes large scale macrodispersion and pore-scale dispersion (hereinafter PSD), and in the large time limit the concentration assumes a steady profile. While for transient plumes the longitudinal dispersion–diffusion plays a crucial role, steady transport is mainly governed by the lateral and vertical macro- and pore-scale dispersion (e.g. Wexler, 1992; Zarlenga and Fiori, 2013b). After the contaminant release, complex biological processes occur; if all the substrates involved are available, natural communities of microorganisms (bacteria and fungi) mineralize organic C, N, P, S and CO2 or other inorganic compounds which are released in the groundwater. The previous reactions constitute the natural attenuation process or bioremediation. Microorganism present in the subsurface are highly versatile, and biodegradation occurs in many different field conditions. We schematize the complex reactions in two categories, whether degradation occurs under aerobic or anaerobic conditions, respectively. This assumptions is rather common, (e.g. Essaid et al., 2011; Gutierrez-Neri et al., 2009; Suthersan and Payne, 2005) and allows to grasp the main features occurring in the natural attenuation processes. Fringe reactions account for the contaminant degradation in aerobic conditions. Those reactions typically take place on the plume fringe where the contaminant, acting as electron donor (subscript D), is mixed with the electron acceptor (subscript A) provided by the fresh, uncontaminated water. Because of the relatively fast kinetics (Chu et al., 2005; Davis et al., 1999) the pore scale dispersion plays a crucial role in controlling the reactions. The latter are summarized through a global bimolecular irreversible reaction f D D þ f A A→f P P

ð1Þ

where P is the product and fi is the stoichiometric coefficient of the i-th specie. Then, core reactions account for the reactions occurring in anoxic conditions, typically in the core plume region (e.g. Fe, Mn reduction and methanogenesis). Although the core reactions are often neglected in many degradation models, recent studies recognize the primary importance of those mechanism on the plume degradation and the need to introduce them, in combination with the fringe model (Amos et al., 2011; Essaid et al., 2011; Lovley et al., 1989). Those reactions are usually represented by the exponential model dC D ¼ −λC D dt

ð2Þ

Each contaminant species is characterized by its own decay coefficient that depends on the particular field condition; for simplicity a global single effective coefficient λ is adopted here.

108

A. Zarlenga, A. Fiori / Journal of Contaminant Hydrology 157 (2014) 106–116

With the above definitions, we present in the following the basic equations describing the reactive transport. The analysis is performed for unsteady conditions, steady transport being a particular case. Assuming a constant porosity, transport at the Darcy scale is governed by the following advection–dispersion-reaction equations

′

ð3Þ

∂C D  2  λt þ v  ∇C D ¼ Dd  ∇ C D −re ∂t ∂C A  2  þ v  ∇C A ¼ Dd  ∇ C A −r ∂t

ð4Þ

Then, we subtract the previous equations and introduce the dimensionless variables t′ = tU/IY,h, x′i = xi/IY,h, v′i = vi/U and the dimensionless concentration difference 



with δ ¼ C D −C A

ð5Þ

obtaining the final equation ∂R D ′ ′ ′2 þ v  ∇ R ¼ d  ∇ R− UIY;h ∂t ′

Dd 2  ∇′ R UIY;h

ð7Þ

where R = 1 over the area A0 and R = 0 outside A0 at t = 0. We emphasize that the variable R matches the boundary conditions of the variables CD(x) and CA(x), i.e. 0

where Dd = (DdL,DdT,DdT) is the pore-scale dispersion tensor, v = (v1,v2,v3) is the local velocity, Ci is the i-th specie concentration and r is the reaction rate accounting for the fringe degradation processes (Eq. (1)). We remind that, because of the randomness of Y, the velocity v and the reaction rate r (which depends on the random concentrations CD,CA) are both random functions. The boundary conditions are of constant electron donor concentration CD = C0D at the injection plane x ∈ A0, constant electron acceptor concentration CA = C0A and zero concentration gradients at infinity. The above formulation can be conveniently simplified by trying to filter out (i) the reaction rate r from the above equations, along the procedure of Rubin (Rubin, 1983), which was also implemented by Cirpka and Valocchi (2007) and Bellin et al. (2011) among others, and (ii) the decay term expressed by λ. For the latter, we substitute CD = fDC∗Dexp(− λt), CA = fAC∗A in the first two equations of Eq. (3), arriving at

R¼

′

v ∇R¼

∂C D 2 þ v  ∇C D ¼ Dd  ∇ C D −f D r−λC D ∂t ∂C A 2 þ v  ∇C A ¼ Dd  ∇ C A − f A r ∂t ∂C P 2 þ v  ∇C P ¼ Dd  ∇ C P þ f P r ∂t

0 δ þ C A =f A 0 C D =f D þ C 0A =f A

Eq. (6) can be neglected, leading to the classic advection dispersion equation for a conservative solute, expressed in dimensionless form; as a consequence, the final equation for steady transport (which is the aim of the work) writes as follows


′  rI Y;h =U eλt IY;h =U −1 C 0D =f D þ C 0A =f A

C D ¼ C D ; C A ¼ 0 → R ¼ 1 ðx∈A0 Þ 0 C D ¼ 0 ; C A ¼ C A → R ¼ 0 ðx→∞Þ

Summarizing, the concentration CD of a solute undergoing fringe and core biodegradation can be assessed from the solution of the ADE (7) for a conservative solute R, the solution of which is much simpler. Despite of the above simplifications, the solution of the ADE in presence of heterogeneous velocity fields is still a challenging task. In order to further simplify the problem and obtain simple analytical solutions, we adopt the findings obtained by previous work and neglect both macrodispersion and PSD in the longitudinal direction, as they do not contribute much to steady transport (see previous discussion and Zarlenga and Fiori, 2013b). This way, one can substitute t = x1/U (i.e. the travel time of the generic solute particle to reach the control plane at distance x1) in the term CD = fDC∗D exp(− λt), leading to C∗D = CDexp(λx1/U)/fD. This point shall be retaken in the next Section. In order to obtain CD from R an analytical expression for the speciation is needed. For simplicity we adopt an instantaneous and irreversible reaction (Borden and Bedient, 1986; Chu et al., 2005; Cirpka and Valocchi, 2007; Gutierrez-Neri et al., 2009; Liedl et al., 2005), for which the electron donor and the electron acceptor cannot coexist at the same time, being CDCA ¼ 0

ð9Þ

Because of the relation between R and CD (Eq. (5)) and the above assumptions, CD is easily obtained by solving the following system (

λx1 =U

CD e ð6Þ

The above is the final equation which describes the fate of the dimensionless concentration difference R. It is seen that the filtering procedure does not completely eliminate the reaction term r, which is still present in the transport equation for the concentration difference. The reaction term vanishes if λ = 0, i.e. in absence of core degradation. Thus, the above procedure suggests that it is not generally possible to solve the fringe and core degradation processes by combining the correspondent solutions for conservative solutes. However, the last term of Eq. (6) is rather small for processes in which the core degradation process is not strong, i.e. when λIY,h/U ≪ 1. In such cases, the last term of

ð8Þ


 0 0 0 =f D −C A =f A ¼ R C D =f D þ C A =f A −C A =f A CDCA ¼ 0

ð10Þ

whose solutions are C D ðxÞ ¼ C 0D



−λx1 =U

e

½Rð1 þ θÞ−θ 0



for RNR  for R≤R

8  0  for RNR C A ðxÞ < 1  ¼ þ1 for R≤ R : 1−R C 0A θ

ð11Þ

ð12Þ

where θ¼

C 0A f D C 0D f A

;



R ¼

f D C 0A ¼ þ f D C 0A

f A C 0D

  1 −1 1þ θ

ð13Þ

A. Zarlenga, A. Fiori / Journal of Contaminant Hydrology 157 (2014) 106–116

are the fringe-reaction dimensionless parameter and the nonreactive concentration R⁎ associated with CD = 0, respectively. The above relations allow determining the CD field, which is the objective of the present work. It depends on the relevant parameter θ, which is the ratio between the background concentration of the electron acceptor and the concentration of the electron donor at the source; basically, biodegradation is more effective for increasing values of θ. The concentration CD depends on the particular solution for the random dimensionless nonreactive concentration variable R, solution of Eq. (7), which leads to CD after substitution of R in Eq. (11). Thus, the crux of the matter is the calculation of R, pertaining to nonreactive transport, as function of the spatially varying velocity field. The solution for R is derived in the following section. 3. Stochastic Lagrangian analysis of nonreactive transport Eq. (7) describes the dimensionless steady concentration R of a nonreactive solute in a heterogeneous formation. Recently, Zarlenga and Fiori (2013b) developed a stochastic approach to analyze the concentration field of a nonreactive solute under steady transport; the solution was derived under the Firstorder approximation, which is valid for low to moderate heterogeneity (Dagan, 1989; Rubin, 2003). As shown there, the local concentration statistics depend on the combined effects of large scale heterogeneity and PSD, as well known from the existing literature on transient transport in heterogeneous formations (Fiori and Dagan, 2000; Kapoor and Gelhar, 1994), (Cirpka et al., 2011; Pannone and Kitanidis, 1999). The interplay between PSD and macroadvection is crucial in order to model mixing and predict the local concentration statistics. In order to better focus on the local mixing processes, Fiori (2001) introduced the definition of Lagrangian concentration. The method attempts at filtering out the large scale advection which does not contribute directly to local mixing. The approach is formalized and solved analytically for unsteady transport by Fiori (2001); it was applied to field experiments, and more recently it was successfully tested against alternative methods for risk assessment problems (Boso et al., 2013). The method displays similarities with the streamline approach introduced by Cirpka et al. (2011, 2012), and it is particularly suited for the prediction of the higher concentrations, although Boso et al. (2013) showed that it can work quite well for the entire range of concentrations in the plume. A brief recapitulation of the solution and its extension to steady transport is provided below, while a more detailed and general description of the approach is found in (Fiori, 2001). The plume is represented as a collection of subparticles defined at the pore scale. The total trajectory Xt of the subparticle originated from x = a (a ∈ A0) is the sum of the advective (X) and diffusive (Xd) terms, i.e. Xt = X + Xd. The basic idea of the Lagrangian concentration approach is to describe the transport process as observed from the advective trajectory of a Darcy-particle P(t;c) = lim Δ → 0Δ−1 ∫ Δ X(t; c + a)da (i.e. of vanishing support volume Δ) originated from x = c, i.e. by switching to a new reference system ξ(x1,x2 − P2(t;c),x3 − P3(t;c)) attached to P (see Fiori, 2001). This way, the dilution process is followed closely for each particle, filtering out part of

109

the spatial variability of trajectories due to macroadvection. The Lagrangian concentration R pertaining to the subparticles for a continuous injection is derived along the developments of Zarlenga and Fiori (2013b; Eq. (1)), with the above change of coordinates. Since we neglect longitudinal macrodispersion (see previous section) the travel time τ of a generic particle is equal to τ = x1/U, and the latter is substituted in Eq. (1) of Zarlenga and Fiori (2013b)leading to Z R¼

A0

h i h i Hðt−x1 =U Þδ ξ2 −W t;2 ðx1 =U; a; cÞ δ ξ3 −W t;3 ðx1 =U; a; cÞ da

ð14Þ where Wt = Xt(t;a) − P(t;c). The Darcy-scale Lagrangian Concentration is obtained after averaging the above expression over the PSD term Xd. For the sake of simplicity we perform this integration together with the integration over the trajectories, which is required in order to estimate the concentration moments (Fiori and Dagan, 2000). Thus, the mean concentration is obtained by averaging Eq. (14) over the displacements Wt and taking the limit t → ∞ Z hRi ¼

A0

f W t;2 ðξ2 ; x1 =U; a; cÞf W t;3 ðξ3 ; x1 =U; a; cÞda

ð15Þ

with f W t;i the pdf of Wt,i (i = 2,3). Fiori (2001)showed that at first order the f W t;i pdfs are Gaussian, with the first and second statistical moments equal to hWt ðx1 ; a−cÞi ¼ a−c W t;ii ðx1 ; a−cÞ ¼ X ii ðx1 =U Þ þ Z ii ðx1 =U; 0Þ þ2α dT x1 −2Z ii ðx1 =U; a−cÞ

ð16Þ

where Xii(t) and Zii(t;r) are the variance and the covariance of the i-th component of the trajectories of two particles injected at a distance r (Dagan, 1989); αdT = DdT/U is the transverse PSD dispersivity. A great simplification of the problem can be obtained assuming in Eq. (16) a − c ≈ 0, which is formally valid for small solute injection areas (L2,L3 ≪ IY,h); this assumption was introduced by Fiori and Dagan (2000) and tested by Tonina and Bellin (2008) for larger injections areas, providing a useful analytical simplification of the trajectory moments (Eq. (16)), that now become hWt ðx1 ; a−cÞi≃0 W t;ii ðx1 ; a−cÞ≃X ii ðx1 =U Þ þ 2α dT x−Z ii ðx1 =U; 0Þ ¼ X t;ii ðx1 =U Þ−Z ii ðx1 =U; 0Þ

ð17Þ

It can be shown (Fiori, 2001) that under the above assumption the choice of c is immaterial, and in the following we shall assume c = 0. Similarly, the variance σ2R = 〈R2〉 − 〈R〉2 can be calculated from Eq. (14). Squaring the latter and taking the limit t → ∞ it yields for 〈R2〉 D E Z 2 R ¼

Z A0

A0

gW t;2 ðξ2 ; ξ2 ; x1 =U; a; bÞgW t;3 ðξ3 ; ξ3 ; x1 =U; a; bÞdadb

ð18Þ where g W t;2 is the joint pdf of two relative displacements Wt,2(x1/U;a,0), Wt,2(x1/U;b,0), and similarly for the vertical component Wt,3; under the First order approach employed here, g W t;2 and g W t;3 are bivariate normal. The joint moment

110

A. Zarlenga, A. Fiori / Journal of Contaminant Hydrology 157 (2014) 106–116

between Wt,2(x1/U;a,0) and Wt,2(x1/U;b,0) is easily calculated as σ W t2 ¼ Z 22 ðx1 =U; 0Þ þ Z 22 ðx1 =U; a−bÞ−Z 22 ðx1 =U; aÞ−Z 22 ðx1 =U; bÞ and under the aforementioned assumption of small plumes (a ≈ b) it becomes σ W t;2 ¼ 0. As a consequence, we can write gW t;2 ðξ2 ; ξ2 ; x1 =U; a; bÞ ¼ f W t;2 ðξ2 ; x1 =U; a; 0Þf W t;2 ðξ2 ; x1 =U; b; 0Þ and similar for the vertical component. Substitution of the latter in Eq. (18) gives to 〈R2〉 = 〈R〉2, leading to σ2R = 0. The procedure is similar to the one developed for transient plumes by Fiori (2001; Sections 2.3– 2.4), and more details are found there. Thus, the variance of R is null and the concentration field R can be treated as a deterministic quantity. We remark again that this result is only valid for small plumes along the Lagrangian concentration approach, i.e. adopting a moving reference system along a Darcy-scale particle. This important result comes from the filtering out of the purely advective fluctuations of the particle trajectory that do not directly contribute to dilution. As pointed out in Fiori (2001), the concentration R calculated by this method basically represents the concentration pertaining to clusters of particles (or “hot spots”) of the solute after the plume shredding and distortion operated by the random velocity field, regardless of the precise transverse location of the high concentration spots. This way, the R field is very much focused on the high concentration values which emerge in the plume after the combined effects of macrodispersion and local dilution. For a constant initial concentration, the final equation for the dimensionless concentration R of a non reactive solute is easily obtained after integration of Eq. (14), adopting a normal distribution for Wt, with moments given by Eq. (17), the final result being R ¼ hRi ¼ ηðξ2 Þηðξ3 Þ

4. Length of the plume The plume length LP is defined as the distance from the injection plane where the electron donor concentration CD is null. Generally speaking, CD is a random variable, leading to a random LP. Recently, Zarlenga and Fiori (2013a) have taken into account the randomness of CD and determined the plume length by a stochastic procedure, leading to a probability density function for LP. Instead, in the present approach we make use of the Lagrangian concentration approach, which leads to a deterministic CD (see previous Section). Hence, the plume length is easily calculated by seeking the distance from the injection plane such that CD = 0. Employing expression (11), the solution of CD = 0 leads to the following implicit solution for LP

RðLP ; 0; 0Þ ¼

f D C 0A ¼ f A C 0D þ f D C 0A

  1 −1 1þ θ

ð22Þ

ð19Þ

with 8 0 0 1 19 > > = 1< B 2ξi þ Li 2ξi −Li C B C erf @ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiA−erf @ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiA ηðξi Þ ¼ > 2: ; 2 2W ðx =U Þ 2 2W ðx =U Þ > t;ii

1

t;ii

1

ð20Þ where Wt,ii(x1/U) = Wt,ii(x1/U,0) for simplicity. An even simpler result can be obtained in the particular case of point-source injection, which is obtained after expansion of Eq. (19) around L2,L3 ≈ 0, obtaining R¼

closed-form expression; the analytic procedure for their calculation is reproduced in the Appendix A. As a final note toward applications, it is well known that the moments Xt,22, Xt,33 calculated by the first order solutions reproduced in the Appendix A may easily underestimate the transverse and vertical dispersion. If the transverse and vertical dispersivities (sum of macrodispersivity and PSD) At,ii are available (e.g. through field experiments or numerical simulations), the above problem can be overcome by following the approach suggested by Fiori (2001) and recently implemented in (Boso et al., 2013), that is to perform the substitution Wt,ii = 2At,iix1(1 − Zii/Xt,ii) in Eq. (19).

" # L2 L3 ξ22 ξ23 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi exp − − 2W t;22 ðx1 =U Þ 2W t;33 ðx1 =U Þ 2π W t;22 ðx1 =U ÞW t;33 ðx1 =U Þ
 ð21Þ L2 =IY;h ; L3 =IY;v ≪1

It is seen that the final expressions (19) and (21) are rather simple and they are equal to the formulas for homogeneous formations, with Wt,ii substituting 2Dd,iit. The moments Wt,iis correspond to the non-ergodic moments discussed in Dagan (1989), and the present approach shows that the above substitution can be formally made only (i) in a reference system moving with the center of mass of a Darcy particle (Lagrangian concentration approach), and (ii) for small plumes. Although the final expressions (19) and (21) are rather simple, the trajectories covariances Xt,ii, Zt,ii from which the moments Wt,ii are calculated do not have a

The above can be easily solved through Eq. (19). Thus, the plume length does not depend on the core reactions parameter λ. However, the parameter λ in Eq. (11) may determine a significant decay of CD far from the injection plane, with concentrations close to zero at distances much shorter than the one predicted by Eq. (22), leading to a smaller length of the plume. Hence, a more useful definition of LP would be the one based on a low threshold C∗D N 0 (e.g. based on the lower threshold concentration measured by the sampling device) for the electron donor concentration, i.e. adopting the condition CD = C∗D; the previous definition, which led to Eq. (22), is recovered for C∗D = 0. This point shall be retaken and illustrated in the next section. 5. Illustration examples and discussion The solution of the Lagrangian concentration field of a steady plume undergoing biodegradation herein proposed is here exemplified through a series of illustration examples. The concentration CD is given by Eq. (11), with R provided by Eq. (19) in which the trajectory moments are given by Eqs. (A.4), (A.8). The (dimensionless) parameters of the problem are the size of the injection area L2′ = L2/IY,h, L3′ = L3/IY,h, the Peclet number Pe = UIY,h/DdT, the anisotropy ratio e = IY,v/IY,h, the logconductivity variance σ2Y and the two dimensionless reaction parameters θ = C0AfD/(C0DfA) and λ′ = λIY,h/U pertaining to the fringe and core biodegradation processes, respectively.

A. Zarlenga, A. Fiori / Journal of Contaminant Hydrology 157 (2014) 106–116

Fig. 1 depicts the dimensionless electron donor concentration CD/C0D as function of the dimensionless distance x1/IY,h from the source, for the case ξ = 0; the latter pertains to the maximal concentration along the transverse plane at x1. The analysis is for a small plume of size L2 = L3 = 0.1IY,h, σ2Y = 1, e = 0.1 and λ = 0 (i.e. in absence of core biodegradation). We explore in the figure the concentration behavior as function of the Peclet number Pe = IY,h/αdT and θ. It is seen that for Pe = ∞ the electron donor concentration stays constant and equal to its initial value C0D. In fact, no local mixing is possible for Pe = ∞ and the electron donor and acceptor do not mix, and thus biodegradation cannot take place. In such a case the length of the plume is infinite, despite of the nonzero macrodispersion coefficients. We remark that this result was obtained by adopting the Lagrangian concentration formulation (Section 3), while a standard Eulerian concentration model (i.e. CD calculated at fixed locations in space) would have predicted a decay of concentration with distance, which was otherwise determined by the uncertainty in the particle locations and not by a true mixing process. The issue will be retaken later when discussing the application to a field case. For finite Peclet, the decay of concentration with distance is inverse proportional to Pe, the latter triggering the mixing between acceptor and donor, and thus biodegradation. The length of the plume (which shall be discussed later) also strongly depends on Peclet. The ratio θ between the initial and the background concentrations of the electron donor and the electron acceptor, respectively, also displays a significant impact on CD, as visible in Fig. 1. As expected, a higher θ ratio enhances biodegradation, leading to a stronger decay of both CD and the plume length. Obviously, the parameter θ has no influence on the concentration when Pe = ∞ as the electron donor and acceptor cannot mix. The effect of the core biodegradation parameter λ is analyzed in Fig. 2, which has similar parameters as Fig. 1, with a constant Pe = 1000; the concentration CD is depicted as function of distance x1 for a few values of λ′ and θ. As discussed in Section 3, the effect of λ manifests in a rather simple exponential decay applied to the solution for the

111

fringe degradation, as dictated by Eq. (11). Thus, the impact of core biodegradation under such model, as illustrated in Fig. 2, is rather expected and straightforward. The issue that we want to discuss here is the length of the plume, along the lines of Section 4. Fig. 2 clearly shows that, although CD is null at a distance which does not depend on λ′, its actual value can become very small because of the exponential decay. As a consequence, the plume length is for all practical purposes smaller than its theoretical value based on Eq. (22), and the plume length is effectively dependent on λ. This is clearly visible in Fig. 2, where for the case θ = 0.05 the plume ends theoretically (i.e. when CD = 0) at x1/IY,h ≅ 10, while for the largest λ′ = 0.5 the actual length is around x1/IY,h ≅ 4. Hence, a plume length definition based on a given threshold for CD, instead of the CD = 0 condition, would be more appropriate, as discussed in Section 4. With the above comment in mind we examine now the “theoretical” length of the plume LP, which does not depend on λ, as previously discussed. We explore first the impact of Peclet and θ on LP; the comparison is given in Fig. 3, the other parameters being the same as the previous figures. As expected, the length of the plume decreases with θ, i.e. when the background concentration of the electron acceptor (e.g. oxygen) increases, for a given electron donor concentration C0D. Conversely, LP increases with Peclet, and it grows unbounded for Pe → ∞ along a power-law like behavior, with exponent roughly equal to 0.9 for all θ cases. Once again, the increase of the plume length with Peclet is a consequence of the decreasing mixing effects operated by local dispersion/ diffusion. The above findings are similar to those obtained in the past for homogeneous formation. The effect of the spatial heterogeneity is generally to increase the local mixing operated by PSD, which is indeed facilitated by the spreading and fingering of the plume. Such process impacts the local concentration and as a consequence the plume length. Fig. 4 depicts the latter as function of the logconductivity variance σ2Y, which epitomizes the level of heterogeneity in the system, for a few values of θ. The figure shows that the plume length generally decreases with increasing levels of heterogeneity, as discussed above; the decrease is stronger

Fig. 1. The contaminant Lagrangian concentration 〈CD〉 evaluated by the Eqs. (11) and (19) as function of the distance x1 (ξ = 0) for σ2Y = 1, L2 = L3 = 0.1IY,h, e = 0.1 for a few values of the Peclet number (different colors) and for a few fringe reaction parameter θ = C0AfD/(C0DfA) (different kinds of line); core degradation is absent (λ = 0).

112

A. Zarlenga, A. Fiori / Journal of Contaminant Hydrology 157 (2014) 106–116

1 0.8 0.6 0.4 0.2 0 0.01

0.1

1

10

100

Fig. 2. The contaminant Lagrangian concentration 〈CD〉 evaluated by the Eqs. (11) and (19) as function of the distance x1 (ξ = 0) for σY2 = 1, L2 = L3 = 0.1IY,h, e = 0.1, Pe = 1000; few values of the core reaction parameters λ (different colors) are showed for two different core reaction parameters θ = 0.5, 0.05 (solid and dashed line).

for lower values of θ because of the increased mixing induced by heterogeneity is more effective when the electron acceptor concentration is much lower than the donor one. This finding further emphasizes the important role played by heterogeneity, and in particular the interplay between large scale advection and local scale dispersion/diffusion, in determining the concentration field of the electron donor concentration and the length of the plume.

6. Application to a field case We apply here the above methodology to the contamination event at the Bemijdi site, Minnesota (USA). The conditions which led to the contamination and the general setup are rather different from those employed in our model, and our goal here is not to provide an exhaustive analysis of the natural remediation processes which occurred at the site but rather to show an example of the model application and its capabilities.

/

In 1979 a buried pipeline sprayed over an area of 6500 m2 an estimated volume of 11 × 106 l of crude oil. After an enormous cleanup effort, about 0.4 × 106 l of remains trapped in the unsaturated zone and floating on the water table constitutes a continuous source of BTEX. The plume had reached the steady state conditions although some changes in the chemical conditions were documented (Baedecker et al., 1993). The site was object of a huge number of studies investigating different chemical and hydrological aspects of the problems (Cozzarelli et al., 2001; Essaid et al., 2011); few geostatistical analyses were conducted in (Essaid et al., 1993) and (Dillard et al., 1997). The general conclusion was that the bioremediation for this site is principally due to the anaerobic processes and hence core degradation. In the present example we model the concentration of benzene downstream the edge of the contaminated area. The model parameters are derived from the data available in the literature (Dillard et al., 1997; Essaid et al., 1993, 1995; Gutierrez-Neri et al., 2009; Suthersan and Payne, 2005), and they are reproduced in Table 1. The only missing (and mostly

,

1000

100

10

= 0.1 = 0.1

1

0.1 100

1000

10000

Fig. 3. The plume length LP, evaluated by the Eqs. (22) and (19), as function of the Peclet number for a few values of the fringe reaction parameter θ; the other parameters are σ2Y = 1, e = 0.1, L2 = L3 = 0.1IY,h.

A. Zarlenga, A. Fiori / Journal of Contaminant Hydrology 157 (2014) 106–116

113

Fig. 4. The plume length LP, evaluated by the Eqs. (22) and (19), as function of the heterogeneity σ2Y for a few values of the fringe reaction parameter θ; the other parameters are Pe = 1000, e = 0.1, L2 = L3 = 0.1IY,h.

relevant) parameter is the transverse/vertical pore-scale dispersivity, and we have adopted the value αdT = 0.5 mm inferred by Fiori and Dagan (1999) for the Cape Cod and Borden sites, which have similar degree of heterogeneity. Since the planar and vertical extension of the source zone are not well defined, we set the injection area as A0 = L2 × L3 = 80 × 2.5 m2 based on the size of the sprayed area and the vertical profiles of the concentration presented in previous works. The length L3 is consistent with the previous studies water table annual variation (see e.g. Essaid et al., 2011). The moments Wt,ii (see Appendix A), which are introduced in Eq. (19), are reproduced in Fig. 5; the moments Xt,ii are also displayed for comparison. In both cases, the transverse Peclet number Pe = UIY,v/Ddt has been employed in the moments solution. The concentration data was taken from the USGS Bemidji Crude-Oil Research Project (available at http://mn.water. usgs.gov/projects/bemidji/index.html); data was collected by several wells along the mean flow direction, at different vertical levels. We remark that the Lagrangian concentration approach followed here is by definition focused on the highest concentrations along the transverse plane, at a distance x1 from the source, as discussed in Section 3. Thus, our model should predict the higher envelope of the concentration data when ξ = 0. Because of the relatively

Table 1 Model parameters for the Bemijdi site. Water flow mean velocity

U = 0.06 m/day

Hydraulic conductivity field variance Horizontal Integral Scale Vertical Integral scale Anisotropy ratio Pore scale dispersivity Lateral source length Vertical source length Electron Donor concentration Electron Acceptor concentration Stoichiometric ratio Core degradation parameter

σ2Y = 0.87 IY,h = 23 m IY,v = 1.25 m e = 0.054 αdT = 0.0005 m L2 = 80 m L3 = 2.5 m C0D = 15 mg/l C0A = 7.68 mg/l fD/fA = 0.32 λ = 0.0007 1/day

large variability of the measured CD at x1 = 0, the initial concentration C0D was taken as the average between the concentration of the three subsequent years 1992–1993– 1994, i.e. C0D ≃ 15 mg/l. The measured concentrations (dots) and the upper envelope predicted by the present stochastic model (lines; for ξ = 0) are represented in Fig. 6. It is seen that the model provides a reasonable prediction of the high concentrations along the longitudinal distance x1, all the approximations notwithstanding. The figure also provides some sensitivity analysis on the PSD component αdT, and in particular for αdT = 1 mm and 0.25 mm as compared with the assumed one αdT = 0.5 mm. Despite the strong influence exerted by the core degradation (i.e. the parameter λ), the impact of the PSD component αdT (which rules the fringe degradation) is not irrelevant. We emphasize that our solution is formally valid for λIY,h/U ≪ 1, and for the Bemijdi site it is λIY,h/U = 0.268. The reasonably good performance of the model makes us believe that the above condition may not be too stringent. Fig. 6 also displays the maximum concentrations CD predicted by the solution for the homogeneous medium and the one for the heterogeneous medium along the classic Eulerian formulation (e.g. Dagan, 1989) (αdT = 0.5 mm in both cases); the two solutions are obtained after substitution of Wt,ii with 2αdTx1 and Xt,ii, respectively. It is seen that the Eulerian solution considerably overestimates biodegradation, for the reasons previously discussed. Instead, the Homogeneous solution underestimates degradation, and its prediction is closer to the one proposed in the present study, mainly because of the rather important role played by core degradation at the Bemijdi site. 7. Summary and conclusions In this paper we presented an analytical framework to assess the contaminant concentration CD of a steady plume undergoing biodegradation. The method is focused on heterogeneous formations, and it embeds both fringe and core degradation. The Lagrangian concentration approach of Fiori (2001) was employed, which is better suited for describing the

114

A. Zarlenga, A. Fiori / Journal of Contaminant Hydrology 157 (2014) 106–116

25

,

=0.5

20 15 ,

10 5

, ,

0 0

20

40

60

80

100

120

140

160

180

200 (

)

Fig. 5. The trajectory moments Xt,ii and Wt,ii as function of the distance from the injection area; the pore scale dispersivity is αdT = 0.5 mm; the other parameters are those reproduced in Table 1.

• Heterogeneity generally enhance mixing and degradation; in fact, the plume shear and distortion operated by the complex, heterogeneous velocity field facilitates local dispersion in diluting the contaminant and mixing it with the electron acceptor. • Fringe degradation is controlled by the parameter θ = C0AfD/ (C0DfA), an increase of which determines a stronger degradation process and shorter plumes. The effect of θ combines with heterogeneity to further enhance biodegradation. Instead, core degradation is ruled by the dimensionless parameter λIY,h/U, which basically superimpose a linear decay model to the fringe model. • If the condition CD = 0 is adopted, the length of the plume calculated is independent on the core reaction parameter λ. The length of the plume strongly depends on Pe, θ and the logconductivity variance σ2Y. However, the decay of concentration operated by core degradation may be such that the condition CD ≅ 0 may be reached well before the distance predicted by the fringe model, thus leading to a shorter plume length.

interplay between the large scale advection caused by heterogeneity and the local dispersion processes. The principal scope of the model is to provide a relatively simple tool for a quick assessment of the contamination level in aquifers, as function of a few relevant, physically based dimensionless parameters, i.e. the size of the injection area L2′ = L2/IY,h, L3′ = L3/IY,v, the Peclet number Pe = UIY,h/Ddt, the anisotropy ratio of the K structure e = IY,v/Iy,h, the logconductivity variance σ2Y and the two dimensionless reaction parameters θ = C0AfD/(C0DfA) and λ′ = λIY,h/U pertaining to the fringe and core biodegradation processes, respectively. The main conclusions of the work can be listed as follows. • The key role played by mixing in the solute degradation is confirmed and quantified by the analytical model. The latter correctly predicts the lack of degradation when Pe → ∞ for fringe degradation, no matter the degree of heterogeneity. The decay of CD, and so the plume length, is proportional to the transverse pore-scale dispersivity αdT ~ Pe−1, which is indeed the parameter ruling mixing and hence degradation.

1986 1987 1988 1990 1992 1993 1995 1996

20 18 16 14 12

=1 =0.5 =0.25 Eulerian Homogeneous

10 8 6 4 2 0 0

20

40

60

80

100

120

140

160

180

200

Fig. 6. Benzene maximum concentration at the Bemidji site (Minnesota) evaluated along the mean flow direction for different values of the PSD αdT (lines) against field data (1986–1996; dots); the classic Eulerian solution for heterogeneous media (red line) and the solution for homogeneous formations (green line) are also reproduced, both for αdT = 0.5 mm.

A. Zarlenga, A. Fiori / Journal of Contaminant Hydrology 157 (2014) 106–116

We emphasize again that the stochastic analytical model herein developed embeds a few simplifying assumptions regarding the nature of flow, the heterogeneous structure of the hydraulic conductivity as well as the types of reactions involved. Hence the model is particularly suited for a screening tool, e.g. to be employed in risk assessment studies, rather than a complete and detailed method. Nevertheless, the physicallybased nature of the model permits some useful insights into the complex interplay between heterogeneity–local dispersionreactions, which would otherwise require more demanding and involved numerical analyses. Acknowledgment

After few manipulations, recasting the problem in cylinpffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffi drical coordinates k1 ¼ ργ; k2 ¼ 1−γ 2 cosϕ; k3 ¼ 1−γ2 sinϕ and performing three analytical integrations, we obtain, finally, the trajectory covariances X ii ðτ Þ 4ePe ¼ 2 σ 2Y I2Y;h π

In this section we derive the analytical expressions for the displacement covariances Wt,ii under the First-Order assumptions. We start from the well-known equation that relates the ^ mn ðκÞ with the Fourier transforms of the velocity field covariance u ^ Y ðκÞ (Dagan, 1989) Fourier transform of the K-field covariance C ðA:1Þ

where k = (k1, k2, k3) are the Fourier space coordinates and δii is the identity matrix. Recasting the problem with the dimensionless variables ρi = ri/IY,h, τ = Ut/IY,h, ki = κiIY, h, Pe = Dd/UIY,h and assuming an axial symmetric K-field covariance (IY,h, IY,h,eIY,h) and an isotropic pore scale dispersion tensor (Dd11 = Dd22 = Dd33 = Dd) the one particle displacement covariance is

ðA:2Þ

The hydraulic conductivity field covariance CY(r) and the PSD tensor αd should be derived by geostatistical analysis and field tests; we apply the general framework traced in (Fiori, 1996) and developed for an exponential CY in Fiori (1996, 2001) to a Gaussian CY in order to reduce the numerical quadratures required ! π r2 2 C Y ðrÞ ¼ σ Y exp − 2 ðA:3Þ ; 4 IY;h ! 3 2 2 ^ ðkÞ ¼ σ 2 8eIY;h exp − κ IY;h C Y Y π ð2πÞ3=2 where r2 = r21 + r22 + e2r23.

Mii ðρ; γ ÞSX ðρ; γ; τÞdρdγ 0

ðA:4Þ

0

with 
1 
13 2 0
0
1−e2 ρ2 1−γ2 1−e2 ρ2 1−γ2 A−I1 @ A5 2π 2π


1 
13 2 0
0
2 2 2 2 2 2 1−e ρ 1−γ 1−e ρ 1−γ @ A @ A5 M 33 ðρ; γ Þ ¼ γ I0 þ I1 2π 2π 24



 2 0 13 2 2 2 1−γ2 τ γ 1−e þ 1 þ e A5 2 SX ðρ; γ; τ Þ ¼
2 exp4−ρ @ þ 2 2 2 2π Pe ρ þ Pe γ

Appendix A. Trajectory and travel time statistical moments

Zτ Zτ Z h
i X mm ðτÞ 1 ^ mm ðkÞ exp −ik1 τ ′ −τ ″ u ¼ σ 2Y I2Y;h ð2πÞ3=2 0 0 

km kn  ′ ″ ′ ″ exp τ −τ  dkdτ dτ ðm; n ¼ 1; 2; 3Þ Pe

Z∞ Z1

2 M 22 ðρ; γ Þ ¼ γ 4I0 @

Research funding for this study was partially provided by PRIN 2010–2011 “Innovative methods for water resources management under hydro-climatic uncertainty scenarios” (2010JHF437). We thank two anonymous reviewers for their constructive comments, which helped us to improve the manuscript.

  3 X 3
X κ nκ j ^ κ κ ^ mn ðκÞ ¼ C Y ðκ Þ u U i U j δim − m i δjn − κ κ i j

115

"

# 
 i ρ2 τ h
2 2 2 2 2 ρ τ−Pe ρ þ Pe Pe þ ρ τ γ Pe ! h
 i 2 2 2 þPe ρ −Pe γ cosðργτÞ−2Peγ sinðργτ Þ exp

ðA:5Þ

where I0 and I1 are modified Bessel function of the first kind. The trajectory moments of two particles originated at a and b at time t0 = 0 are expressed as function of their travel times τ and τ′
 ′ Z ii τ; τ σ 2Y I2Y;h

¼

1 ð2πÞ3=2

Zτ Zτ Z

 h

^ ii ðkÞ exp −i k1 τ′ −τ″ u

0 0 i þk2 ða2 −b2 Þ þ k3 ða3 −b3 ÞÞ 

k k
′ ″ ′ ″ dkdτ dτ exp m n τ þ τ Pe

ðm; n ¼ 1; 2; 3Þ ðA:6Þ

Neglecting the dependence on a–b and using the cylindrical coordinates, we obtain the final expressions
 Z ii τ; τ′ σ 2Y I2Y;h

¼

2ePe2 π2

Z∞ Z1 0


 ′ Mii ðρ; γÞSZ ρ; γ; τ; τ dρdγ

ðA:7Þ

0

with

 0 1
 1−γ 2 γ 2 1−e2 þ 1 þ e2 ′ 2 A SZ ρ; γ; τ; τ ¼
2 exp@−ρ 2 2 2 2π ρ þ Pe γ ! ! "
 ρ2 τ ρ2 τ ′ ′ cos ργτ cosðργτÞ− exp − 1− exp − Pe Pe
1 0 #

 ρ2 τ þ τ′ A cos ργ τ′ −τ þ exp@− ðA:8Þ Pe

116

A. Zarlenga, A. Fiori / Journal of Contaminant Hydrology 157 (2014) 106–116

References Amos, R.T., Bekins, B.A., Delin, G.N., Cozzarelli, I.M., Blowes, D.W., Kirshtein, J.D., 2011. Methane oxidation in a crude oil contaminated aquifer: delineation of aerobic reactions at the plume fringes. J. Contam. Hydrol. 125, 13–25. http://dx.doi.org/10.1016/j.jconhyd. 2011.04.003. Baedecker, M.J., Cozzarelli, I.M., Eganhouse, R.P., Siegel, D.I., Bennett, P.C., 1993. Crude oil in a shallow sand and gravel aquifer, III, Biogeochemical reactions and mass balance modeling in anoxic groundwater. Appl. Geochem. 8, 6. http://dx.doi.org/10.1016/0883-2927(93)90014-8. Bellin, A., Severino, G., Fiori, A., 2011. On the local concentration probability density function of solutes reacting upon mixing. Water Resour. Res. 47, 1. http://dx.doi.org/10.1029/2010WR009696. Borden, R.C., Bedient, P.B., 1986. Transport of dissolved hydrocarbons influenced by oxygen-limited biodegradation: 1. Theoretical development. Water Resour. Res. 22, 13. http://dx.doi.org/10.1029/ WR022i013p01973. Boso, F., de Barros, F.P.J., Fiori, A., Bellin, A., 2013. Performance analysis of statistical spatial measures for contaminant plume characterization towards risk-based decision making. Water Resour. Res. 49. http:// dx.doi.org/10.1002/wrcr.20270. Chu, M., Kitanidis, P.K., McCarty, P.L., 2005. Modeling microbial reactions at the plume fringe subject to transverse mixing in porous media: when can the rates of microbial reaction be assumed to be instantaneous? Water Resour. Res. 41, 6. http://dx.doi.org/10.1029/2004WR003495. Cirpka, O., Valocchi, A.J., 2007. Two-dimensional concentration distribution for mixing-controlled bioreactive transport in steady state. Adv. Water Resour. 30, 6–7. http://dx.doi.org/10.1016/j.advwatres.2006.05.022. Cirpka, O.A., de Barros, F.P.J., Chiogna, G., Nowak, W., 2011. Probability density function of steady state concentration in two-dimensional heterogeneous porous media. Water Resour. Res. 47, 11. http://dx.doi.org/10.1029/ 2011WR010750. Cirpka, O.A., Rolle, M., Chiogna, G., de Barros, F.P.J., Nowak, W., 2012. Stochastic evaluation of mixing-controlled steady-state plume lengths in two dimensional heterogeneous domains. J. Contam. Hydrol. 138. http://dx.doi.org/10.1016/j.conhyd.2012.05.007. Cozzarelli, I.M., Bekins, B.A., Baedecker, M.J., Aiken, G., Eganhouse, R.P., Tuccillo, M.E., 2001. Progression of natural attenuation processes at a crude oil spill site: I. Geochemical evolution of the plume. J. Contam. Hydrol. 53, 3–4. http://dx.doi.org/10.1016/S0169-7722(01)00174-7. Dagan, G., 1989. Flow and Transport in Porous Formations. Springer-Verlag, New York. Dagan, G., Cvetkovic, V., 1996. Reactive transport and immiscible flow in geological media. I. General theory. Proc. R. Soc. A 452, 1945. http:// dx.doi.org/10.1098/rspa.1996.0016. Davis, G.B., Barber, C., Power, T.R., Thierrin, J., Patterson, B.M., Rayner, J.L., Qinglong, W., 1999. The variability and intrinsic remediation of a BTEX plume in anaerobic sulphate-rich groundwater. Journal of Contaminant Hydrology 36, 3. http://dx.doi.org/10.1016/S0169-7722(98)00148-X. Dillard, L.A., Essaid, H.I., Herkelrath, W.N., 1997. Multiphase flow modeling of a crude-oil spill site with a bimodal permeability distribution. Water Resour. Res. 33, 7. http://dx.doi.org/10.1029/97WR00857. Domenico, P.A., 1987. An analytical model for multidimensional transport of a decaying contaminant species. J. Hydrol. 91, 1–2. http://dx.doi.org/ 10.1016/0022-1694(87)90127-2. Essaid, H.I., Herkelrath, W.N., Hess, K.M., 1993. Simulation of fluid distributions observed at a crude oil spill site incorporating hysteresis, oil entrapment, and spatial variability of hydraulic properties. Water Resour. Res. 29, 6. http://dx.doi.org/10.1029/93WR00370. Essaid, H.I., Bekins, B.A., Godsy, E.M., Warren, E., Baedecker, M.J., Cozzarelli, I.M., 1995. Simulation of aerobic and anaerobic biodegradation processes at a crude oil spill site. Water Resour. Res. 31, 12. http://dx.doi.org/ 10.1029/95WR02567. Essaid, H.I., Bekins, B.A., Herkelrath, W.N., Delin, G.N., 2011. Crude oil at the Bemidji site: 25 years of monitoring, modeling, and understanding. Ground Water 49, 5. http://dx.doi.org/10.1111/j.1745-6584.2009. 00654.x.

Farhadian, M., Vachelard, C., Duchez, D., Larroche, C., 2008. In situ bioremediation of monoaromatic pollutants in groundwater: a review. Bioresour. Technol. 99, 13. http://dx.doi.org/10.1016/j.biortech.2007.10.025. Fiori, A., 1996. Finite Peclet extensions of Dagan's solutions to transport in anisotropic heterogeneous formations. Water Resour. Res. 32, 1. http:// dx.doi.org/10.1029/95WR02768. Fiori, A., 2001. The Lagrangian concentration approach for determining dilution in aquifer transport: theoretical analisis and comparison with field experiment. Water Resour. Res. 37, 12. http://dx.doi.org/10.1029/ 2001WR000228. Fiori, A., Dagan, G., 1999. Concentration fluctuations in transport by groundwater: comparison between theory and field experiments. Water Resour. Res. 35, 1. http://dx.doi.org/10.1029/98WR01862. Fiori, A., Dagan, G., 2000. Concentration fluctuations in aquifer transport: a rigorous first-order solution and applications. J. Contam. Hydrol. 45, 1–2. http://dx.doi.org/10.1016/S0169-7722(00)00123-6. Gutierrez-Neri, M., Ham, P.A.S., Schotting, R.J., Lerner, D.N., 2009. Analytical modelling of fringe and core biodegradation in groundwater plumes. J. Contam. Hydrol. 107, 1–9. http://dx.doi.org/10.1016/j.jconhyd.2009.02.007. Ham, P.A.S., Schotting, R.J., Prommer, H., Davis, G.B., 2004. Effects of hydrodynamic dispersion on plume lengths for instantaneous bimolecular reactions. Adv. Water Resour. 27, 8. http://dx.doi.org/10.1016/ j.advwatres.2004.05.008. Kapoor, V., Gelhar, L.W., 1994. Transport in three-dimensionally heterogeneous aquifers: 1. Dynamics of concentration fluctuations. Water Resour. Res. 30, 6. http://dx.doi.org/10.1029/94WR00076. Liedl, R., Valocchi, A.J., Dietrich, P., Grathwohl, P., 2005. Finiteness of steady state plumes. Water Resour. Res. 41, 12. http://dx.doi.org/10.1029/ 2005WR004000. Lovley, D.R., Baedecker, M.J., Lonergan, D.J., Cozzarelli, I.M., Phillips, E.J.P., Siegel, D.I., 1989. Oxidation of aromatic contaminants coupled to microbial iron reduction. Nature 339 (6222), 297–299. http://dx.doi.org/10.1038/ 339297a0. Maier, U., Grathwohl, P., 2005. Numerical experiments and field results on the size of steady state plumes. J. Contam. Hydrol. 85, 1–2. http:// dx.doi.org/10.1016/j.jconhyd.2005.12.012. Pannone, M., Kitanidis, P.K., 1999. Large-time behavior of concentration variance and dilution in heterogeneous formations. Water Resour. Res. 35, 3. http://dx.doi.org/10.1029/1998WR900063. Rubin, J., 1983. Transport of reacting solutes in porous media: relation between mathematical nature of problem formulation and chemical nature of reactions. Water Resour. Res. 19, 5. http://dx.doi.org/10.1029/ WR019i005p01231. Rubin, Y., 2003. Applied Stochastic Hydrogeology. Oxford University Press, New York. Severino, G., Tartakovsky, D.M., Srinivasan, G., Viswanathan, H., 2012. Lagrangian models of reactive transport in heterogeneous porous media with uncertain properties. Proc. R. Soc. A 468, 2140. http:// dx.doi.org/10.1098/rspa.2011.0375. Suthersan, S.S., Payne, F.C., 2005. In Situ Remediation Engineering. CRC Press, Boca Raton, Florida. Tonina, D., Bellin, A., 2008. Effects of pore-scale dispersion, degree of heterogeneity, sampling size, and source volume on the concentration moments of conservative solutes in heterogeneous formations. Adv. Water Resour. 31, 2. http://dx.doi.org/10.1016/j.advwatres.2007.08.009. Van Breukelen, B.M., Griffioen, J., Roling, W.F.M., van Verseveld, H.W., 2004. Reactive transport modeling of biogeochemical processes and carbon isotope geochemistry inside a landfill leachate plume. J. Contam. Hydrol. 70, 3–4. http://dx.doi.org/10.1016/j.jconhyd.2003.09.003. Wexler, E.J., 1992. Analytical solutions for one-, two- and three-dimensional solute transport in ground-ware systems, Techniques of WaterResources Investigations of the 402 United States Geological Survey (USGS). US Government Print office. Zarlenga, A., Fiori, A., 2013a. Stochastic modelling of the length of steady plumes undergoing bioremediation. Procedia Environ. Sci. (ISSN: 18780296) 19, 633–642. http://dx.doi.org/10.1016/j.proenv.2013.06.072. Zarlenga, A., Fiori, A., 2013b. Steady plumes in heterogeneous porous formations: a stochastic Lagrangian approach. Water Resour. Res. 49, 864–873. http://dx.doi.org/10.1002/wrcr.20106.