## Abstract

Most stochastic analyses of reactive transport in physically and geochemically heterogeneous aquifers have focused on the analysis of a single reactive species. Here we conduct the stochastic analysis of multicomponent competitive monovalent cation exchange. Transport equations for dissolved cations are coupled with nonlinear cation exchange terms, which, for chemical equilibrium, are described by mass-action law expressions. These equations can be effectively decoupled by assuming that the weighted sum of cation concentrations is constant. The weight of each cation is equal to the reciprocal of its selectivity. Randomness of cation exchange capacity (CEC) leads to random retardation factors. Analytical expressions for effective retardation factors, longitudinal macrodispersivities, and concentration spatial moments are derived for a chemical system made of three monovalent cations (Na^{+}, K^{+}, and Cs^{+}) using the stochastic analytical solution of Miralles-Wilhelm and Gelhar (1996). Our results indicate that effective retardation factor, *R _{C,i}*, spatial moments, and macrodispersivities of K

^{+}are significantly different from those of Na

^{+}. Effective retardation factors asymptotically attain their mean values after a transient phase of cation-dependent duration. They strongly depend on the correlation between log-permeability (log

*K*) and CEC. Pre-asymptotic effective retardation factor values for a negative correlation are smaller than the mean value, regardless of the value of the coefficient of variation of CEC (

*CV*

_{CEC}). The smaller (larger) the variance of log

*K*,

^{+}is smaller (larger) than that of K

^{+}for a negative (positive) correlation structure. The macrodispersivity of K

^{+}increases with

*k*

_{1}and decreases with the Na

^{+}/Cs

^{+}selectivity,

*k*

_{2}, for a positive correlation. The first-order spatial moment of Na

^{+}is greater than that of K

^{+}, but is smaller than that of a nonreactive species at the asymptotic phase. Second-order spatial moments of cations are smaller than that of a nonreactive species for a positive correlation structure, but are larger for a negative correlation.

## INTRODUCTION

A vast amount of research has been conducted in the last two decades to analyze the impact of geological heterogeneity on contaminant transport (Bosma et al., 1993; Brusseau, 1998; Dentz et al., 2000a, 2000b; Espinoza and Valocchi, 1997; Gelhar, 1993; Ginn et al., 1995; Kapoor et al., 1997; Miralles-Wilhelm and Gelhar, 1996; Mojid and Vereecken, 2005; Pang et al., 2003; Simmons et al., 1995; Srivastava et al., 2002, 2004). Experimental evidence and theoretical results indicate that the transport of reactive solutes in porous media can be significantly influenced by the spatial variability in physical and geochemical parameters, such as hydraulic conductivity, porosity, and sorption coefficient. Ptak et al. (2004) summarized recent developments on field tracer tests in heterogeneous porous media and stochastic modeling of flow and transport.

Contaminant mixtures are usually considered to contain two or more distinct compounds that chemically interact with one another as well as with solids in a highly nonlinear manner. One of the most challenging aspects is the derivation of effective parameters at field scale for reactive species that migrate through heterogeneous porous formations. Miralles-Wilhelm et al. (1997) argued that the assumption of equal dispersivities for substrate and oxygen could be inadequate because they found longitudinal macrodispersivities for the contaminant and dissolved oxygen to be significantly affected by biodegradation. In comparison with the large amount of stochastic studies for a single reactive solute transport with adsorption, much less attention has been paid to multicomponent competitive cation exchange reactive transport in physically and geochemically heterogeneous aquifers. Yang and Samper (2004) and Samper and Yang (2004) analyzed the behavior of multicomponent cation exchange in a heterogeneous column by Monte-Carlo simulations. They pointed out that apparent dispersivities are different for each reactive cation.

In this paper, we evaluate the effects of hydrodynamic and geochemical heterogeneities on the transport of a multicomponent set of chemical species undergoing competitive monovalent cation exchange. A synthetic case with three monovalent cation exchange reactions is considered. K^{+} and Na^{+} are injected into a physically (hydraulic conductivity) and geochemically (cation exchange capacity [CEC]) heterogeneous aquifer in which Cs^{+} is the initially dominant cation occupying the exchange positions of solid phase. Na^{+} and K^{+} competitively exchange with Cs^{+} and enter into the solid phase. Cation exchange in groundwater is a surface reaction in which mass transfer between solution and solid commonly proceeds according to nonlinear mass-action equations. Cation transport equations, which are generally coupled due to nonlinear cation exchange sink/source terms, can be effectively decoupled when the total weighted dissolved cation concentration (*Q* = *C*_{1} + *C*_{2}/*k*_{1} + *C*_{3}/*k*_{2}) is constant, where *C _{i}* is dissolved cation concentration, and

*k*is selectivity coefficient. Retardation factors for each cation are derived in terms of CEC and selectivity coefficients. Randomness of retardation factors stems from the randomness of CEC. Mean and variance of retardation factors are calculated from those of CEC. The original problem is decoupled into three transport equations with random groundwater velocity and different random retardation factors. Analytical expressions of effective retardation factor, longitudinal macrodispersivity, and concentration spatial moments derived by Miralles-Wilhelm and Gelhar (1996) are employed to analyze the stochastic behavior of Na

_{i}^{+}and K

^{+}.

We describe here the mathematical formulation of the problem, and then we present the analytical expressions for effective retardation factors, macrodispersivities, and concentration spatial moments for Na^{+} and K^{+}. Effective retardation factors, macrodispersivities, and concentration spatial moments of Na^{+} and K^{+} are evaluated in terms of correlation structure, variance, and correlation length of log *K* and selectivity coefficients.

## MATHEMATICAL STATEMENT OF STOCHASTIC CATION TRANSPORT AND EXCHANGE

### Cation Transport Equations

*N*dissolved cations undergoing cation exchange in porous media can be expressed as: where

*C*is concentration of the

_{i}*i*th dissolved cation (mol/L);

*W*is concentration of the

_{i}*i*th exchanged cation in the solid phase (mol/L);

**D**is the local dispersion tensor (L

^{2}/T); and

**v**is the water velocity vector (L/T).

### Cation Exchange Reactions

^{+}, K

^{+}, and Cs

^{+}can be expressed as where (X-Na), (X-K), and (X-Cs) denote exchanged cations. It should be noted that these three reactions may take place simultaneously; however, they are linearly dependent. In fact, the third one can be obtained by adding the first two equations. Therefore, only two out of three possible exchange reactions are sufficient to describe all exchange reactions. Application of mass-action laws to reactions 2 and 3 leads to: where

*k*

_{1}and

*k*

_{2}are selectivity coefficients of Na

^{+}to K

^{+}and Na

^{+}to Cs

^{+}, respectively, and β

_{Na}, β

_{K}, and β

_{Cs}are equivalent fractions of exchanged Na

^{+}, K

^{+}, and Cs

^{+}in the solid phase, respectively, which add to one:

*W*, through where CEC is cation exchange capacity (meq/100 g of solids), 𝛉 is porosity (−),

_{i}*z*is charge (in this case,

_{i}*z*= 1 ∀

_{i}*i*), ρ

_{d}is dry density of medium (kg solid/dm

^{3}), and

*i*= Na

^{+}, K

^{+}, and Cs

^{+}.

### Random Fields

*C*, and exchanged cations,

_{i}*W*, as well as groundwater velocity,

_{i}**v**, in equation 1 are modeled as random fields. These random fields are driven by the inherent physical and geochemical heterogeneities of aquifers. Hydraulic conductivity,

*K*, is typically found to be log-normally distributed, so that in order to describe its spatial variability, it is convenient to work with its logarithm:

*f*(

*x*) can be expressed as the sum of its ensemble mean

*f̄*and a random perturbation

*f*′: Here

*f*′ is assumed to be a statistically homogeneous random field, described by a three-dimensional anisotropic exponential autocovariance function with correlation lengths λ

_{i}(

*i*= 1, 2, 3) and a variance

*K*according to where

*a*and

*b*are the intersect and the slope of the linear correlation of CEC and log-conductivity,

*f*, and η is a zero-mean random residual accounting for the fact that some of the variance of CEC is not related to the variability of

*f*. It should be noticed that

*b*> 0 indicates a positive correlation,

*b*< 0 indicates a negative correlation, and

*b*= 0 means no correlation between CEC and

*f*.

In our analysis, CEC and *f* are assumed to be perfectly correlated, and therefore η is equal to 0.

## SOLUTION METHODOLOGY

### Decoupled Cation Transport Equations

Contrary to the linear sorption problems usually discussed in the literature, the problem described in the previous section has nonlinear chemical source/sink terms involving concentrations of several components. Therefore, the first step is to linearize source/sink terms so that transport equations can be decoupled.

*W*) in equation 1 can be expressed in terms of solute concentrations by solving equations 5, 6, 7, and taking into account equation 8. According to the derivation given in Appendix A, exchanged concentrations,

_{i}*W*

_{Na},

*W*

_{K}, and

*W*

_{Cs}can be written as where

In the case study analyzed here, it is assumed that Cs^{+} occupies initially all exchange sites. After the injection of Na^{+} and K^{+} into the system, these cations are sorbed to the solid phase, while Cs^{+} is released to the aqueous solution. Dissolved concentrations of Na^{+} and K^{+} are much smaller than that of Cs^{+}, so Na^{+} and K^{+} are trace cations compared to Cs^{+}. The focus of this paper is the stochastic analysis of Na^{+} and K^{+}.

*Q*is constant, one finds where

*L*(

*C*) = ∇ · (

_{i}**D**∇

*C*) –

_{i}**v**∇

*C*, and

_{i}It should be noted that the two equations in equation 17 are uncoupled stochastic partial differential equations similar to those usually found in the literature. *R*_{Na} and *R*_{K} are retardation factors for Na^{+} and K^{+}, respectively, while CEC/100*Q* and CEC/100*Qk*_{1} are the distribution coefficients of Na^{+} and K^{+}, respectively.

### Validity of Constant *Q*

Our solution assumes that the weighted sum of cation concentrations, *Q*, is constant. The weight of each cation is equal to the reciprocal of its selectivity. Jin and Ye (1999) used a similar assumption when deriving an analytical solution for binary cation exchange.

*Q*is discussed by Appelo and Postma (1993). They present calculations corresponding to a set of divalent cations including Ca

^{+2}, Mg

^{+2}, and Sr

^{+2}. The latter is taken as a trace cation element, and its distribution coefficient is calculated according to: It should be noted that the denominator in this expression is

*Q*. They compared the

*K*values computed with this expression with values measured by Johnston et al. (1985) at an experimental site. Both measured and computed values agreed well and showed that

_{d}*K*(and therefore

_{d}*Q*) remained constant for a wide interval of Sr

^{+2}concentrations, ranging from 10

^{−6}to 1 mg/L. These results provide additional evidence for the validity of our assumption of constant

*Q*for trace cations, satisfying the condition of

*C*≪

_{i}/k_{i}*Q*.

The validity of our assumption has also been tested numerically by solving numerically for nonlinear reactive cation exchange with the reactive transport code CORE^{2D} (Samper et al., 2003; Molinero et al., 2004; Dai and Samper, 2004). The case used to test the assumption corresponds to uncorrelated random conductivity and CEC in a rectangular 40 × 10 m domain. A random realization of conductivity and CEC was generated using Monte-Carlo simulation. Both hydraulic conductivity and CEC had a correlation length of 2.8 m. Variances of log *K* and log CEC were both equal to 1. Na^{+}, K^{+}, and Cl^{−} were instantaneously injected at the center of the domain (see Fig. 1, Animation 1). Initially, the concentration of Cs^{+} was 0.01 mol/L. Other parameters of this numerical simulation can be found in ^{01Table 1}. Figure 1 shows contour plots of the spatial distributions of Na^{+} and K^{+} in the aquifer 45 d after the pulse injection of Na^{+} and K^{+} into the system. For comparison purposes, the plot also shows the spatial distribution of a conservative species, Cl^{−}. It can be seen clearly that the displacement of the center of mass for each species is different. Reactive species are retarded compared to the conservative species, Cl^{−}. K^{+} suffers more retardation than Na^{+} because its selectivity, *k*_{1}, is less than 1.

Temporal fluctuations of *Q* around its initial value, *Q*_{0}, were evaluated at a selected node located downstream from the injection point. Relative fluctuations, computed as 100 (*Q* − *Q*_{0} / Q_{0})%, were < 0.5% (see Fig. 2).

### Effective Retardation Coefficients, Macrodispersivities, and Spatial Moments of Na^{+} and K^{+}

^{+}as the reference cation, retardation factors of other cations, such as K

^{+}, can be written in terms of those of Na

^{+}. From equations 18 and 19, it follows that:

*R̄*is related to statistics of CEC through

_{i}*CV*= [σ

_{Ri}_{CEC}(

*P*/

*k*)]/[1 + (

_{i}*P*/

*k*)] with

_{i}*k*= 1 for

_{i}*i*= 1, which corresponds to Na

^{+}. The right-hand side of this expression can be rewritten in terms of the coefficient of variation of CEC,

*CV*

_{CEC}, and

*R̄*:

_{i}It can be seen that in general *CV _{Ri}* depends on both

*CV*

_{CEC}and

*R̄*. For a strongly sorbing cation (

_{i}*R̄*≫ 1), both coefficients of variation coincide. However, a cation having a retardation factor close to 1 has

_{i}*CV*much smaller than that of CEC.

_{Ri}^{+}or K

^{+}are given by where γ is the flow factor, while the so-called effective retardation factors,

*R*, are given by where

_{C,i}*R̄*and σ

_{i}_{R,i}are given by equations 21 and 22;

*t*=

_{D,i}*vt*/λ

_{1}

*R̄*is dimensionless time; and

_{i}*b*is a parameter that defines the correlation structure between the retardation factor of each cation and log

_{i}*K*. The values of

*b*

_{Na}and

*b*

_{K}can be calculated from equations 11, 21, and 22 by taking into account equation 23, such that

As stated by Miralles-Wilhelm and Gelhar (1996), the field scale retardation factor has a macrokinetic term that is caused by the cross-correlation of retardation and the log hydraulic conductivity field and by the local retardation factor field. It should be noted that for a large value of time, effective retardation factors are equal to the mean retardation factors. Macrodispersivities in equation 25 also include a transient term.

^{+}and K

^{+},

*X̄*, are given by while longitudinal, second-order spatial moments can be computed from:

^{+}to that of K

^{+}at the asymptotic phase is given by:

It can be seen that at the asymptotic phase, this ratio depends on the selectivity coefficient (*k*_{1}), *P*, and the mean CEC.

^{+}to that of K

^{+}at the asymptotic phase is given by

^{+}to that of K

^{+}at the asymptotic phase can be written as:

One can see that this ratio depends on *k*_{1}, the mean CEC, and the correlation structure between log *K* and CEC.

## RESULTS

A synthetic case inspired from the flow configuration and heterogeneity conditions of the Borden site (Miralles-Wilhelm and Gelhar, 1996; Sudicky, 1986; Roberts et al., 1986) is employed here to evaluate the behavior of multicomponent competitive monovalent cation exchange in heterogeneous aquifers by using the stochastic approach outlined in the previous section. Both hydraulic conductivity and CEC are assumed to be random processes having correlation lengths of 2.8 m along the mean direction of groundwater flow. The variance of log *K* is 0.29, and the average velocity of groundwater is 0.1 m/d. CEC has a mean value of 0.1 meq/100 g (typical of a sandy aquifer) and a coefficient of variation of 70%. Selectivity coefficients *k*_{1} and *k*_{2} are equal to 0.2 and 0.08, respectively, which correspond to published values for selectivities of Na^{+} to K^{+} and Na^{+} to Cs^{+} (Appelo and Postma, 1993). Correlation between log *K* and CEC is defined by parameter *b* in equation 12, which takes values of −0.1 (negative correlation), 0 (uncorrelated), and 0.1 (positive correlation). ^{01Table 1} summarizes the main parameters for this synthetic case.

It is worth noting that retardation factors, spatial moments, and macrodispersivities of each cation depend on the variance and correlation length of log *K*, the correlation structure between log *K* and CEC, and selectivity coefficients. Effective retardation factors, macrodispersivities, and spatial moments of each cation during both pre-asymptotic and asymptotic phases are analyzed next.

### Effective Retardation Factors

According to equation 26, the effective retardation factor, *R _{C,i}* for the

*i*th cation shows a transient trend that depends on the coefficient of variation of CEC (

*CV*

_{CEC}), the variance of

*f*(

*b*), and the mean retardation factor (

*R̄*).

_{i}*R*depends on a cation-specific dimensionless time,

_{C,i}*t*, which measures the time needed for a cation to migrate a distance equal to the correlation length of log

_{D,i}*K*.

Figure 3 shows the effective retardation factors of Na^{+} and K^{+} as a function of dimensionless time for two values of *CV*_{CEC} (30% and 70%). One can see that *R _{C,i}* depends on time and reaches the asymptotic mean value, which equals 1.125 for Na

^{+}and 1.61 for K

^{+}. Dimensionless times needed to reach the mean values are different for each cation. Such times are approximately equal to 3 for Na

^{+}and 4 for K

^{+}. As stated by Miralles-Wilhelm and Gelhar (1996), the effective retardation factor is significantly influenced by the correlation structure of log

*K*and retardation. One can see in Figure 3 that for a negative correlation structure (

*b*< 0), pre-asymptotic values of

*R*are smaller than mean values, regardless of the value of

_{C,i}*CV*

_{CEC}. The reverse is true for a positive correlation structure (

*b*> 0).

Curves in Figure 3 illustrate also that pre-asymptotic retardation factors, *R _{C,i}*, for

*CV*

_{CEC}= 30% are always greater than those corresponding to

*CV*

_{CEC}= 70%. Therefore, the smaller the

*CV*

_{CEC}, the larger the effective retardation coefficient. Figure 4 illustrates the sensitivity of effective retardation factors to changes in

*R*, for a negative correlation structure and

_{C,i}*R*values reach large values at

_{C,i}*t*= 0.

### Macrodispersivities

Miralles-Wilhelm et al. (1997) claimed dispersivities of reactive species in heterogeneous systems may be different. Next we show that our results confirm their statement. According to equation 25, cation macrodispersivities show a transient pre-asymptotic phase and depend on *CV*_{CEC},

*b*,

*R*, and λ

_{C}_{1}.

Figure 5 shows the macrodispersivities of Na^{+} and K^{+} together with that of a nonreactive species as a function of dimensionless time (*t _{D,K}*) for positive and negative correlation structures. One can see that macrodispersivities for Na

^{+}and K

^{+}are different from each other and depend on the correlation structure between log

*K*and CEC. Macrodispersivities of all species increase with increasing dimensionless time and asymptotically attain constant values, which are equal to 0.77 m for Na

^{+}and 1.26 m for K

^{+}for a negative correlation structure. Macrodispersivities are smaller for a positive correlation structure: 0.46 m for Na

^{+}and 0.19 m for K

^{+}. It should be noted that the asymptotic macrodispersivity of a conservative species is equal to 0.6 m. Macrodispersivity of the nonreactive species is smaller (greater) than that of reactive species for a negative (positive) correlation structure. The macrodispersivity of Na

^{+}for a negative correlation structure is ∼1.5 times that for a positive correlation structure. A negative correlation structure results in a greater macrodispersivity of reactive species than a positive correlation structure.

Figure 5 also shows the ratio of the macrodispersivities of Na^{+} to K^{+} as a function of dimensionless time for three correlation structures. Such ratios slightly decrease with time for all three correlation structures and reach constant values. The ratio is greater than two for a positive correlation structure but less than one (∼0.7) for a negative correlation structure (see Fig. 5). At the asymptotic phase, the ratio is 0.7 for a negative correlation, around 1 for an uncorrelated structure, and ∼ 2.4 for a positive correlation structure.

Stochastic studies of heterogeneous porous media have focused on the major role played by the variance of log *K*. This parameter measures the degree of physical heterogeneity. To evaluate the effect of this parameter on the behavior of reactive species involved in competitive cation exchange, macrodispersivities at the asymptotic phase were computed for a wide range of values of the variance of log *K* from 0.001 to 1. Figure 6 shows macrodispersivities of Na^{+}, K^{+}, and a nonreactive species as a function of the variance of log *K* at the asymptotic phase for negative and positive correlation structures. Macrodispersivities of all species increase with increasing variance of log *K*. The macrodispersivity of the nonreactive species is greater (smaller) than that of Na^{+} and K^{+} for a positive (negative) correlation structure. For

^{+}and 0.12 m for K

^{+}for a positive correlation structure. They increase to 1.4 m for Na

^{+}and 2.7 m for K

^{+}for a negative correlation structure.

Figure 6 also shows the ratio of the macrodispersivities of Na^{+} to K^{+} as a function of the variance of log *K* at the asymptotic phase. One can see that the ratio is nearly constant and does not depend on the variance of log *K* for uncorrelated and negative correlation structures. From equation 32, one can see that for *b* = 0, the macrodispersivity ratio is equal to 1. For a positive correlation structure, the ratio of the macrodispersivity of Na^{+} to K^{+} increases from 2 to 4.5 when

The coefficient of variation of CEC is the parameter that measures the degree of geochemical heterogeneity. Macrodispersivities are not sensitive to *CV*_{CEC} because the case considered here assumes a perfect relationship between *f* and CEC without a contribution from η, the random residual accounting for the variability of CEC not related to the variability of log *K*. In the case of perfect correlation between *f* and CEC, the coefficient of variation of CEC has almost no effect on the macrodispersivities of Na^{+} and K^{+}. The contribution of η to the macrodispersivity could be significant depending on the effective retardation factor and the correlation length. Future studies should address the case of imperfect correlation between CEC and *f*.

According to equation 25, cation macrodispersivities at the asymptotic phase increase linearly with increasing correlation length.

Cation selectivity coefficients are parameters that control the distribution of mass between solution and solid in cation exchange reactions. From equation 32, it follows that macrodispersivities of Na^{+} and K^{+} are equal when *k*_{1} is 1. Figure 7 shows the macrodispersivities of reactive and nonreactive species as a function of selectivity coefficient, *k*_{1}, at the asymptotic phase. Macrodispersivity of the nonreactive species is 0.6 m. Macrodispersivity of Na^{+} does not depend on *k*_{1} and is equal to 0.46 m for a positive correlation and 0.77 m for a negative correlation. For a negative correlation structure, the macrodispersivity of K^{+} decreases with increasing *k*_{1} and asymptotically approaches the macrodispersivity of a nonreactive species. It is greater than that of Na^{+} if *k*_{1} is less than 1 and smaller otherwise. For a positive correlation structure, it increases with increasing *k*_{1} and asymptotically reaches the macrodispersivity of a nonreactive species.

The effect of *k*_{2} on macrodispersivities is illustrated in Figure 8. One can see that cation macrodispersivities increase with increasing selectivity coefficient, *k*_{2}, for a negative correlation structure. The opposite holds true for a positive correlation structure. Macrodispersivities of Na^{+} and K^{+} reach asymptotically identical values for a negative correlation structure and for a large value of *k*_{2}, but the difference in macrodispersivities of Na^{+} and K^{+} becomes more significant for a positive correlation structure.

### Spatial Moments

According to equations 28 and 29, spatial moments for the *i*th cation depend on *CV*_{CEC},

*b*,

*R*, and λ

_{C}_{1}. The results of our analysis of spatial moments are consistent with published results, which indicate that the first-order spatial moment of a single reactive species does not depend on the correlation structure, while the second-order spatial moment does. A positive correlation structure leads to a decrease in the second-order spatial moment, whereas a negative correlation structure causes an increase in the second-order moment (Bellin et al., 1993; Bosma et al., 1993; Hu et al., 1995).

Reactive species are retarded and move more slowly than a nonreactive species because Na^{+} and K^{+} are exchanged with Cs^{+} and get sorbed onto the solid. For this reason the first-order spatial moment of K^{+} is smaller than that of Na^{+}, which in turn is smaller than that of the conservative species. Figure 9 shows the second-order spatial moments of Na^{+}, K^{+}, and a nonreactive species as a function of dimensionless time. One can see that the second-order spatial moment of a nonreactive species is smaller (greater) than that of reactive species for a negative (positive) correlation structure. Differences in second-order spatial moments of Na^{+} and K^{+} are more significant for a positive correlation structure than for a negative correlation structure.

Figure 10 shows asymptotic first- and second-order spatial moments as a function of selectivity coefficient, *k*_{1}, for both positive and negative correlation structures. First-order spatial moments of reactive species are smaller than those of nonreactive species. However, one can see that the first-order spatial moments of Na^{+} and the nonreactive species do not depend on the selectivity coefficient. The first-order moment of K^{+}, however, increases with increasing selectivity coefficient and reaches a value equal to that of Na^{+} when *k*_{1} = 1. Therefore, the ratio of the first-order spatial moments of Na^{+} and K^{+} increases with increasing *k*_{1} and is equal to 1 when *k*_{1} = 1. It can be seen that the differences in second-order spatial moments between Na^{+} and K^{+} are greater for a positive correlation structure than they are for negative correlation.

## CONCLUSIONS

We have presented here the stochastic analysis of multicomponent competitive monovalent cation exchange. Coupled transport equations for dissolved cations have been decoupled by assuming that the weighted sum of cation concentrations, *Q*, is constant. Derived expressions for cation retardation factors depend on CEC and selectivity coefficients. Randomness of CEC leads to random retardation factors. We have derived analytical expressions for effective retardation factors, longitudinal macrodispersivities, and concentration spatial moments for a chemical system made of Na^{+}, K^{+}, and Cs^{+} using the stochastic analytical solution of Miralles-Wilhelm and Gelhar (1996). The validity of the assumption of constant *Q* has been tested in a synthetic case by numerically solving the reactive transport of the three cations. The maximum relative change in *Q* was found to be <0.5%. Effective retardation factors, spatial moments, and macrodispersivities of Na^{+} and K^{+} and their ratios have been analyzed in terms of the correlation structure of log-conductivity (log *K*) and CEC, the variance and correlation length of log *K* and selectivity coefficients. The analysis of the results leads to the following major conclusions:

Effective retardation factors, spatial moments, and macrodispersivities of K

^{+}differ significantly from those of Na^{+}.Cation effective retardation factors are time dependent and asymptotically attain their mean values. Times needed to reach mean values are different for each cation. Effective retardation factors are influenced by the correlation structure between log

*K*and CEC: the smaller the*CV*_{CEC}value, the larger the effective retardation factor. In the case of perfect correlation between log*K*and CEC, the coefficient of variation of CEC has a minor effect on the macrodispersivities of Na^{+}and K^{+}. The smaller (larger) the variance of log*K*,\({\sigma}\ _{\mathit{f}}^{2}\), the greater (smaller) the effective retardation factor for a negative (positive) correlation.Macrodispersivities of reactive species at the asymptotic phase increase with increasing variance of log

*K*and depend on the correlation structure. The macrodispersivity of a nonreactive species is greater than that of reactive species for positive correlation, while it is smaller than that of reactive species for negative correlation. The macrodispersivity of Na^{+}is smaller than that of K^{+}for a negative correlation structure. The opposite holds true for the positive correlation structure. Macrodispersivity of K^{+}increases with*k*_{1}and decreases with*k*_{2}for the positive correlation structure, while it decreases with*k*_{1}and increases with*k*_{2}for the negative correlation structure.The first-order spatial moment of Na

^{+}is greater than that of K^{+}but smaller than that of the nonreactive species at the asymptotic phase.The second-order spatial moment of nonreactive species is greater than that of a reactive species for a positive correlation structure, but smaller than that for a negative correlation structure. At the asymptotic phase, the second-order spatial moment of Na

^{+}is greater than that of K^{+}for a positive correlation structure. The opposite holds true for a negative correlation structure. Differences in second-order spatial moments of Na^{+}and K^{+}are more significant for a positive correlation structure than for a negative correlation structure.

The method outlined in this paper can be applied to any number of monovalent cation exchange reactions as long as the weighted sum of dissolved cation concentrations, *Q*, is constant. Conditions for the validity of this assumption have been provided. Future work is needed to test the results of the stochastic analysis presented in this paper with both Monte-Carlo simulations and field data.

## APPENDIX A. DERIVATION OF EQUATIONS 14, 15, AND 16

_{Na}in terms of cation concentrations and selectivity coefficients, such that

Substituting A4, A5, and A6 back into equation 8, and taking into account that we are dealing with monovalent cations (*z _{i}* = 1), gives equations 14, 15, and 16.

This work was supported by the Spanish Nuclear Waste Company (ENRESA) through contracts 770045 and 770125, the European Union through projects NF-PRO (FI6W-CT-2003-02389) and FUNMIG (FP6-516514), the Spanish Ministry of Science and Technology (CICYT Project HID98-282), Xunta de Galicia (Incentive PGIDT04PXIC11801PM), and the University of La Coruña through a research scholarship awarded to the second author. We thank the two reviewers for their constructive comments and suggestions, which resulted in an improved version of the paper.