HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in ISI Web of Science Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (18) Citing Articles via Google Scholar Google Scholar Articles by Lichtner, P. C. Articles by Steefel, C. I. Search for Related Content PubMed Articles by Lichtner, P. C. Articles by Steefel, C. I. GeoRef GeoRef Citation Agricola Articles by Lichtner, P. C. Articles by Steefel, C. I. Related Collections Vadose Zone Processes and Chemical Transport Multicomponent Transport Models

ORIGINAL RESEARCH PAPERS

Role of Competitive Cation Exchange on Chromatographic Displacement of Cesium in the Vadose Zone beneath the Hanford S/SX Tank Farm

^a Los Alamos National Laboratory, SM-30 Bikini Atoll Rd., Los Alamos, NM 87545
^b Pacific Northwest National Laboratory, P.O. Box 999, Richland, WA 99352
^c Lawrence Berkeley National Laboratory, 1 Cyclotron Rd., MS 90-1116, Berkeley, CA 94720
^d Lawrence Livermore National Laboratory, 7000 East Ave., Livermore, CA 94550-9234
^* Corresponding author (lichtner{at}lanl.gov).

Received 16 June 2003.

	ABSTRACT

TOP ABSTRACT INTRODUCTION CHEMICAL PROCESSES HEAT AND MASS CONSERVATION... NUMERICAL SIMULATIONS OF CESIUM... CONCLUSION AND IMPLICATIONS APPENDIX REFERENCES

 
Migration of radionuclides under the SX-tank farm at the Hanford nuclear waste complex involves interaction of variably water saturated sediments with concentrated NaOH–NaNO₃–NaNO₂ solutions that have leaked from the tanks. Constant K_d models for describing radionuclide retardation are not valid under these conditions because of strong competition for sorption sites by abundant Na⁺ ions, and because of dramatically changing solution compositions with time as the highly concentrated tank fluid becomes diluted as it mixes with infiltrating rainwater. A mechanistic multicomponent sorption model is required that can account for effects of competition and spatially and temporally variable solution compositions. To investigate the influence of the high ionic strength tank fluids on Cs⁺ migration, numerical calculations are performed using the multiphase-multicomponent reactive transport code FLOTRAN. The computer model describes reactive transport in nonisothermal, variably saturated porous media including both liquid and gas phases. Pitzer activity coefficient corrections are used to describe the high ionic strength solutions. The calculations take into account multicomponent cation exchange based on measured selectivity coefficients specific to the Hanford sediments. Solution composition data obtained from Well 299-W23-19, documenting a moderately concentrated leak from the SX-115 tank, are used to calibrate the model. In addition to exchange of cations Na⁺, K⁺, Ca²⁺, and Cs⁺, aqueous complexing and a kinetic description of precipitation and dissolution of calcite are also included in the calculations. The fitted infiltration rate of 0.08 m yr^–1, and fitted cation exchange capacity of 0.05 mol kg^–1 are consistent with measured values for the Hanford sediments. A sensitivity analysis is performed for Na⁺ concentrations ranging from 5 to 20 m to investigate the mobility of Cs⁺ interacting with a highly concentrated background electrolyte solution believed to have been released from the SX-108/SX-109 tanks. The calculations indicate that during the initial period of the tank leak when Cs⁺ is associated with high Na⁺ concentrations, there is little retardation of the Cs⁺ plume. However, as time increases the Na⁺ and Cs⁺ profiles become chromatographically separated due to differences in their selectivity coefficients and dilution of the tank leak plume with infiltrating rainwater. Eventually the two species become separated spatially, and Cs⁺ becomes highly retarded and remains essentially fixed in the sediments by cation exchange. For the 20 m Na⁺ simulated tank leak, the sorbed Cs⁺ profile is in close agreement with data obtained from the slant borehole and consistent with the estimated tank supernatant concentration. The simulations suggest that natural attenuation processes should result in strong fixation of Cs⁺ in the vadose zone in spite of the release of high Na⁺ concentrations during a tank leak event.

Abbreviations: FES, frayed edge sites • 1D, one-dimensional • 2D, two-dimensional • 3D, three-dimensional

	INTRODUCTION

TOP ABSTRACT INTRODUCTION CHEMICAL PROCESSES HEAT AND MASS CONSERVATION... NUMERICAL SIMULATIONS OF CESIUM... CONCLUSION AND IMPLICATIONS APPENDIX REFERENCES

 
THE HANFORD TANK FARM COMPLEX is situated on the Columbia River plateau in a semiarid region in south-central Washington. The site served as a Pu production facility for nuclear weapons from 1944 to the end of the Cold War era. Storage tanks at the site contain waste fluid streams from the Pu extraction process. Numerous leaks have occurred at the Hanford Site from single-shelled tanks containing high-level radioactive waste (it has been estimated that 67 of 149 single-shell tanks have leaked). The largest areas of vadose zone contamination occur in the S-SX Waste Management Area located in the 200 Area of the Hanford site, consisting of the 241-S and 241-SX tank farms, near Tanks S-104, SX-108/SX-109, and SX-115. Several of these tanks contain high ionic strength fluids with the general composition NaOH + NaNO₃ + NaNO₂, which because of the high Na⁺ concentrations can adversely affect the migration of radionuclides by reducing their retardation. Competitive cation exchange effects, involving a limited number of exchange sites, between radionuclides (e.g., Cs⁺) and cations (primarily Na⁺) in the background electrolyte solution leaked from the tanks are expected to enhance radionuclide mobility. This study focuses on the SX tank farm and in particular tanks SX-115 and SX-108/SX-109. The purpose of this paper is to investigate the role high ionic strength fluids might have played in enhancing the migration of Cs⁺ leaked from the tanks. A companion study (Pruess et al., 2002) investigated the unique hydrologic system created by the tanks. We focus on competitive cation exchange effects between Cs⁺ and Na⁺ as well as other minor cations that may lead to enhanced mobility during and following a tank leak event.

Leaked ¹³⁷Cs has been observed in boreholes at the SX tank farm at depths greater than anticipated based on previous estimates that relied on transport models representing Cs⁺ mobility through a single distribution coefficient and hence neglecting competitive cation exchange effects on retardation (Serne et al., 2001b). The depth of penetration of Cs⁺ appears to be different depending on the tank waste composition. The SX-115 tank shows little migration of Cs⁺, whereas for the SX-108/SX-109 tanks Cs⁺ appears to have migrated several tens of meters beneath the tanks. This is consistent with the observation, if competitive cation exchange is the cause for the more rapid migration, that fluid leaked from the SX-115 tank is much more dilute compared with the SX-108/SX-109 tanks.

In order to better understand the potential hazard Cs⁺ could pose to the environment and also to facilitate remediation efforts to clean up the site, it is important to understand the mechanisms controlling the mobility of Cs⁺ in the Hanford vadose zone. In the environment of the SX tank farm, the concentrations of both Cs⁺ and Na⁺ were highly variable with time and depth within the vadose zone. Leaks from the tanks took place over relatively short periods of time, releasing short duration pulses (weeks to months or years) of variable composition tank effluents into the vadose zone. Under such conditions, the mobility of Cs⁺ also would have been highly variable. During a pulse release with high Na⁺ concentration, Cs⁺ mobility would be greatly enhanced. However, even under such conditions a slight difference in the Na⁺ and Cs⁺ distribution coefficients would exist, with Cs⁺ still slightly more retarded than Na⁺. Following cessation of the leak and dilution by infiltrating rainwater, chromatographic separation between the Na⁺ and Cs⁺ plumes would increase. Eventually, given sufficient time and travel distance, the two plumes would become completely separated, with Cs⁺ left behind in a dilute solution with a greatly reduced mobility. The picture that emerges from these qualitative considerations is that during the initial release of fluid from the tank, Cs⁺ would migrate rapidly as long as it was in the presence of high Na⁺ concentrations. But after the leak ceases and chromatographic separation sets in, Cs⁺ would become highly retarded. How far Cs⁺ could actually migrate through the vadose zone would appear to depend on a number of factors including the leak composition, rate, and duration, surface infiltration rate, and sediment sorptive properties.

We use the two-phase, reactive flow and transport model FLOTRAN (Lichtner, 2001) to describe the interaction of the leaked tank fluid with host sediments and ambient water in the vandose zone. The effect of high ionic strength fluid, ranging from 1 to 20 m Na⁺ concentration, leaked from tanks on the mobility of Cs⁺ through competitive cation exchange is investigated. A simplified, one-dimensional (1D) model is considered. The model incorporates flow of liquid water, water vapor, air, and heat, in addition to transport of multicomponent reactive solutes. The Pitzer activity coefficient algorithm for high ionic strength fluids is used in the model calculations.

	CHEMICAL PROCESSES

TOP ABSTRACT INTRODUCTION CHEMICAL PROCESSES HEAT AND MASS CONSERVATION... NUMERICAL SIMULATIONS OF CESIUM... CONCLUSION AND IMPLICATIONS APPENDIX REFERENCES

 
Conservation equations describing transport of radionuclides leaked from storage tanks at the Hanford waste complex must account for the high ionic strength tank fluids. For the most concentrated tanks SX-108 and SX-109, these fluids are estimated to have high pH values on the order of 14, and high Na⁺ concentrations on the order of 20 m (Lichtner and Felmy, 2003). Under such circumstances it might be expected that Cs⁺, an exchangeable cation, could be highly mobile, with Na⁺ easily displacing Cs⁺.

Chemical reactions involved in interaction of the tank leak with ambient vadose zone water and sediments include various homogeneous complexing reactions taking place within the aqueous phase, mineral precipitation and dissolution reactions, and sorption reactions, specifically cation exchange in the case of retardation of Cs⁺. Laboratory studies of the interaction of highly basic tank fluids with S-SX sediments indicate that significant mineral alteration and precipitation (e.g., of zeolites) occurs as a result of tank fluid–sediment interaction (McKinley et al., 2001). This result led to speculation that a large, deep mineral alteration areole should exist surrounding the leaked tanks that would impact contaminant behavior. However, significant sediment alteration was not observed in samples collected from the slant borehole (Serne et al., 2001b). The slant borehole was drilled at an approximate angle of 30° from the vertical, passing directly beneath the center of the SX-108 tank to a depth of 43.9 m (total borehole length is 52.2 m) (Serne et al., 2001b). X-ray diffraction measurements indicated no evidence of secondary mineral formation, although "faint indications of caustic alteration" were found approximately 9 m beneath the base of the SX-108 tank (Serne et al., 2001b).

In addition, the laboratory experiments indicated that the high base-treated sediment had a surprisingly small impact on Cs⁺ exchange with the altered sediment (Zachara et al., 2001). Thus, although the mineralogy of the clay fraction changed, this had minimal effect on the overall sorptivity of Cs⁺. Thus, apparently, the high base reaction did not significantly affect the high affinity frayed edge sites associated with micas. The residence time of the leaked fluid in the sediment before it becomes diluted by mixing with ambient water in the vadose zone and infiltrating water may be too short for significant alteration of aluminosilicate minerals to occur, in spite of its caustic nature (Serne et al., 2001b).

Pitzer Activity Coefficient Algorithm
For relatively low ionic strength fluids, defined as I 0.1 m, where I denotes the ionic strength of the fluid given by

[1]

with the motility of the kth species denoted by m_k with valence z_k, the extended Debye–Hückel activity coefficient k defined by

[2]

is adequate for computing activity coefficient corrections. In this expression A(T), B(T), and (T) refer to temperature-dependent coefficients, and å_k denotes the Debye radius. The activity of the solvent, in this case water, is unity for an ideal dilute solution. In the Debye–Hückel algorithm, the activity coefficient is a function of ionic strength I and is the same for species with identical valencies and Debye radii å_k.

However, for the high ionic strength fluids of primary interest in this study, the Debye–Hückel model is inadequate. For such fluids an approach such as the Pitzer model is needed (Pitzer, 1981). In this approach the activity coefficients are expressed in a virial expansion of the form

[3]

where _k refers to the individual ion activity coefficient (Pitzer, 1981; Felmy, 1995), and ⁰_k denotes a modified form of the Debye–Hückel activity coefficient. The expansion coefficients, B_kk'(I) and C_kk'k'', must be determined through fits to experimental measurements for a range of pressure and temperature conditions. The activity of water is calculated from the relation

[4]

where R denotes the gas constant, T is temperature, W_w refers to the formula weight of water, and Ø denotes the osmotic coefficient of water (see Pitzer, 1981).

Aqueous Complexing and Mineral Reactions
Homogeneous reactions taking place within the aqueous phase and heterogeneous aqueous–gaseous reactions can be written in the general form

[5]

and

[6]

for the ith complex or gaseous species A_i, where species A_j refers to a primary species, and _ji denotes the stoichiometric reaction matrix. Subscripts i and j are reserved for primary and secondary species, respectively, and subscript k is reserved for either primary or secondary species in the following. Primary species, selected from the set of aqueous species, serve as independent basis species to write the reactions. Their choice is arbitrary, and primary and aqueous secondary species may be freely interchanged so long as the resulting reactions (Eq. [5]) remain linearly independent (Lichtner et al., 1996). Concentrations of secondary species in local chemical equilibrium are obtained as nonlinear functions of the concentrations of primary species C^l_j from the law of mass action as

[7]

with equilibrium constant K_i associated with the ith secondary species and th phase.

Mineral reactions have the similar form

[8]

expressed in terms of primary species, with stochiometric reaction matrix _jm associated with mineral M_m. Mineral reactions involve mass transfer between the aqueous and solid phases and their rates are described through a kinetic rate law of the form

[9]

provided by transition state theory. In this expression, k_m denotes the kinetic rate constant, a_m the mineral surface area per bulk volume of porous medium, K_m refers to the equilibrium constant for reaction Eq. [8], and Q_m denotes the ion activity product defined by

[10]

The quantity in parentheses in Eq. [9] is referred to as the affinity factor and vanishes at equilibrium. The sign of the affinity factor determines whether precipitation or dissolution occurs, with a positive rate designating precipitation corresponding to reaction Eq. [8] proceeding from left to right. If a mineral is not present at some particular point in space, it may precipitate if it becomes supersaturated. In this case it is necessary to provide an initial surface area associated with the mineral. Obviously, a mineral cannot dissolve if it is not present, and the rate in this case is zero.

Cation Exchange
Retardation of Cs⁺ by the Hanford sediments involves a number of distinct exchange sites (Zachara et al., 2002; Steefel et al., 2003). In these studies multiple exchange sites were proposed to capture the observed behavior of multicomponent exchange reactions involving cations Cs⁺, Na⁺, K⁺, and Ca²⁺ on Hanford sediments. Zachara et al. (2002) proposed two distinct sites to represent batch experiments for Cs⁺–Na⁺ exchange and up to five sites to describe Cs⁺–Ca²⁺ exchange. Steefel et al. (2003), on the basis of batch and column experiments, proposed three distinct sites corresponding to two frayed edge sites (FES) on weathered micas, which exhibit high affinity for Cs⁺ compared with Na⁺, and one planar site associated with expandable clays with lower affinity for Cs⁺ compared with Na⁺, and one planar site associated with expandable clays with lower affinity for Cs⁺. The planar sites account for more than 99% of the cation exchange capacity of bulk Hanford sediment.

In addition to these studies, more recent work has considered the role of nonideality of the sorbed concentration (Liu et al., 2003a), temperature dependence of selectivity coefficients (Liu et al., 2003b), and desorption kinetics (Liu et al., 2003c). Steefel et al. (2003) also considered solid phase nonideality effects but used an empirical relation, in contrast to the more rigorous approach taken by Liu et al. (2003a).

None of the studies mentioned above considered Na⁺ concentrations at the levels expected in leaks from the SX-108/SX-109 tanks, estimated to be on the order of 20 m (Lichtner and Felmy, 2003). Steefel et al. (2003) considered Na⁺ concentrations in the range 0.01 to 5 M NaNO₃. Zachara et al. (2002) and Liu et al. (2003a) used NaNO₃ ranging from 0.01 to 7 M, well below the expected leak concentrations. In this study, the effects of nonideality, temperature, and desorption kinetics on exchange are neglected for simplicity and because they have yet to be extended to the ionic strengths relevant for this study. To some extent the effects of nonideality and temperature may cancel one another (Liu et al., 2003a, 2003b), with increasing temperature leading to a reduction and increasing ionic strength an increase in retardation.

At high NaNO₃ electrolyte concentrations, two effects control retardation of Cs⁺: competitive cation exchange and aqueous complexing reaction with NO^–₃. A set of competitive exchange reactions involving Cs⁺, Na⁺, K⁺, and Ca²⁺ may be written as the half reactions

[11a]

[11b]

[11c]

[11d]

where the subscript refers to a particular exchange site denoted by X. Each half reaction is associated with a selectivity coefficient k_j referring to the jth cation, assumed to be a primary species for convenience (although this is not essential). The overall exchange reaction (e.g., exchange of Na⁺ and Cs⁺) is obtained by combining two half reactions balancing on the exchange site X^– to give

[12]

More generally, for two cations A^zj+_j and A^zk+_k with valencies z_j, z_k, the overall exchange reaction has the form

[13]

The corresponding mass action equation has the form

[14]

where K_ji denotes the selectivity coefficient for the overall exchange reaction related to the half-reaction selectivity coefficients k_j, k_k, by the expression

[15]

and _j, _k, and _j, _k, refer to aqueous and solid activity coefficients, respectively. For the Gaines–Thomas formulation the quantities _j refer to equivalent mole fractions defined as

[16]

where the sorbed concentration S_j for the jth cation on site is referenced to bulk volume of sediment, and where refers to the concentration of exchange sites related to the sorbed concentrations S_j by the expression

[17]

The exchange site concentration, in units of eq/bulk dm³, is related to the sediment cation exchange capacity Q (mol kg^–1) by the expression

[18]

where denotes the porosity of the sediments and the sediment particle density is denoted by _s. It is tacitly assumed that the site concentration is independent of the water saturation state of the porous medium. This assumption is based on the suppositions that liquid water is the wetting phase and therefore even for partially saturated conditions, mineral surfaces are assumed to be wetted, and all exchange sites are in contact with liquid water.

The local retardation factor _j(r, t) at a point in space r and time t is defined by

[19]

with the local distribution coefficient K^D_j defined as

[20]

where ^l_j denotes the total concentration of the jth primary species (see Eq. [30]). The distribution coefficient gives the ratio of the sorbed to aqueous phase concentration averaged over a control volume. The distribution coefficient as defined in Eq. [20] is dimensionless. It is, in general, not constant but varies spatially and temporally, and is directly proportional to the local sorbed concentration and inversely proportional to the local water saturation state, porosity, and total aqueous concentration. It is related to the more conventionally defined distribution coefficient ^D_j, with units of (mol/g^–1 solid)/(mol/cm^–3 liquid), by the expression

[21]

Ignoring heterovalent Ca²⁺ exchange and considering only monovalent exchange of Cs⁺, Na⁺, and K⁺, the local distribution coefficient has the explicit form

[22]

In the Pitzer formulation of activity coefficients, if it is assumed that no explicit complexes form, then ^l_j = C^l_j. In this case the effects of complexing are included in the activity coefficients ^l_j, with a value <1 indicating attractive ion interactions. This is manifested in a reduction in the distribution coefficient, similar to the effect of explicitly including complexes in which case ^l_j > C^l_j. The same considerations apply to heterovalent exchange.

To investigate the mobility of Cs⁺ on solution composition, a simplified calculation of the distribution coefficient for exchange of cations Na⁺ and Cs⁺ is considered in a background NaNO₃ electrolyte solution using the three-site sorption model developed for the Hanford sediments by Steefel et al. (2003) at 25°C. In addition, the effect of temperature on the distribution coefficient is shown comparing the Steefel et al. (2003) model with the two-site model presented by Liu et al. (2003b) using selectivity coefficients measured at 65°C for the 20 m Na⁺ concentration. Cation exchange capacities reported by Steefel et al. (2003) of 2.285 x 10^–5, 2.62 x 10^–4, and 0.12 mol kg^–1 are used. The species CsNO₃, considered by Steefel et al. (2003) in fitting the Cs⁺ exchange isotherm, is not included explicitly in this study because the Pitzer model already includes binary interaction parameters for Cs⁺ and NO^–₃ (Pitzer, 1991). However, available data are limited to relatively low ionic strengths (Pitzer, 1991). Also not included is addition of a term proportional to ionic strength in the sorbed activity coefficients (Steefel et al., 2003). Here ideality is assumed for sorbed concentrations. At 20 m Na⁺, nitratine (NaNO_3(s)) is slightly supersaturated at 25°C, and close to equilibrium at 65°C. The selectivity coefficients for Na⁺ and Cs⁺ used in the calculations are listed in Table 2. Fig. 1 shows the Cs⁺ distribution coefficient, K^D_Cs+, plotted against the concentrations of Cs⁺. A porosity of 30% and grain density of 2.65 g cm^–3 is used in obtaining the retardation. As can be seen from the figure, the Cs⁺ distribution coefficient ranges over six orders of magnitude. Cesium becomes more strongly retarded as its concentration decreases and as the concentration of Na⁺ decreases. However, as the concentration of Na⁺ increases to high values, Cs⁺ retardation rapidly decreases. The contributions from the three different exchange sites are visible in the figure as plateaus separated by transition regions. The high affinity, low abundance, frayed edge sites give the major contribution to the distribution coefficient at low Cs⁺, followed by the medium affinity frayed edge sites, and low affinity planar sites.

View larger version (28K): [in this window] [in a new window]   Fig. 1. Distribution coefficient for Cs+ exchanging with Na+ on Hanford sediments plotted as a function of aqueous Cs+ concentration for different Na+ concentrations indicated in the figure using the Steefel et al. (2003) selectivity coefficients at 25°C listed in Table 2. Also shown for the 20 M Na+ case is the distribution coefficient calculated at 65°C but using Steefel et al. (2003) selectivity coefficients determined at 25°C (up triangle), and compared with selectivity coefficients measured by Liu et al. (2003b) at 65°C (down triangle).

	HEAT AND MASS CONSERVATION EQUATIONS

TOP ABSTRACT INTRODUCTION CHEMICAL PROCESSES HEAT AND MASS CONSERVATION... NUMERICAL SIMULATIONS OF CESIUM... CONCLUSION AND IMPLICATIONS APPENDIX REFERENCES

 
Calculations are performed using the reactive flow and transport computer code FLOTRAN (Lichtner, 2001), a two-phase, multicomponent, reactive flow and transport numerical model that describes nonisothermal, flow of mass and heat in variably water saturated porous media, and (enhanced) binary diffusion of water vapor and air. The flow equations are coupled to equations describing multicomponent reactive transport through field variables pressure, temperature, liquid and gas velocities, and water saturation state. The reactive transport equations account for such processes as aqueous speciation, reactions with minerals, cation exchange and surface complexation, and colloid-facilitated transport. The transport equations, in turn, can be coupled to the flow equations through changes in material properties such as porosity, permeability, tortuosity, and fluid density, although this possibility is not considered in the simulations carried out below.

In terms of the general forms of homogeneous and heterogeneous reactions described above, the computer code FLOTRAN (Lichtner, 2001) solves the following transient, two-phase, mass conservation equations for the solvent water (w), air (a), heat, a set of primary solute species, and minerals

[23a]

[23b]

[23c]

[23d]

and

[23e]

providing a total of N_p + N_min + 3 partial differential equations for N_p primary species, N_min minerals, p_g, s_g, and p_a, the partial pressure of air. Temperature is obtained from the condition of thermodynamic equilibrium in a two-phase system. In these equations the Darcy fluxes for liquid and gas phases q, ( = l, g) are given by

[24]

where k_sat denotes water saturated permeability, k^r(s_l) denotes the relative permeability of the th phase, in general a function of water saturation, g denotes the acceleration of gravity, and denotes the mass density and µ the viscosity of phase . The quantities , D, n, U = H – p/n, and H denote the tortuosity, diffusion/dispersion coefficient, molar density, internal energy, and enthalpy, respectively, of phase , _r and c_r refer to the rock density and heat capacity, respectively, and X_k refers to the mole fractions of water and air (k = 1, 2) in phase satisfying

[25]

s refers to liquid water saturation of phase = l, g with

[26]

Liquid and gas pressures are related by the capillary pressure p_c(s_l):

[27]

where the capillary pressure is a function of water saturation. Vapor pressure lowering is incorporated through the Kelvin equation

[28]

where p_v denotes the vapor pressure, and p_sat(T) refers to the saturation pressure of pure water. The effective gaseous diffusion coefficient has the form

[29]

where D⁰_g refers to the diffusion coefficient in a pure gas phase, _g denotes an enhancement factor, and p_ref and T_ref refer to reference pressure and temperature. The quantities S_w,a, S_e, and S_j represent source–sink terms accounting for the tank leak, and providing release of mass and heat.

The quantities _j, _j denote the total concentration and flux, respectively, in liquid and gas phases defined by

[30]

and

[31]

respectively, where the sum over i runs over secondary species, with the solute flux J_k for individual aqueous and gaseous species defined by

[32]

This formulation of aqueous diffusion based on a single diffusion coefficient is highly simplistic for high ionic strength electrolyte solutions (Felmy and Weare, 1991a, 1991b). However, a more complex description is beyond the scope of this study and would not be expected to change the results significantly over the relatively short times involved and large length scales of tens of meters.

The mineral mass transfer equation, Eq. [23e] is expressed in terms of the mineral volume fraction _m. The quantity _m denotes the mineral molar volume. Assuming that minerals form with zero intrinsic porosity, the total porosity is related to the mineral volume fractions by the expression

[33]

For the relatively short time spans of tens of years involved, significant changes in porosity are not expected.

The mass and energy conservation equations are solved sequentially. First the water, air, and heat equations are solved over a single time step using a fully implicit backward Euler algorithm. This is followed by solving the solute transport equations over, generally, a smaller time step and interpolating the velocity, pressure, temperature, and water saturation fields passed to the transport equations, at the intermediate times. Again a fully implicit backward Euler algorithm is used. Finally, the mineral mass transfer equations are solved using an explicit finite difference procedure using the kinetic reaction rates for minerals obtained from the solution to the transport equations. This latter simplification is possible because of the slow rates of reaction associated with solids resulting in the approximate decoupling of the mineral and solute conservation equations over a single time step.

It is clear that a substantial increase in the fluid density can occur at high solute concentrations. However, this does not necessarily imply that the fluid moves faster as a result. The fluid mobility, given as the ratio of density to viscosity, determines the rate of movement in response to a pressure gradient. Using data taken from Isono (1984) for concentrated NaNO₃ solutions for density and viscosity, the mobility ratio relative to pure water (/µ)_l/(/µ)_w is plotted in Fig. 2 as a function of temperature. As can be seen from the figure, the mobility decreases with increasing NaNO₃ concentration. Although not apparent from the figure, this is a consequence of the increase in viscosity with increased concentration overwhelming the corresponding increase in density, thereby resulting in a net reduction in mobility.

View larger version (19K): [in this window] [in a new window]   Fig. 2. Mobility ratio (/µ)l/(/µ) for 0.1, 1, 4, and 8 m NaNO3 solutions compared with pure water plotted as a function of temperature. Points refer to data taken from Isono (1984).

	NUMERICAL SIMULATIONS OF CESIUM TRANSPORT IN THE HANFORD VADOSE ZONE

TOP ABSTRACT INTRODUCTION CHEMICAL PROCESSES HEAT AND MASS CONSERVATION... NUMERICAL SIMULATIONS OF CESIUM... CONCLUSION AND IMPLICATIONS APPENDIX REFERENCES

In what follows 1D simulations are carried out along a vertical column from the ground surface to the water table. In the simulations, 136 nodes in the vertical direction are used with equal spacing of 0.5 m. The temperature is held fixed at the top and bottom of the column with values 12.8 and 17.0°C, respectively. A steady-state infiltration rate is imposed at the top of the column. The infiltration rate is used as a fit parameter as discussed below. At the bottom of the column, a constant pressure of 1.4 x 10⁵ Pa is imposed, with fully saturated conditions representing the water table. Fluid is injected at a depth corresponding to the base of the tank. Combined with steady surface infiltration, the interaction of fluid with sediments is simulated to describe the migration of Cs⁺ as it is released from the tank. To obtain an appropriate leak rate for the 1D simulations, the model is calibrated against observations from Well 266-W23-19, which is located next to the SX-115 tank.

The stratigraphic sequence used in the calculations is listed in Table 3. We did not include the possible presence of a thin compacted layer beneath the tanks in these simulations. The compacted layer has been estimated to be approximately 0.5 to 1 m thick. However, the physical properties of the compact layer (porosity, permeability, capillary parameters, and density) are presently unknown. Material properties for soil, hydraulic, and thermal parameters used in the calculations are listed in Table 4. Parameters and n refer to the phenomenological van Genuchten equation relating capillary pressure and water saturation (van Genuchten, 1980). An aqueous diffusion coefficient of D_l = 10^–9 m²/s^–1 and unit tortuosity are used in the simulations. Dispersion was not included in the simulations (but see below for an estimate of numerical dispersion). Effective binary gas diffusion is incorporated with a gaseous diffusion coefficient of 2.13 x 10^–5 m² s^–1 and a tortuosity of 0.5. Vapor pressure lowering resulting from capillary suction was included in the simulations as described by Kelvin's equation. An even more important effect, but not yet incorporated into the simulations, could be the effect of salts on vapor pressure.

View this table: [in this window] [in a new window]   Table 4. Parameters for material and thermal properties for intrinsic rock density s, heat capacity c, thermal conductivity , porosity , residual water saturation sr, van Genuchten parameters and , and vertical water saturated permeability sat. Data taken from Khaleel and Freeman (1995), Khaleel et al. (2001), and Pruess et al. (2002).

View larger version (24K): [in this window] [in a new window]   Fig. 3. Transient and steady-state water saturation profiles resulting from a tank leak with volume and duration given in the text with a steady surface infiltration rate of 0.08 m yr4. Indicated in the figure are the positions of the different stratigraphic units listed in Table 3.

As suggested here, the model parameters for leak rate and duration, average infiltration rate, solution composition of the leaking fluid, and sediment cation exchange capacity are rather tightly constrained by field observations. A stepwise procedure is used to fit the various model parameters, rather than attempting to fit all parameters simultaneously. A simplified chemical system is used, compared with estimated tank compositions, consisting of eight primary species Na⁺, K⁺, Ca²⁺, Mg²⁺, NO^–₃, Cl^–, CO_2(g), and H⁺. Reactions include several aqueous complexing reactions; cation exchange reactions involving Na⁺, K⁺, Ca²⁺, and Mg²⁺; and precipitation and dissolution of calcite using a kinetic rate law of the form given in Eq. [9]. For exchange of Mg²⁺, the selectivity coefficients were used as fit parameters because data are not available for this species. The resulting fit values are only approximate. Mg²⁺ is more weakly sorbed compared than Ca²⁺. Because the major effect of ion interaction is captured in the Pitzer model through activity coefficient corrections to species activities, only a few select aqueous complexes need be considered explicitly as listed in Table 1.

An infiltration rate of q_l = 0.08 m yr^–1 gave the best fit to the NO₃ profile, consistent with Gee et al. (1992). With this value the approximate travel time to the water table (t l_l/q_l, l = depth to water table, _l, = spatially averaged water saturation and porosity) for a conservative tracer released from the surface without the tank present is found to be 100 yr, and 68 yr when released from a depth coinciding with the base of the tank in the presence of the leak. Based on a grid spacing z = 0.5 m, an infiltration rate of 0.08 m yr^–1 introduces numerical dispersion of approximately D_num = q_lz/(2s_l) 10^–9 m² s^–1, comparable to diffusion.

In addition to the leak duration and volume and infiltration rate, the solution composition of the leak can be constrained by field data derived from Serne et al. (2001a). The concentration of NO^–₃ in the leak is adjusted to match the observed peak concentration. Likewise the pH of the leak composition is adjusted to match the observed pH peak at depth. The Na² concentration is determined by charge balance.

The cation exchange capacity is adjusted to approximately fit the peak in the observed Na⁺ profile. If the exchange capacity becomes too large, more Na⁺ is sorbed and the resulting peak in Na⁺ in solution becomes too low. An optimal value of 0.05 mol kg^–1 is found. This value is consistent with values obtained by Zachara et al. (2002) for Hanford sediments. Zachara et al. (2002) found a range of values depending on the saturating ion with values 0.0426 mol kg^–1 for Na⁺, 0.0825 mol kg^–1 for K⁺, and 0.0469 mol kg^–1 for Ca²⁺. Using isotopic exchange of ²²Na, Steefel et al. (2003) obtained a value of 0.046 mol kg^–1 with a 2 mM NaNO₃ solution in a batch experiment, but obtained a larger value of 0.12 mol kg^–1 with a 1 m NaNO₃ solution in a column experiment. The initial and source fluid compositions are given in Tables 5 and 6, respectively. The source term consists of a moderately concentrated 0.75 m NaNO₃ solution. Calcite is reported present in minor amounts of both Well 299-W23-19 and the slant borehole (Serne et al., 2001a, 2001b), and is taken to have an initial volume fraction of 0.01 in the calculations. The calcite effective rate constant, defined as the product of the rate constant in mol cm^–2 s^–1 times the specific surface (cm² cm^–3), was fit to approximately match the pH profile observed in Well 299-W23-19 resulting in a fit value of 5 x 10^–13 mol cm^–3 s^–1.

View larger version (31K): [in this window] [in a new window]   Fig. 4. Calculated solute and pH profiles for an elapsed time of 35 yr compared with field observations (symbols) taken from Well 299-W23-19 at the SX-115 tank (Serne et al., 2001a).

View larger version (23K): [in this window] [in a new window]   Fig. 5. Calculated Mg2+ and K+ profiles for an elapsed time of 35 yr compared with field observations (symbols) taken from Well 299-W23-19 at the SX-115 tank (Serne et al., 2001a).

Effect of Background Electrolyte Concentration on Cesium Mobility
Cation exchange processes are characterized by chromatographic separation between the exchanging cations with distance and time. The extent of separation depends on the affinity of each cation for the exchange surface, the exchange capacity of the medium, and length of time and travel distance involved. In certain cases it is possible to predict the mobility of each cation based on its distribution coefficient, defined as the ratio of sorbed to aqueous concentration. However, in order to use the distribution coefficient approach to estimate mobility, it is necessary that the bulk solution composition remain relatively constant with distance and time. This requirement is most certainly not met at Hanford for situations in which high ionic strength fluids have leaked from the tanks, as is evident from inspection of Fig. 1. As the background electrolyte solution becomes more and more dilute, the Cs⁺ distribution coefficient increases markedly.

For the simplifying case of a dominating single exchange site and assuming aqueous complexing effects are included in the Pitzer activity coefficients , according to Eq. [22] the ratio of distribution coefficients corresponding to different cations is independent of the concentration of either cation and simply equal to the ratio of their respective selectivity coefficients

[34]

For Cs⁺–Na⁺ exchange, this ratio ranges from 10^7.25 to 10^1.99, depending on the exchange site, according to the values listed in Table 2. Thus, even for the case involving high Na⁺ concentrations, the ratio of distribution coefficients is the same, implying that a plume containing several cations will eventually separate into several plumes corresponding to each cation, in relation to their relative selectivity coefficients.

To investigate the effect of high Na⁺ concentrations on the mobility of Cs⁺, varying source term concentrations of Na⁺ were used. The solutions are charged balanced with NO^–₃. For the 20 m Na⁺ solution, nitratine is slightly supersaturated. However, at the elevated temperatures of the SX-108/SX-109 tanks, nitratine becomes close to equilibrium or undersaturated. As noted in Fig. 1, the distribution coefficient for the 20 m Na⁺ solution at 25°C lies closer to the 65°C solution using the Liu et al. (2003b) selectivity coefficients. Reducing the NO^–₃ concentration would tend to increase Cs⁺ retardation by increasing its aqueous activity coefficient (equivalent to reducing the concentration of the aqueous complex CsNO₃). The injection interval was set to 1 yr corresponding roughly to the SX-108/SX-109 tank leaks (Jones et al., 2000), keeping the total mass of fluid injected the same as in the calibration problem. All other variables were kept fixed to the values used in the calibration problem. No attempt was made to mimic the actual temperature distribution in the 1D simulation. However, it was necessary to speciate the tank fluid at 75°C in order to obtain a reasonable pH and chemical composition. The temperature of the leaked fluid dissipated rapidly following injection. Radioactive decay of ¹³⁷Cs⁺ also was not taken into account in the calculations. It should be noted that the total Cs inventory contains the isotopes ¹³³Cs and ¹³⁵Cs, in addition to ¹³⁷Cs.

Results of the simulations are presented in Fig. 6 through 10 , showing aqueous and sorbed Cs⁺ concentrations and aqueous Na⁺ concentrations plotted as a function of depth for 35, 50, and 75 yr. Sodium concentrations ranged from 0.75 m, corresponding to the calibration problem, to 20 m, which could have leaked from SX-108/SX-109 tanks (Lichtner and Felmy, 2003). As is apparent from the figures, with increasing Na⁺ concentration, Cs⁺ mobility increases. For lower concentrations of Na⁺, however, Cs⁺ is more retarded compared with Na⁺, resulting in rapid chromatographic separation of the two ions. As a consequence, Cs⁺ is left behind in a dilute solution in which it is highly retarded. This is apparent in Fig. 6, 7, and 8. However, for more concentrated solutions, Cs⁺ penetrates more deeply into the vadose zone as can be seen in Fig. 9 and 10. Even in these cases Cs⁺ lags behind Na⁺ and is more strongly sorbed with depth.

View larger version (24K): [in this window] [in a new window]   Fig. 6. Calculated total Na+, Cs+, and sorbed SCs+ profiles as a function of depth at the indicated times ranging from 1 to 75 yr for a 0.75 m NaNO3 source concentration corresponding to Tank SK-115.

View larger version (29K): [in this window] [in a new window]   Fig. 10. Calculated total Na+, Cs+, and sorbed SCs+ profiles as a function of depth at the indicated times for 20 m NaNO3 source concentration. The symbols refer to measured Cs+ concentrations from the slant borehole (Serne et al., 2001b).

View larger version (25K): [in this window] [in a new window]   Fig. 7. Calculated total Na+, Cs+, and sorbed SCs+ profiles as a function of depth at the indicated times for 5 m NaNO3 source concentration.

View larger version (26K): [in this window] [in a new window]   Fig. 8. Calculated total Na+, Cs+, and sorbed SCs+ profiles as a function of depth at the indicated times for 10 m NaNO3 source concentration.

View larger version (27K): [in this window] [in a new window]   Fig. 9. Calculated total Na+, Cs+, and sorbed SCs+ profiles as a function of depth at the indicated times for 15 m NaNO3 source concentration.

Initially, the Cs⁺ pulse released from the tank migrates downward. However, as time increases the pulse appears to migrate upward against the downward flow. This happens as the Na⁺ pulse becomes separated from the Cs⁺. As this occurs Cs⁺ becomes more strongly sorbed, and more Cs⁺ is taken up by the solid decreasing the amount of Cs⁺ in solution, thereby giving the appearance that Cs⁺ is moving against the direction of flow. For the 20-m simulation, the aqueous Cs⁺ concentration ceases to change after approximately 35 yr. The sorbed Cs⁺ concentration remains stationary after only roughly 10 yr have elapsed.

Also plotted in Fig. 10 is the observed Cs^+</