Attribution: You must attribute the work in the manner specified by the author or licensor (but no in any way that suggests that they endorse you or your use of the work).Noncommercial ‒ you may not use this work for commercial purpose.No Derivative works ‒ You may not alter, transform, or build upon this work.Sharing ‒ Individual scientists are hereby granted permission, without fees or further requests to GSA, to use a single figure, a single table, and/or a brief paragraph of text in other subsequent works and to make unlimited photocopies of items in this journal for noncommercial use in classrooms to further education and science.

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, RC,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 (CVCEC). The smaller (larger) the variance of log K,

\({\sigma}\ _{\mathit{f}}^{2}\)
, the greater (smaller) the effective retardation factor for a negative (positive) correlation. Cation macrodispersivities for a positive correlation structure are smaller than that of a nonreactive species and increase with increasing
\({\sigma}\ _{\mathit{f}}^{2}\)
. On the other hand, for a negative correlation structure they are larger than that of a nonreactive species. Macrodispersivity of Na+ is smaller (larger) than that of K+ for a negative (positive) correlation structure. The macrodispersivity of K+ increases with k1 and decreases with the Na+/Cs+ selectivity, k2, 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 = C1 + C2/k1 + C3/k2) is constant, where Ci is dissolved cation concentration, and ki 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+ 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

Transport equations for N dissolved cations undergoing cation exchange in porous media can be expressed as:  
formula
where Ci is concentration of the ith dissolved cation (mol/L); Wi is concentration of the ith exchanged cation in the solid phase (mol/L); D is the local dispersion tensor (L2/T); and v is the water velocity vector (L/T).

Cation Exchange Reactions

Cation exchange takes place when free dissolved cations exchange with interlayer cations of clay minerals. This process can be described as an equilibrium reaction between dissolved and exchanged cations. According to the Gaines-Thomas convention (Appelo and Postma, 1993), cation exchange reactions of Na+, K+, and Cs+ can be expressed as  
formula
 
formula
 
formula
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:  
formula
 
formula
where k1 and k2 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:  
formula
Equivalent fractions are related to exchanged cation concentrations in equation 1, Wi, through  
formula
where CEC is cation exchange capacity (meq/100 g of solids), 𝛉 is porosity (−), zi is charge (in this case, zi = 1 ∀ i), ρd is dry density of medium (kg solid/dm3), and i = Na+, K+, and Cs+.

Random Fields

Concentrations of dissolved cations, Ci, and exchanged cations, Wi, as well as groundwater velocity, 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:  
formula
The log conductivity field f(x) can be expressed as the sum of its ensemble mean and a random perturbation f′:  
formula
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
\({\sigma}\ _{\mathit{f}}^{2}\)
(Gelhar, 1993).
Several researchers (Cassel et al., 2000; Jacques et al., 1999; Kuzyakova et al., 2001; Matschonat and Vogt, 1996; Pucci et al., 1997; Wirth, 2001) have reported that CEC of soils can be regarded as a stochastic field showing distinct spatial correlation as well as cross-correlation with log K according to  
formula
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.
The mean,graphic, and fluctuation, CEC′, of CEC are given by  
formula
 
formula

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.

Concentrations of exchanged species (Wi) 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, WNa, WK, and WCs can be written as  
formula
 
formula
 
formula
where  
formula

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+.

By substituting equations 14 and 15 back into equation 1, and assuming that Q is constant, one finds  
formula
where L(Ci) = ∇ · (DCi) – vCi, and  
formula
 
formula

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. RNa and RK are retardation factors for Na+ and K+, respectively, while CEC/100Q and CEC/100Qk1 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.

The validity of constant 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:  
formula
It should be noted that the denominator in this expression is Q. They compared the Kd 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 Kd (and therefore 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 Ci /kiQ.

The validity of our assumption has also been tested numerically by solving numerically for nonlinear reactive cation exchange with the reactive transport code CORE2D (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, k1, is less than 1.

Temporal fluctuations of Q around its initial value, Q0, were evaluated at a selected node located downstream from the injection point. Relative fluctuations, computed as 100 (QQ0 / Q0)%, were < 0.5% (see Fig. 2).

Effective Retardation Coefficients, Macrodispersivities, and Spatial Moments of Na+ and K+

Taking Na+ 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:  
formula
According to the definition of cation retardation factors given in equations 18 and 19, randomness of CEC causes retardation factors to become random. Therefore, the mean and variance of cation retardation factors can be computed from those of CEC  
formula
and  
formula
where graphicand
\({\sigma}\ _{CEC}^{2}\)
are the mean and variance of CEC, respectively; and  
formula
From equations 21 and 22, it follows that the coefficient of variation of i is related to statistics of CEC through CVRi = [σCEC(P/ki)]/[1 + (P/ki)graphic] with ki = 1 for 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, CVCEC, and i:  
formula

It can be seen that in general CVRi depends on both CVCEC and i. For a strongly sorbing cation (i ≫ 1), both coefficients of variation coincide. However, a cation having a retardation factor close to 1 has CVRi much smaller than that of CEC.

According to the results of Miralles-Wilhelm and Gelhar (1996), the macrodispersivities of Na+ or K+ are given by  
formula
where γ is the flow factor, while the so-called effective retardation factors, RC,i, are given by  
formula
where i and σR,i are given by equations 21 and 22; tD,i = vt1i is dimensionless time; and bi is a parameter that defines the correlation structure between the retardation factor of each cation and log K. The values of bNa and bK can be calculated from equations 11, 21, and 22 by taking into account equation 23, such that  
formula

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.

First-order spatial moments of Na+ and K+, , are given by  
formula
while longitudinal, second-order spatial moments can be computed from:  
formula
The ratio of the first-order spatial moment of Na+ to that of K+ at the asymptotic phase is given by:  
formula

It can be seen that at the asymptotic phase, this ratio depends on the selectivity coefficient (k1), P, and the mean CEC.

The ratio of the second-order spatial moment of Na+ to that of K+ at the asymptotic phase is given by  
formula
Similarly, the ratio of the macrodispersivity of Na+ to that of K+ at the asymptotic phase can be written as:  
formula

One can see that this ratio depends on k1, 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 k1 and k2 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, RC,i for the ith cation shows a transient trend that depends on the coefficient of variation of CEC (CVCEC), the variance of f (

\({\sigma}\ _{\mathit{f}}^{2}\)
), the correlation structure (b), and the mean retardation factor (i). RC,i depends on a cation-specific dimensionless time, tD,i, which measures the time needed for a cation to migrate a distance equal to the correlation length of log K.

Figure 3 shows the effective retardation factors of Na+ and K+ as a function of dimensionless time for two values of CVCEC (30% and 70%). One can see that RC,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 RC,i are smaller than mean values, regardless of the value of CVCEC. The reverse is true for a positive correlation structure (b > 0).

Curves in Figure 3 illustrate also that pre-asymptotic retardation factors, RC,i, for CVCEC = 30% are always greater than those corresponding to CVCEC = 70%. Therefore, the smaller the CVCEC, the larger the effective retardation coefficient. Figure 4 illustrates the sensitivity of effective retardation factors to changes in

\({\sigma}\ _{\mathit{f}}^{2}\)
. During the pre-asymptotic phase, retardation factors, RC,i, for a negative correlation structure and
\({\sigma}\ _{\mathit{f}}^{2}\)
= 0.29 are greater than those corresponding to
\({\sigma}\ _{\mathit{f}}^{2}\)
= 1. Therefore, the smaller the value of
\({\sigma}\ _{\mathit{f}}^{2}\)
, the greater the effective retardation factor. However, the opposite holds true for a positive correlation structure. For a positive correlation structure and a large
\({\sigma}\ _{\mathit{f}}^{2}\)
, all RC,i values reach large values at 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 CVCEC,

\({\sigma}\ _{\mathit{f}}^{2}\)
, b, RC, and λ1.

Figure 5 shows the macrodispersivities of Na+ and K+ together with that of a nonreactive species as a function of dimensionless time (tD,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

\({\sigma}\ _{\mathit{f}}^{2}\)
= 1, macrodispersivities are equal to 0.35 m for Na+ 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

\({\sigma}\ _{\mathit{f}}^{2}\)
increases from 0.001 to 1.

The coefficient of variation of CEC is the parameter that measures the degree of geochemical heterogeneity. Macrodispersivities are not sensitive to CVCEC 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 k1 is 1. Figure 7 shows the macrodispersivities of reactive and nonreactive species as a function of selectivity coefficient, k1, at the asymptotic phase. Macrodispersivity of the nonreactive species is 0.6 m. Macrodispersivity of Na+ does not depend on k1 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 k1 and asymptotically approaches the macrodispersivity of a nonreactive species. It is greater than that of Na+ if k1 is less than 1 and smaller otherwise. For a positive correlation structure, it increases with increasing k1 and asymptotically reaches the macrodispersivity of a nonreactive species.

The effect of k2 on macrodispersivities is illustrated in Figure 8. One can see that cation macrodispersivities increase with increasing selectivity coefficient, k2, 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 k2, 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 ith cation depend on CVCEC,

\({\sigma}\ _{\mathit{f}}^{2}\)
, b, RC, and λ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, k1, 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 k1 = 1. Therefore, the ratio of the first-order spatial moments of Na+ and K+ increases with increasing k1 and is equal to 1 when k1 = 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:

  1. Effective retardation factors, spatial moments, and macrodispersivities of K+ differ significantly from those of Na+.

  2. 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 CVCEC 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.

  3. 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 k1 and decreases with k2 for the positive correlation structure, while it decreases with k1 and increases with k2 for the negative correlation structure.

  4. 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.

  5. 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

Equation 5 can be rewritten so that βk is on the left-hand side:  
formula
Similarly, from equation 6, one can get the expression for βCs:  
formula
By substituting A1 and A2 back into equation 7, one obtains:  
formula
This first-order equation can be solved for βNa in terms of cation concentrations and selectivity coefficients, such that  
formula
Substitution of A4 into A1 and A2 provides explicit expressions for βk and βCs, respectively:  
formula
 
formula

Substituting A4, A5, and A6 back into equation 8, and taking into account that we are dealing with monovalent cations (zi = 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.