## Abstract

We consider multi-component reactive transport in heterogeneous porous media with uncertain hydraulic and chemical properties. This parametric uncertainty is quantified by treating relevant flow and transport parameters as random fields, which renders the governing equations stochastic. We adopt a stochastic Lagrangian framework to replace a three-dimensional advection–reaction transport equation with a one-dimensional equation for solute travel times. We derive approximate expressions for breakthrough curves and their temporal moments. To illustrate our general theory, we consider advective transport of dissolved species undergoing an irreversible bimolecular reaction.

## 1. Introduction

Advection, molecular diffusion, hydrodynamic dispersion and (bio)chemical transformations are some of the key mechanisms that affect fate and transport of solutes in fluids flowing through porous media. Predictions of each of these phenomena, and of their combined effects are plagued with uncertainty. We focus on the impact of parametric uncertainty (uncertainty about parameter values characterizing macroscopic hydraulic and chemical properties of a porous medium); while equally important, quantification of uncertainty associated with the validity of macroscopic flow (Darcy's law) and transport (advection–dispersion models) equations lies outside the scope of the present analysis. To be specific, we ground the subsequent discussion in applications to transport in natural geological materials. Yet, the conceptual approaches and computational challenges addressed in the present study are equally applicable to manufactured and biological porous media.

Understanding the long-term migration of dissolved contaminants in the subsurface is critical to protecting the environment and the world's water supplies. Although subsurface hydrologists have been very successful at quantifying and predicting the bulk movement of subsurface water, predicting the movement of dissolved contaminants has proved to be more difficult. The transport of contaminants in natural aquifers is caused by advection due to groundwater flow driven by pressure gradients, molecular diffusion due to concentration gradients, hydrodynamic dispersion due to mixing of flow paths of varying velocities as a result of heterogeneities in the medium and chemical interactions between multiple reacting species with the soil matrix as well as themselves. Chemical and physical heterogeneities are inherently present in the porous medium on scales ranging from millimetres to kilometres. Transport modelling is further hampered by data sparsity, which leads to large epistemic uncertainties. The result is that subsurface properties cannot often be characterized at the scale of interest, and are considered to vary randomly.

Hydraulic conductivity and other properties of natural aquifers typically exhibit significant spatial variability over a wide range of scales. This natural heterogeneity gives rise to transport behaviours that are often non-Fickian (Neuman & Tartakovsky 2009). Transport of conservative (non-reacting) solutes in uncertain (random) velocity fields has been the subject of numerous studies, many of which are summarized in the research monographs by Dagan (1989), Cushman (1997) and Rubin (2003), to name just a few. Such stochastic analyses were carried out within either Eulerian (Dentz & Tartakovsky 2008) or Eulerian–Lagrangian (Neuman 1993) or purely Lagrangian (Cvetkovic & Dagan 1994) frameworks. The latter approach is employed in the present analysis.

Incorporation of chemical reactions into such stochastic models requires one to cope with uncertainty in conditions and parameters that control these reactions (e.g. reaction rate constants, distribution coefficient, the level of mixing and the surface areas of the reacting species). Srinivasan *et al.* (2007) discuss various types of uncertainty, both epistemic and aleatory, as pertains to geochemical reactions. Parametric uncertainty reflects our incomplete knowledge of the values of these parameters. This type of uncertainty refers to a discrepancy between parameter values obtained from laboratory and field experiments. A typical example is a difference of several orders of magnitude between field and laboratory estimates of reaction rate constants (Lichtner & Tartakovsky 2003 and references therein). Reactions may be between several dissolved species or between dissolved species and the solid matrix comprising a porous medium (e.g. reversible or irreversible sorption). Coupling the heterogeneity of the sorption parameters with that of the flow field has been the topic of several studies (Bellin *et al.* 1993; Berglund 1995; Cvetkovic *et al.* 1998). Recent assessments of parametric uncertainty in, and stochastic modelling of, reactive transport heterogeneous porous media include Sanchez-Vila & Rubin (2003), Dentz & Berkowitz (2005), Cirpka *et al.* (2008), Tartakovsky *et al.* (2009) and Broyda *et al.* (2010).

Our analysis is closely related to the Lagrangian models of reactive transport in randomly heterogeneous porous media proposed by Cvetkovic *et al.* (1998) and Cvetkovic & Dagan (1994, 1996). In accordance with these models, we disregard pore-scale dispersion and molecular diffusion because their effects on average large-scale characteristics of transport (e.g. travel times, and spatial and temporal concentration moments) are small when compared with advection. Our goal is to extend these analyses by accounting for the transport of multiple reacting species that undergo homogeneous and/or heterogeneous chemical reactions.

In §2, we provide a mathematical formulation of advective multi-component reactive transport in heterogeneous porous media, whose hydraulic and chemical properties are uncertain. Section 3 contains a derivation of our general approach, which enables one to express solute breakthrough curves (BTCs) at a control plane and their temporal moments in terms of the statistics of uncertain hydraulic conductivity and reaction rate constants. The approach results in closed-form analytical expressions, which makes it suitable for both parameter identification procedures (inverse modelling) and probabilistic risk assessment (Tartakovsky 2007). In §4, we demonstrate the salient features of the general approach by applying it to transport of solutes undergoing bimolecular chemical reactions. Section 5 consists of a summary of the key results.

## 2. Problem formulation

### (a) Flow model

Steady three-dimensional fluid flow in a porous medium with hydraulic conductivity *K* and porosity *ϑ* can be described by combining the Darcy law and mass conservation,
2.1
where **q**(**x**) is the Darcy velocity and *h*(**x**) is the hydraulic head. In heterogeneous porous media, both conductivity *K*(**x**) and porosity *ϑ*(**x**) vary in space (). Our analysis is motivated by transport in subsurface environments, wherein *K*(**x**) often varies by orders of magnitude and is heavily under-sampled. It is common (Dagan 1989; Rubin 2003 and references therein) to characterize this parametric uncertainty by treating *K*(**x**) as a random field. Despite some reservations (Gómez-Hernández & Wen 1998, Winter & Tartakovsky 2000, 2002), the standard practice is to assume that log-conductivity is a second-order stationary (statistically homogeneous) multi-variate Gaussian field with constant mean and variance , and a two-point correlation function *ρ*_{Y}(*r*) with *r*=|**r**|≡|**x**−**y**|. To be concrete, we choose
2.2
where *I*_{h} and *I*_{v} are horizontal and vertical integral scales, respectively. This statistical anisotropy is a salient feature of most subsurface environments, which reflects morphological processes leading to their formation (Dagan 1989; Tartakovsky & Neuman 1998; Indelman *et al.* 1999).

The spatial variability of porosity *ϑ*(**x**) is often much smaller than that of *K*(**x**) and, hence, is routinely ignored (Dagan 1989; Rubin 2003 and references therein). We follow this trend by treating *ϑ* as a deterministic constant. Finally, we restrict our analysis to flows that are driven by an externally imposed hydraulic head gradient **J**=(*J*,0,0)^{T} and are uniform in the mean, .

### (b) Transport model

#### (i) Eulerian description

A continuum-level Eulerian description of advective transport of a mixture of *N* dissolved species, with volumetric concentrations *C*_{i}(**x**,*t*), is provided by an advection–reaction equation (ARE),
2.3
Here, **C**≡(*C*_{1},…,*C*_{N}) is a set of *N* concentrations, and **V**(**x**) is the (random) average macroscopic velocity defined as **V**≡**q**/*ϑ*. Its ensemble mean is denoted by . The set of *N* functions ** Ψ**≡(

*Ψ*

_{1},…,

*Ψ*

_{N}) represents homogeneous and heterogeneous reactions between the species dissolved in the fluid and between these species and the solid matrix comprising the porous medium, respectively. Reaction rates in these reaction laws are uncertain and treated as random constants.

#### (ii) Lagrangian description

Suppose that at *t*≥0, a mixture of *N* dissolved species enters the system through an injection area *A*_{0} located in the plane. Both the initial concentrations of each species and (instantaneous or time-varying) injection rates are given. They can be either certain or uncertain.

Consider a point **a**∈*A*_{0} at *t*=0. Let **X**(*t*;**a**)=(*X*_{1},*X*_{2},*X*_{3})^{T} denote the position of the solute particle **a** at later times *t*>0, i.e. its trajectory. Because the flow velocity **V** is random, the particle trajectory **X**(*t*;**a**) satisfies a stochastic ordinary differential equation
2.4
Following Cvetkovic & Dagan (1994), we derive a Lagrangian alternative to (2.3) by replacing the Eulerian coordinates **x**=(*x*_{1},*x*_{2},*x*_{3})^{T} with their Lagrangian counterparts ** ξ**=(

*ξ*

_{1},

*ξ*

_{2},

*ξ*

_{3})

^{T}defined as 2.5 Here,

*τ*(

*x*

_{1};

**a**) denotes the time it takes the solute particle

**a**to travel from the injection plane (

*x*

_{1}=0) to the control plane (

*x*

_{1}>0). It is defined as a solution

*t*=

*τ*of

*x*

_{1}=

*X*

_{1}(

*t*;

**a**). The functions

*η*and

*ζ*are defined as

*η*(

*x*

_{1};

**a**)=

*X*

_{2}(

*τ*;

**a**) and

*ζ*(

*x*

_{1};

**a**)=

*X*

_{3}(

*τ*;

**a**), such that a streamline originating at the point

**a**is described by equations

*x*

_{2}=

*η*(

*x*

_{1};

**a**) and

*x*

_{3}=

*ζ*(

*x*

_{1};

**a**).

Molecular diffusion and hydrodynamic dispersion act to displace solute particles from one streamline to another. Because these transport mechanisms are neglected in the present analysis, a fluid particle (e.g. **a**) stays on the same streamline, so that *ξ*_{2}=*ξ*_{3}=0. Therefore, rewriting (2.3) in the ** ξ** coordinate system (2.5) yields a transport equation
2.6
which applies to any stream-tube originating from

**a**∈

*A*

_{0}. The concentrations

**C**=(

*C*

_{1},…,

*C*

_{N}) are now functions of the Lagrangian variables

*τ*and

*t*(although for simplicity we keep the same notation). One implication of the Lagrangian formulation (2.6) is that the three-dimensional Eulerian equations (2.3) have been recast in a one-dimensional form because the original three-dimensional nature of the flow and transport is entirely encapsulated in the travel time

*τ*.

A main goal of the subsequent analysis is to express the ensemble statistics of travel time *τ* and corresponding BTCs in terms of the statistics of the input parameters, i.e. log-conductivity *Y* (**x**) and reaction rates. This goal amounts to propagating parametric uncertainty through the modelling process, thus quantifying predictive uncertainty.

### (c) Modelling assumptions

Our analysis is based on the following assumptions.

—

*Porous media are assumed to be mildly heterogeneous.*The assumption of mild heterogeneity refers to the requirement that the variance of log-conductivity be small, . It has implications for both the mathematical developments in §3*a*and the validity of Fickian models of transport on a large (e.g. field) scale. The former employ as a perturbation parameter. The latter has to do with the impact of heterogeneity on a complex interplay of mixing, spreading and chemical reactions (Bellin*et al.*2011; Le Borgne*et al.*2011). The range of applicability of the perturbation solution presented in §3*a*can be extended to highly heterogeneous environments (large ) by employing the random domain decomposition of Winter & Tartakovsky (2000, 2002). Analysis of non-Fickian transport models lies outside the scope of the present study, which deals with the effects of parametric (rather than model or structural) uncertainty.—

*Local-scale molecular diffusion and hydrodynamic dispersion are neglected.*We assume that the flow velocity**V**is sufficiently large for molecular diffusion to be negligible, and define the ARE (2.3) on a spatial scale that is sufficiently small to neglect local-scale hydrodynamic dispersion. These transport mechanisms give rise to (transverse) mixing that can significantly affect reactive transport when the characteristic initial plume size is smaller than the horizontal integral scale of the hydraulic conductivity, and a subsurface environment is highly heterogeneous (Bellin*et al.*2011). Quantities most affected by transverse spreading are higher moments of the probability density function (PDF) of solute concentration (Bellin*et al.*2011). We neglect local-scale dispersion because our analysis both assumes that porous media are mildly heterogeneous and deals with the first moment (BTCs). This issue is discussed further in Tartakovsky*et al.*(2009) and references therein.—

*Uncertain reaction rate constants are treated as random variables.*Reaction rate constants of homogeneous reactions are affected by thermodynamics and kinetics. Because uncertainty in such rate constants=(*κ**κ*_{1},…,*κ*_{M}) largely stems from measurement and/or interpretive errors, which are roughly the same in all samples collected throughout a porous medium, the treatment ofas random variables is warranted. Reaction rate constants of heterogeneous reactions also depend on the (highly variable) surface area of a porous matrix. Uncertainty in spatial variability of these constants is typically quantified by treating them as a spatially correlated random field*κ*(*κ***x**) (Tartakovsky*et al.*2009; Tartakovsky & Broyda 2011 and references therein). In this study, we are concerned with a typical situation in which only a very few sparse measurements ofare available so that any statistics beyond their mean and variance become meaningless. This practical situation calls for the treatment of*κ*(*κ***x**) as a random variable whose mean and variance are determined from the measurements.—

*Fluid flow is unaffected by chemical reactions*. In general, chemical processes (e.g. mineral dissolution and precipitation) can change hydraulic properties of a porous medium (e.g. its porosity and permeability). However, these changes typically occur on time scales that are orders of magnitude larger than the time scale associated with the chemical reactions under consideration. Consequently, we assume that chemical reactionsdo not affect the flow velocity*Ψ***V**.—

*Flow is uniform in the mean*. The assumption of constant mean flow velocity is made to simplify the presentation, and is not intrinsic to the proposed approach. Accounting for spatial variability of the mean velocity in Lagrangian models of reactive transport is relatively straightforward. Consider, for example, the ease with which the Tartakovsky*et al.*(2009) solution derived for mean uniform flow was generalized by Broyda*et al.*(2010) to accommodate radial flow conditions.

## 3. General theory

In the absence of molecular diffusion and local hydrodynamic dispersion, the *x*_{1}-component of the mass flux of the *i*th species is given by *C*_{i}(**x**,*t*)*V* _{1}(**x**). Because the flow is uniform in the mean, the ensemble mean of the flux-averaged concentration can be defined as
3.1
Here and below, we use the notation and interchangeably to denote the conditional ensemble mean of a random quantity over all possible realizations of random log-conductivity *Y* (**x**), for a given realization of the reaction rates. The ergodicity hypothesis postulates the equivalency between the ensemble mean and the spatial average (in our case, over the injection domain *A*_{0}),
3.2
When plotted as a function of time *t*, the mean flux-average concentration in (3.2) is called a BTC at a control plane *x*_{1}. It represents the total concentration of the *i*th species in the fluid that is extracted (after complete mixing) over a large area of the control plane. Note that the BTCs defined by (3.2) are random because they depend on the uncertain (random) reaction rates.

Suppose that the set of chemical reactions ** Ψ** is controlled by a set of

*M*uncertain (random) reaction constants

**=(**

*κ**κ*

_{1},…,

*κ*

_{M}) with a joint PDF . Let denote a PDF of the travel time

*τ*that is related to the statistics of

*g*

_{Y}, the PDF of log-conductivity

*Y*, below. Then, , the ensemble mean of the concentration averaged over all possible realizations of

**, is given by 3.3**

*κ*### (a) Travel-time probability density function

Up to first order in , the PDF can be expressed in terms of the variance of the horizontal component of the trajectory of non-reactive particles, , and the corresponding longitudinal macro-dispersion coefficient, 3.4 as (Cvetkovic & Dagan 1996) 3.5a where 3.5b and 3.5c Far away from the injection plane, the variance grows linearly with time , so that the term reduces to the lognormal PDF predicted by Dagan (1989, eqn (5.8.9)). Thus, it follows from (3.5) that travel times do not obey lognormal PDFs, except as an approximation far away from the injection plane.

The probability of a solute crossing the control plane *x*_{1} no later than time *t*, i.e. the probability of the arrival time *τ* not exceeding *t*, is obtained from (3.5) as
3.6

The variance of the horizontal component of the trajectory of non-reactive particles can be obtained via a perturbation expansion in (Dagan 1989). Instead, we follow Dagan & Cvetkovic (1993) and Severino *et al.* (2005) by approximating the variance with
3.7a
where
3.7b
and
3.7c
In geological formations, the anisotropy ratio 0≤*λ*≤1. It is easy to verify that *χ*(1)≤*χ*(*λ*)≤*χ*(0) with *χ*(0)=1 and *χ*(1)=8/15.

Severino *et al.* (2005) demonstrated that this approximation of is in close agreement with its perturbation-based counterparts (Dagan 1989). They found that the difference between the two expressions is less than 0.5 per cent for *λ*<1 and less than 3 per cent for *λ*=1.

The travel time PDF is now given by (3.5) with (3.7) and (3.4). The variance of this distribution has a simple form (appendix A)
3.8
It sheds light onto the uncertainty mechanisms affecting travel-time predictions. Specifically, the predictive uncertainty vanishes, i.e. , when *β*≪1. This occurs if either the travel distance between the injection and control plane is much smaller than the horizontal correlation length of hydraulic conductivity (*x*_{1}/*I*_{h}≪1) or a porous medium is perfectly stratified (). On the other hand, the predictive uncertainty grows linearly with the travel distance, , when *β*≫1. This scenario corresponds to either travel distances that are much larger than the horizontal correlation length of hydraulic conductivity (*x*_{1}/*I*_{h}≫1) or poorly correlated porous media (). As discussed earlier, the latter regime (*β*≫1) allows for the approximation , i.e. the travel-time distribution is approximately lognormal.

### (b) Temporal moments

In addition to the mean BTCs of the *i*th dissolved species, , their temporal moments
3.9
are of practical interest (Cirpka & Kitanidis 2000). We will analyse the behaviour of the normalized moments
3.10a
and
3.10b

## 4. Example: irreversible bimolecular reactions

To demonstrate the salient features of the general theory described above, we consider advective transport of dissolved species *A* and *B* that undergo an irreversible elementary chemical reaction to form a product *C*,
4.1
Despite its relative simplicity, the reaction (4.1) describes a number of geochemical processes (Srinivasan *et al.* 2007). A rate equation for the reaction (4.1) is
4.2
where *C*_{γ} denotes the concentrations of the species *γ* (*γ*=*A*,*B*,*C*), and *κ* is the rate coefficient. Because the latter depends on the surface area of the adsorbent, in porous media applications, it is highly uncertain and is treated as a random variable. Substituting (4.2) into (2.6) yields
4.3
The Lagrangian equations (4.3) are defined on the semi-infinite domain *x*_{1}>0, and are subject to the initial and boundary conditions
4.4
where and are, respectively, the initial and boundary concentration of the species *γ*.

### (a) Solute concentrations

Let us introduce a conservative concentration component *u*(*τ*,*t*)≡*C*_{A}(*τ*,*t*)−*C*_{B}(*τ*,*t*). It follows from (4.3) and (4.4) that *u*(*τ*,*t*) satisfies subject to the initial and boundary conditions
4.5
Solving this boundary-value problem yields
4.6
where *H*(*x*) is the Heaviside step function. Substituting *C*_{A}=*u*+*C*_{B}, with *u* given by (4.6), into (4.3) and using the method of characteristics to solve the resultant ordinary differential equation, we obtain
4.7
4.8
and
4.9
where *s*≡*τ*−*t*, and .

### (b) Breakthrough curves

Substituting (4.7)–(4.9) into (3.3), we obtain the BTCs of the species *γ* (*γ*=*A*,*B*,*C*),
4.10a
where
4.10b
4.10c
and
4.10d
and the functions (*γ*=*A*,*B*,*C*) are obtained from their respective counterparts *Γ*^{in}_{γ}(*τ*,*t*) by replacing *s*≡*τ*−*t* with −*s* in (4.10b)–(4.10d).

The BTCs in (4.10) can be thought of as a superposition of two transport features. The first term in (4.10a) represents the contribution of the streamlines along which the solute particles at time *t* have already reached the control plane (i.e. *τ*>*t*). The second term in (4.10a) represents the contribution of the remaining streamlines along which the solute particles have not reached the control plane (i.e. *τ*<*t*).

Evaluation of the BTCs in (4.10) requires the knowledge of the PDFs of travel times *τ* and reaction rates *κ*. While the former, , is given by (3.5), (3.7) and (3.4), the latter, , is to be specified as an input parameter reflecting the uncertainty about *κ*. Below, we analyse reactive transport with both certain (deterministic) and uncertain (random) reaction rates.

In both cases, the initial concentrations of the two reactants are set to
4.11
where . Here, *c*^{in}_{A} and *α*_{in} are known (deterministic) constants, and the constants *τ*_{a} and *τ*_{b} are defined by the interval *x*∈[*x*_{a}≡*Uτ*_{a},*x*_{b}≡*Uτ*_{b}], in which the initial concentrations of the reactants are non-zero. The initial concentration of the reactant product is set to *C*^{in}_{C}=0.

#### (i) Deterministic reaction rate

We start by assuming that *κ* is known with certainty (deterministic), i.e. where *δ*(⋅) is the Dirac delta function. For the sake of concreteness, we consider (4.4) with homogeneous boundary conditions .

Under these conditions, (4.10) admits a closed-form analytical solution,
4.12a
where
4.12b
4.12c
and
4.12d
with *G* given by (3.6). The factor `P`(*t*;*x*_{1}) represents the impact of uncertainty in hydraulic conductivity on estimation of the BTCs, whereas the factor *χ*_{γ}(*t*) quantifies the impact of the geochemical reaction. It follows from (4.12) that the time *t*^{☆} taken by the product concentration to reach the concentration of a limiting reactant (say, *B*) is given by . This result is useful for validating more evolved numerical codes.

To facilitate a physical interpretation of the BTCs in (4.12), we start by considering transport of conservative solutes. This case corresponds to *κ*=0, which yields *b*(*t*)=*α*_{in} and for the non-reacting species *γ*=*A*,*B*. For uniform initial concentrations throughout the flow domain, i.e. for *τ*_{a}=0 and , (4.12d) yields `P`(*t*;*x*_{1})=1−*G*(*t*;*x*_{1}). Thus, the normalized BTC, , represents the probability of a solute particle that originated at the injection plane *x*_{1}=0 at time *t*=0 not crossing the control plane *x*_{1} by time *t*, i.e. the probability that *τ*≥*t*. Identical interpretations of BTCs were reported by Dagan & Nguyen (1989) for these transport conditions, by Rodriguez-Iturbe & Rinaldo (1997) for transport in rivers and by others. Under the more general conditions represented by (4.12), chemical reactions act to reduce the ‘probability’ of *τ*≤*t* (as quantified by the normalized BTCs ) by the factor *χ*_{γ}.

Figure 1 exhibits the temporal behaviour of `P`(*t*;*x*_{1}) at control planes *x*_{1}=5*I*_{h}, 15*I*_{h} and 20*I*_{h}, for several values of the log-conductivity variance . When values of the hydraulic conductivity are certain (), `P`(*t*;*x*_{1}) has two jump discontinuities at *t*=*x*_{1}/*U*−*τ*_{b} and *t*=*x*_{1}/*U*−*τ*_{a}. Uncertainty in the hydraulic conductivity () smoothes `P`(*t*;*x*_{1}), leading to asymmetric curves. The asymmetry increases with and decreases with the normalized distance between the injection and control planes, *x*_{1}/*I*_{h}. According to (4.12), this translates into a large proportion of early arrival times for *x*_{1}/*I*_{h}=5, which grows significantly as becomes larger (figure 1*a*). As the travel distance increases (*x*_{1}/*I*_{h}=15 and 20), the likelihood of encountering early arrival times decreases (figure 1*b*,*c*) owing to the plume's self-averaging as it samples more and more (statistically homogeneous or stationary) heterogeneities. Although not shown here, the impact of the anisotropy ratio *λ* on the temporal behaviour of `P`(*t*;*x*_{1}) is discernible only for *x*_{1}/*I*_{h}≪1. As, in most practical applications, one is interested in travel distances that exceed the correlation length of hydraulic conductivity (*x*_{1}/*I*_{h}≥1), the effect of *λ* is of minor importance.

To quantify the impact of chemical reactions on the BTCs in (4.12), i.e. to explore the parameter space controlling the temporal variability of *χ*_{γ}(*t*), we introduce the Damköhler number *Da*=*κc*^{in}_{A}*I*_{h}/*U* as the ratio of the reaction (1/*t*_{r}=*κc*^{in}_{A}) and advection (*t*_{h}=*U*/*I*_{h}) time scales. Then, time dependence of *χ*_{γ}(*t*) in (4.12) arises from the function where *t*_{d}=*t*/*t*_{h} is the dimensionless time. Note that mass conservation introduces the following constraints on *χ*_{γ}(*t*_{d}) (*γ*=*A*,*B*,*C*): at all times *t*_{d}, *χ*_{A}(*t*_{d})−*χ*_{B}(*t*_{d})=1−*α*_{in} and *χ*_{A}(*t*_{d})+*χ*_{C}(*t*_{d})=1.

Figure 2 demonstrates the temporal variability of *χ*_{γ}(*t*_{d}) (*γ*=*A*,*B*,*C*) in the advection-dominated (*Da*=10^{−2}), reaction-dominated (*Da*=10^{2}) and intermediate (*Da*=1) regimes. In all cases, we set *α*_{in}=1/2, which renders *B* a limiting reactant. In the advection-dominated regime (*Da*=10^{−2}), the reaction is negligible (*χ*_{γ}'s remain approximately constant) until *t*_{d}≥10. By this time, the plume migrates far enough (distances that are two orders of magnitude larger than *I*_{h}) to sample the medium's heterogeneity, i.e. solute transport becomes ergodic. In the reaction-dominated regime (*Da*=10^{2}), the reaction occurs before the plume samples the medium's heterogeneity. Thus, the value of the Damköhler number *Da* determines whether the reaction took place and/or solute transport became ergodic before the plume crosses the control plane.

The BTCs at *x*_{1}=5*I*_{h} computed with the first equation in (4.12) are shown in figure 3. The impact of the chemical reaction is seen within the time interval 10^{−3}≤*t*_{d}≤10. Thus, for *t*<*t*_{d}, the species *A* and *B* behave as conservative tracers. For *t*>*t*_{d}, the species *B* does not exist anymore, and the remaining species *A* and *C* migrate without undergoing chemical tranformations.

#### (ii) Random reaction rates

Let the reaction rate *κ* be uncertain, distributed uniformly on the interval with mean and half-band *Δ*, i.e. for and =0 for or . We consider boundary conditions representing the injection of the reactant *A* and *B* during a finite time interval 0<*t*≤*T*_{0},
4.13
where and are constants. In the simulations reported below, we set , *λ*=0.1 and .

Substituting the uniform into (4.10b)–(4.10d) leads to
4.14a
and
4.14b
where *s*=*τ*−*t* and
4.14c
The functions (*γ*=*A*,*B*,*C*) are obtained from their respective counterparts *Γ*^{in}_{γ}(*τ*,*t*) by replacing *s* with −*s* in (4.14). Note that taking the limit , one recovers the results obtained in §3 for deterministic *κ*.

Combining (4.14) and (4.10) yields the following expressions for the BTCs:
4.15a
where
4.15b
and
4.15c
Here, , ,
4.15d
and is the coefficient of variation of *κ*. The function is obtained by replacing *α*_{in} with *α*_{0} and *u*_{in} with *u*_{0} in the expression (4.15d).

Figure 4 depicts the temporal variability of the BTCs at the control plane *x*_{1}=5*I*_{h}, in the reaction-dominated regime (*Da*=10^{2}, where the Damköhler number is now defined as ). Comparison of the curves corresponding to the coefficient of variation *ξ*_{κ}=0.01 and 1.0 reveals that uncertainty in the reaction rate constant affects the BTC predictions only at times . In other words, the chemical reaction affects the concentrations of solutes crossing the control plane *x*_{1}/*I*_{h}=5 only during this time interval, regardless of whether the reaction rate is certain (figure 3) or uncertain (figure 4). This is because for *t*_{d}<10^{−2}, the reaction has not yet started, whereas for *t*_{d}>1, the reaction is already exhausted.

The impact of uncertainty in *κ* can be quantified in terms of the global relative difference, over a time interval 0≤*t*≤*T*, between evaluated for *ξ*_{κ}=1.0 () and *ξ*_{κ}=0.01 (),
4.16
The BTCs of species *B*, which is the limiting reactant, are most sensitive (*Δ*_{B}>100%) to the parametric uncertainty. The corresponding global differences for the BTCs of species *A* and *C* are *Δ*_{A}≈20 per cent and *Δ*_{C}≈55 per cent, respectively.

The results presented above can be augmented to answer the following question. For a given value of the Damköhler number *Da*, what is the maximum size of a flow domain (i.e. the maximum position of the control plane *x*_{1}=*ηI*_{h}, where *η* is a constant) within which transport remains reaction dominated? Let denote a characteristic reaction length. Then, the control plane coordinate can be written as *x*_{1}=*ηDa*ℓ_{r}. According to the very definition of ℓ_{r}, the transport is reaction-dominated inside flow domains whose horizontal coordinates are . This yields a condition . Hence, for *Da*≫1, transport is dominated by reactions within a region whose horizontal extent is much smaller than the horizontal correlation length *I*_{h} (*η*≪1). For *Da*≪1, the control plane must be located at a distance that is much larger than the horizontal integral scale *I*_{h} (*η*≫1) for the reaction's impact to be noticeable.

Figure 5 illustrates this analysis by exhibiting the BTCs corresponding to *Da*=10^{−2}, at two control planes: *x*_{1}=100*I*_{h} (*ηDa*=1) and *x*_{1}=5*I*_{h} (*ηDa*=0.05). The BTCs at *x*_{1}=100*I*_{h} are affected by uncertainty in the reaction rate constant, with the global differences of *Δ*_{A}≈14 per cent, *Δ*_{B}≈36 per cent and *Δ*_{C}≈33 per cent. At *x*_{1}=5*I*_{h}, the profiles of and are identical—up to the factors of 1−*α*_{in}, with *α*_{in} representing the mass partition coefficients for *A* and *B*, respectively—to their respective counterparts in figure 1*a* (the curves corresponding to ). The reaction product *C* is not yet produced by this time.

### (c) Temporal moments of breakthrough curves

Temporal moments of BTCs are useful in practical applications, in which the complete evaluation of BTCs is not feasible owing to an insufficient number of samples. In this section, we analyse the behaviour of the first three moments (3.10) for *Da*=0.1. According to the discussion above, in this transport regime, uncertainty about the reaction rate constant should affect only the BTCs of the limiting reactant *B* at control planes *x*_{1}≥10*I*_{h}.

Figure 6 shows the normalized first moments (temporal means) *T*_{γ}/*t*_{h}, the normalized second moments (temporal variances) and the third moments (skewness) *S*_{γ} as functions of the dimensionless position of the control plane, *x*_{1}/*I*_{h}. The BTC moments of the reacting species *A*, computed for the coefficients of variation *ξ*_{κ}=0.01 and 1.0, visually coincide, as do their counterparts for the reaction product *C*. At the same time, the BTC moments of the limiting reactant *B* are visibly affected by the degree of uncertainty in the reaction rate constant (values of *ξ*_{κ}). This behaviour reflects the fact that while the uncertainty influences at all times until the reaction is completed, its effect on and is limited to the time interval *t*_{d}≈10÷100. Outside this interval, they behave like conservative tracers. The latter observation is made explicit in figure 6*b* via comparison with the normalized variance of travel times of a conservative tracer (3.8). For distances *x*_{1}/*I*_{h}≥10, the slope of the latter curve is nearly identical to those for *A* and *C*.

Reactions manifest themselves in higher order moments, as one can clearly see from the non-zero skewness (figure 6*c*). Although for the illustrated curves, skewness is positive, it may attain negative values as well, when non-zero feeding (i.e. boundary) conditions are accounted for.

## 5. Summary and concluding remarks

We generalized the stochastic Lagrangian framework of Cvetkovic & Dagan (1994) to model multi-component reactive transport in heterogeneous porous media with uncertain hydraulic and chemical properties. We derived general expressions for the solute BTCs and their temporal moments. The proposed framework differs from that proposed by Cvetkovic *et al.* (1998) who dealt with the spatial variability of sorption parameters. In the present paper, we quantified parametric uncertainty in predictions of transport of two aqueous solutes undergoing irreversible chemical reactions. This methodology can be applied to multiple aqueous species participating in multiple reversible or irreversible reactions. We derived analytical solutions for a dissolution model that is frequently encountered in the field. The solutions were obtained by transforming a three-dimensional Eulerian domain (**x**,*t*) into a one-dimensional (*τ*,*t*) Lagrangian-domain (*τ*,*t*), where *τ* is the stochastic travel time of a particle released from an injection plane and recovered at a given control plane at *t*=*τ*.

Our analysis leads to the following major conclusions.

— The concentrations of reactants and their product are simultaneously influenced by the distribution of the travel time

*τ*and a geochemical irreversible reaction.— The variance of the time it takes particles to travel from the injection plane to the control plane quantifies uncertainty in the advective component of the solute transport.

— The Damköhler number

*Da*, a dimensionless parameter defined as the ratio of the advection and reaction time scales, determines the relative importance of various sources of parametric uncertainty. For large*Da*, i.e. when reactions are much faster than advection, uncertainty in the reaction rate constants dominates predictive uncertainty. For small*Da*, this source of parametric uncertainty is negligible.— Our analysis of the BTCs of the reactants and their product concentrations reveals that for given

*Da*the control plane should be located at*x*_{1}=*I*_{h}/*Da*(where*I*_{h}is the integral scale of log-conductivity) in order to detect the reaction regime.— In this regime, the moments of the limiting reactant are more sensitive than those of either the reactant present in excess or the product. This is due to the fact that heterogeneity plays a bigger role in their transport.

## Acknowledgements

This study was supported in part by subcontract no. 80227-001-10 of the LANL IGPP Programme and by the DOE Office of Science, Advanced Computing Research (ASCR) programme in Applied Mathematical Sciences. The first author acknowledges support from ‘Programma di scambi internazionali per mobilità di breve durata’, Naples University, Italy.

## Appendix A. Derivation of travel-time variance

We derive an approximate (closed-form) expression for the variance of travel times. Consistent with the approximation that has led to the longitudinal variance of the particle trajectory for non-reactive solutes (Dagan & Cvetkovic 1993), we approximate the longitudinal velocity covariance by
A1
Severino *et al.* (2005) demonstrated that the variance coincides with the first-order (in ) approximation of . The approximation (A1) also preserves the horizontal integral scale (Dagan & Cvetkovic 1993).

Up to the first order in , the travel-time variance can now be computed as (Cvetkovic & Dagan 1996) A2 Substituting (A1) into (A2) leads to (3.8).

## Footnotes

↵† On leave from University of Naples.

- Received June 12, 2011.
- Accepted October 25, 2011.

- This journal is © 2011 The Royal Society