- Soil Science Society of America
Migration of radionuclides under the SX-tank farm at the Hanford nuclear waste complex involves interaction of variably water saturated sediments with concentrated NaOH–NaNO3–NaNO2 solutions that have leaked from the tanks. Constant Kd 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+, Ca2+, 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.
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 + NaNO3 + NaNO2, 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 137Cs 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.
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 1with the motility of the kth species denoted by mk with valence zk, the extended Debye–Hückel activity coefficient γk defined by 2is 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 3where γk refers to the individual ion activity coefficient (Pitzer, 1981; Felmy, 1995), and γ0k denotes a modified form of the Debye–Hückel activity coefficient. The expansion coefficients, Bkk′(I) and Ckk′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 4where R denotes the gas constant, T is temperature, Ww 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 5and 6for the ith complex or gaseous species Aπi, where species Aj 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. ) 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 Clj from the law of mass action as 7with equilibrium constant Kπi associated with the ith secondary species and πth phase.
Mineral reactions have the similar form 8expressed in terms of primary species, with stochiometric reaction matrix νjm associated with mineral Mm. Mineral reactions involve mass transfer between the aqueous and solid phases and their rates are described through a kinetic rate law of the form 9provided by transition state theory. In this expression, km denotes the kinetic rate constant, am the mineral surface area per bulk volume of porous medium, Km refers to the equilibrium constant for reaction Eq. , and Qm denotes the ion activity product defined by 10The quantity in parentheses in Eq.  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.  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.
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 Ca2+ 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+–Ca2+ 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 NaNO3. Zachara et al. (2002) and Liu et al. (2003a) used NaNO3 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 NaNO3 electrolyte concentrations, two effects control retardation of Cs+: competitive cation exchange and aqueous complexing reaction with NO−3. A set of competitive exchange reactions involving Cs+, Na+, K+, and Ca2+ may be written as the half reactions 11a11b11c11dwhere 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 12More generally, for two cations Azj+j and Azk+k with valencies zj, zk, the overall exchange reaction has the form 13The corresponding mass action equation has the form 14where 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 15and γ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 16where 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 17The exchange site concentration, in units of eq/bulk dm3, is related to the sediment cation exchange capacity Qσ (mol kg−1) by the expression 18where φ 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 19with the local distribution coefficient KDj defined as 20where Ψlj denotes the total concentration of the jth primary species (see Eq. ). 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.  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 K̃Dj, with units of (mol/g−1 solid)/(mol/cm−3 liquid), by the expression 21
Ignoring heterovalent Ca2+ exchange and considering only monovalent exchange of Cs+, Na+, and K+, the local distribution coefficient has the explicit form 22In the Pitzer formulation of activity coefficients, if it is assumed that no explicit complexes form, then . In this case the effects of complexing are included in the activity coefficients γlj, 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 Ψlj > Clj. 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 NaNO3 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 × 10−5, 2.62 × 10−4, and 0.12 mol kg−1 are used. The species CsNO3, 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−3 (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 (NaNO3(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, KDCs+, 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.
For the 20 m Na+ case, the figure also shows the distribution coefficient calculated at 65°C using Steefel et al. (2003) selectivity coefficients determined at 25°C with selectivity coefficients measured by Liu et al. (2003b) at 65°C. The cation exchange capacity corresponding to the first and third sites are used in Liu et al. (2003b) case. The difference in distribution coefficients between the 25 and 65°C calculations using the same selectivity coefficients from Steefel et al. (2003) is a result of changes in speciation of NaNO3(aq) and Cs+ activity coefficient with temperature. The distribution coefficient calculated using Liu et al. (2003b) selectivity coefficients coincides with the distribution coefficient at 25°C using the Steefel et al. (2003) coefficients for Cs+ concentrations in the range 2 × 10−4 to 2 × 10−3 mol/L−1, but is somewhat reduced at lower and increased at higher Cs+ concentrations.
HEAT AND MASS CONSERVATION EQUATIONS
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 23a23b23c23dand 23eproviding a total of Np + Nmin + 3 partial differential equations for Np primary species, Nmin minerals, pg, sg, and pa, 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 24where ksat denotes water saturated permeability, krπ(sl) 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 cr 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 25sπ refers to liquid water saturation of phase π = l, g with 26Liquid and gas pressures are related by the capillary pressure pc(sl): 27where the capillary pressure is a function of water saturation. Vapor pressure lowering is incorporated through the Kelvin equation 28where pv denotes the vapor pressure, and psat(T) refers to the saturation pressure of pure water. The effective gaseous diffusion coefficient has the form 29where D0g refers to the diffusion coefficient in a pure gas phase, λg denotes an enhancement factor, and pref and Tref refer to reference pressure and temperature. The quantities Sw,a, Se, and Sj 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 30and 31respectively, where the sum over i runs over secondary species, with the solute flux Jπk for individual aqueous and gaseous species defined by 32This 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 V̅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 33For 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 NaNO3 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 NaNO3 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.
NUMERICAL SIMULATIONS OF CESIUM TRANSPORT IN THE HANFORD VADOSE ZONE
The SX tank farm consists of 15 single shell tanks with different heat loads and chemical compositions, many of which have leaked radioactive fluid into the vadose zone. Tanks were emplaced at a depth of 15.4 m below ground surface. Tank height is 13.6 m, with a radius of 11.8 m. Spacing between tank centers is 30.4 m. After emplacement the tanks were backfilled with excavated sediment and topped with a gravel layer 1.8 m thick. Because of the finite number of tanks, their variable heat loads, and the varying positions of suspected leaks from the tanks, it is not possible in a rigorous sense to reduce the size of the problem by symmetry considerations. Nevertheless, to gain insight into the processes affecting Cs+ mobility within the vadose zone, simplifications are necessary to reduce the system to a manageable size for which relatively rapid simulations can be performed that capture the basic physics and chemistry involved. Use of a 1D model is clearly a gross simplification compared to two- or three-dimensional (2D or 3D) models and therefore should not be considered quantitatively accurate. In particular, a 1D model is unable to account for lateral spreading of the contaminant plume. This is a particularly important property of the Hanford sediments caused by capillary barriers due to heterogeneous layering of fine- and coarse-grained sediment (Pruess et al., 2002). As a consequence, for example, the leak volume and duration used in a 1D model must be obtained by calibration to field data. It is also important to keep in mind that it may not be reliable to extrapolate the 1D model to future times far beyond the calibration period. However, because of the vastly reduced computational effort required for a 1D model compared with 2D or 3D models, the 1D model allows sensitivity analyses to be performed more easily with incorporation of multiphase, multicomponent chemistry. Clearly this problem would be ideally suited to high performance computing techniques to be carried out in future work. Alternatively, such simulations may be considered more generally as streamtubes taken from a 3D flow field (although no attempt is made to do this here). In this case the spatial coordinate represents the distance along the streamtube. However, steady-state conditions cannot be assumed because of the pulse release of fluid from the tank.
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 × 105 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 Dl = 10−9 m2/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 × 10−5 m2 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.
Initial flow conditions and moisture profile are determined by first running the model without the tank in place. A steady-state profile for moisture content, temperature, and infiltration is obtained which is then used with the tank in place. The initial steady-state water saturation profile is shown in Fig. 3 along with profiles at indicated times following release of fluid from the tank using the material properties and van Genuchten parameters listed in Table 4. After approximately 10 yr the moisture profile has returned to its steady-state value. Two zones of elevated water saturation, located at a depth of approximately 43 to 48 m, coincide with the Plio-Pleistocene and Upper Ringold Gravel stratigraphic units. The water table is located at a depth of approximately 62.4 m. In these simulations the tank temperature is fixed at 25°C.
Calibration to Well 299-W23-19
Extensive data for concentration of solute species with depth are available from Well 299-W23-19, located in the vicinity of the SX-115 tank (Serne et al., 2001a). The timing of the leak for SX-115 tank is thought to have occurred in 1965, and leak volume and duration are estimated to be 60000 gallons for a 2-wk period, based on historical inventory estimates (Jones et al., 2000). The position of the leak is thought to be close to the edge of the tank near the well. Because the geometry of the system is inherently 3D, it is not possible to use the actual leak volumes and rates in the 1D simulation. To obtain these parameters for the 1D model, it is necessary to calibrate the model against field observations.
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+, Ca2+, Mg2+, NO−3, Cl−, CO2(g), and H+. Reactions include several aqueous complexing reactions; cation exchange reactions involving Na+, K+, Ca2+, and Mg2+; and precipitation and dissolution of calcite using a kinetic rate law of the form given in Eq. . For exchange of Mg2+, the selectivity coefficients were used as fit parameters because data are not available for this species. The resulting fit values are only approximate. Mg2+ is more weakly sorbed compared than Ca2+. 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.
From the observed profile for NO−3, a conservative species that does not take part in exchange, several parameters can be determined. The width of the NO−3 profile is used to fix the leak duration, and the depth of the peak concentration determines both the infiltration and leak rate. A time span of 35 yr is used corresponding to the time of the observations. For the leak rate a value of 8.4 × 10−4 kg s−1 is obtained, and for the leak duration a value of 1.213 × 106 s, or just more than 2 wk, in agreement with the estimated value. With these values a total leak volume of 247.6 gallons is released in the 1D model. It is difficult to compare the fitted value with the much larger value of 60000 gallons believed to have leaked from the SX-115 tank (Jones et al., 2000). In the model calculation a cross-sectional area of 1 m2 is used. To match the estimated leak of 60000 gallons, a cross-sectional area of 60000/247.6 ≃ 242 m2 would be needed. While the area of the leak at the tank is not known, the leak would not necessarily need to occur from a single point source. The cross-sectional area of the tank itself is π(11.8 m)2 = 437.4 m2, or 1.8 times larger than the area needed to match the estimated leak volume in the model calculation. In addition, lateral spreading is neglected in the 1D model.
An infiltration rate of ql = 0.08 m yr−1 gave the best fit to the NO3 profile, consistent with Gee et al. (1992). With this value the approximate travel time to the water table (t ≃ ls̅lφ̅/ql, l = depth to water table, s̅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 Dnum = qlΔz/(2slφ) ≃ 10−9 m2 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−3 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 Na2 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 Ca2+. Using isotopic exchange of 22Na, Steefel et al. (2003) obtained a value of 0.046 mol kg−1 with a 2 mM NaNO3 solution in a batch experiment, but obtained a larger value of 0.12 mol kg−1 with a 1 m NaNO3 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 NaNO3 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 (cm2 cm−3), was fit to approximately match the pH profile observed in Well 299-W23-19 resulting in a fit value of 5 × 10−13 mol cm−3 s−1.
The results of the fit to the Serne et al. (2001a) date from Well 299-W23-19 are shown in Fig. 4 and 5 . Concentration profiles for species NO−3, Na+, Ca2+, and pH are shown in Fig. 4 along with field observations from Well 299-W23-19 (Serne, 2001a). While generally speaking the numerical simulation appears to capture the main features of the field data there are some differences. The simulation does not succeed in describing the broad peak in the Na+ concentration profile, nor does it capture the smaller peak at greater depth. The upper portion of the NO−3 is not well described either. There is a small displacement between the Na+ and NO−3 peaks, indicating that Na+ is slightly retarded compared to NO−3 which acts like a tracer. Note that the peak concentrations are less than the source term values by almost a factor of two for NO−3 resulting from dilution effects, and roughly a factor of 7.5 for Na+ resulting from sorption of Na+ on mineral surfaces as well as dilution. The pH profile was rather remarkably reproduced. The sharp rise in pH (at 21 m) was caused by dissolution of calcite at the leading edge of the zone where secondary calcite formed. The pH gradually decreases as calcite precipitates and begins to increase again as the reaction slows. Calcite disappears and stops reacting at a depth of approximately 47.5 m.
Figure 5 shows comparisons of calculated and observed profiles for K+ and Mg2+. For both of these species, there is a drop in concentration in the region of high Na+, followed by an increase in concentration further downstream produced by displacement of the cations from the sediment by Na+. The Mg2+ profile was fit by adjusting the initial and injected Mg2+ concentrations, and setting the Mg2+ selectivity coefficients so that Mg2+ is more weakly sorbed compared with Ca2+. The K+ profile was fit by adjusting only the initial and injected K+ concentrations. The resulting fit values are listed in Tables 5 and 6.
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.  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 34For Cs+–Na+ exchange, this ratio ranges from 107.25 to 101.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−3. 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−3 concentration would tend to increase Cs+ retardation by increasing its aqueous activity coefficient (equivalent to reducing the concentration of the aqueous complex CsNO3). 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 137Cs+ also was not taken into account in the calculations. It should be noted that the total Cs inventory contains the isotopes 133Cs and 135Cs, in addition to 137Cs.
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.
The relation between the aqueous and sorbed Cs+ concentration profiles and Na+ is noteworthy. As the Na+ source concentration increases, the sorbed Cs+ concentration penetrates more deeply into the vadose zone. Near the base of the tank close to the source of the leak, Cs+ is highly retarded once the Na+ pulse advances beyond the injection point. For this reason, Cs+ is not completely displaced from its point of release. The aqueous Cs+ concentration is rapidly attenuated with depth as more Cs+ is sorbed on the sediments after the Na+ plume has passed. The sorbed Cs+ profile formed a quasistationary state as time increased and remained essentially fixed in the vadose zone in the presence of dilute ambient water.
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+ concentration from the slant borehole (Serne et al., 2001b). The calculated sorbed Cs+ profile shows remarkable agreement with the observed profile, especially considering that no additional adjustments were made to the model other than increasing the ionic strength of the tank leak and increasing the leak duration to 1 yr. The measured 137Cs+ concentration is corrected to include the isotopes 133, 135, and 137 to give the total Cs+ concentration using the relation Cs+total = 4.75(±0.18) 137Cs+ obtained by Evans et al. (2001) for the total Cs+ inventory in terms of 137Cs+. The depth of penetration of the Cs+ plume is reasonably well reproduced; however, the peak concentration is somewhat too high. Probably the close agreement should not be taken too seriously, given the various approximations made in the simulations. Slight changes in the cation exchange capacity, effect of elevated temperatures, and the extrapolation of exchange data to high Na+ concentration could certainly alter the results. Nevertheless they suggest that the model could be capturing the gross behavior of the system. It should be noted that the model does not capture the deeper, smaller peak of Cs+ observed at a depth of roughly 40 m. This could be a result of several causes, such as two separate leak events with the second leak displacing Cs+ emplaced during the first leak, contamination of the borehole during drilling, lateral spreading of the Cs+ plume, or contamination from neighboring tanks (Serne et al., 2001b). It is also quite possible that Cs+ could have migrated deeper than would be apparent from the slant borehole observations (Serne et al., 2001b).
To confirm agreement with the sorbed Cs+ in the slant borehole, it would be necessary to carry out a detailed comparison with the full suite of species concentrations that were measured including NO−3 and Na+, similar to the comparison made for Well 299-W23-19 associated with the SX-115 tank. Such a detailed comparison is beyond the scope of the present contribution and is reserved for a future study.
To investigate the possibility of advancing the Cs+ further, two additional runs were performed for the 20 m Na+ case. First, the effect of elevated temperature was investigated using the selectivity coefficients given in Liu et al. (2003b) at 70°C for exchange of Na+ and Cs+. The result was only a slight advancement of the Cs+ front by several meters. Second, the effect of increased K+ concentration in the tank leak was investigated. An unrealistic value of 0.5 m K+ resulted in advancement of the Cs+ front again by only several meters. The effect of nonideality of the sorbed concentration was not investigated, but presumably this effect would tend to increase Cs+ retardation (Liu et al., 2003a). Short of dramatically reducing the sediment cation exchange capacity, it was not possible to increase the penetration of Cs+ deeper into the vadose zone by any appreciable amount.
CONCLUSION AND IMPLICATIONS
By utilizing field data from Well 299-W23-19 in the neighborhood of the SX-115 tank, it was possible to calibrate flow parameters and source term representing a tank leak for a 1D reactive flow and transport model. The model was then used to analyze the mobility of Cs+ with varying Na+ concentrations in the background electrolyte solution associated with the SX-108 tank fluid. It was found that only for very high Na+ concentrations (>5 m) was it possible for Cs+ to advance to greater depths in the vadose zone. However, the aqueous Cs+ concentration was rapidly attenuated with depth, as more Cs+ exchanged with the sediments following passage of the Na+ plume. The sorbed Cs+ profile became essentially fixed in the vadose zone with very high retardation in the presence of dilute ambient water. The case with 20 m Na+, consistent with the estimated SX-108/SX-109 tank supernatant (Lichtner and Felmy, 2003), concentration gave the best fit to field observations at the slant borehole for the sorbed Cs+ concentration. If correct, this result suggests that Cs+ contamination did not advance significantly farther than what was observed in the slant borehole. Contrary to what one might expect intuitively, dilute infiltrating rainwater, rather than enhancing migration of Cs+, actually serves to reduce Cs+ mobility by creating conditions favorable for strong Cs+ retardation. With a half-life of 30.2 yr, roughly 300 yr are required for 137Cs+ to decay to negligible levels, which appears to be more than adequate for 137Cs+ to decay to harmless levels within the vadose zone.
It would appear, given the close agreement between the predicted Cs+ retardation and the observed behavior of Cs+ in the slant borehole, that a useful exercise would be to combine the complexity of the nonisothermal flow model presented in Pruess et al. (2002) extended to 3D with the high ionic strength and multicomponent cation exchange model considered here. This would allow a fully mechanistic approach to be applied to the SX tank farm and enable a more complete understanding to be obtained of the detailed mechanisms involved in radionuclide mobility. The new understanding could be applied to other sites at Hanford and elsewhere within the DOE complex.
åk, Debye radius (Å)
aw, Activity of water
𝒜gi, Designation for ith gaseous species in gas phase (g)
𝒜li, Designation for ith secondary species in liquid water (l)
𝒜j, Designation for jth primary species
A(T), Debye-Hückel parameter [(kg H2O/mol)1/2]
am, Specific surface area of mth mineral (cm−1)
ḃ(T), Debye b-dot parameter (kg H2O mol−1)
B(T), Debye–Hückel parameter [(kg H2O/mol)1/2/Å]
ℬkk′(I), Pitzer model expansion coefficient (kg H2O mol−1)
cr, Bulk rock capacity (kJ mol−1 K−1)
𝒜kk′k″, Pitzer model expansion coefficient [(kg H2O/mol)2]
Cπi, Concentration of ith secondary species in πth phase (mol L−1)
Clj, Concentration of jth primary species in aqueous phase (mol L−1)
Dπ, Diffusion–dispersion coefficient associated with πth fluid phase (m2 s−1)
εσj, Sorbed equivalent mole fraction of jth cation
g, Acceleration of gravity (m s−2)
Hw, Enthalpy of pure water (kJ mol−1)
Hπ, Enthalpy of fluid phase π (kJ mol−1)
I, Ionic strength of an aqueous solution (mol kg H2O−1)
Im, Kinetic reaction rate for mth mineral (mol dm−3 s−1)
Jπi, Flux of ith secondary species in phase π (mol m−2 s−1)
Jπj, Flux of jth primary species in phase π (mol m−2 s−1)
ksat, Water saturated permeability (m2)
krπ (sl), Relative permeability of πth phase
kσj, Exchange half-reaction selectivity coefficient for jth cation
km, Kinetic rate constant for mth mineral (mol m−2 s−1)
kπi, Equilibrium constant for ith secondary species reaction in πth phase
Kσjk, Selectivity coefficient for exchange of the jth and kth cations on sites σ
KDj, Local dimensionless distribution coefficient for jth cation
K̃Dj, Local dimensioned distribution coefficient for jth cation (L g−1)
Km, Equilibrium constant for mth mineral reaction
mk, molality of kth species [mol (kg H2O)−1]
ℳm, Designation for mth mineral
nπ, Molar density of πth fluid phase (mol m−3)
Nm, Number of minerals
Np, Number of primary species
pa, Air partial pressure (Pa)
pc, Capillary pressure (Pa)
pg, Gas pressure (Pa)
Pl, Liquid pressure (Pa)
Pv, Water vapor partial pressure (Pa)
ql, Liquid Darcy flux (m s−1)
qg, Gas Darcy flux (m s−1)
Qm, Ion-activity product associated with mth mineral
Qπ, Cation exchange capacity associated with sites σ [eq kg−1]
r, Position vector (m)
ℛ, Constant (J K−1 mol−1)
ℛj, Local retardation factor of jth cation
sl, Liquid saturation
sg, Gas saturation
Sσj, Sorbed concentration of jth cation on exchange sites σ (mol dm−3)
t, Time (s)
τπ, Tortuosity of fluid phase π
T, Temperature (°C)
T0, Reference temperature (°C)
Uw, Internal energy of pure water (kJ mol−1)
Uπ, Internal energy of fluid phase π (kJ mol−1)
V̅m, Molar volume of mth mineral (cm3 mol−1)
Ww, Formula weight of water (g mol−1)
Xσ, Cation exchange site with valence −1
Xπk, Mole fraction of kth species in phase π
zk, Valence of kth species
α, van Genuchten capillary parameter (m−1)
γi, Activity coefficient of ith secondary species
γj, Activity coefficient of jth primary species
κ, Thermal conductivity (W m−1)
λ, van Genuchten parameter
λj, Sorbed activity coefficient of jth cation
μw, Viscosity of pure water (Pa s)
νlji, Stoichiometric coefficient matrix for homogeneous reactions
νgji, Stoichiometric coefficient matrix for gaseous reactions
νjrn, Stoichiometric coefficient matrix for mineral reactions
ωσ, Concentration of exchange sites designated by σ (mol dm−3)
φm, Volume fraction of mth mineral
ρr, Rock density (kg m−3)
ρs, Sediment grain density (kg m−3)
ρπ, Density of πth fluid phase (kg m−3)
Ψπj, Total concentration of jth primary species in phase π (mol L−1)
Ωπj, Total flux of jth primary species in phase π (mol m−2 s−1)
Ø, Osmotic coefficient of water
Super- and Subscripts
i, Secondary species
j, Primary species
k,k′,k″, Arbitrary aqueous species
l, Liquid phase
g, Gas phase
σ, Designation for exchange type
π, Designation of πth phase
The authors would like to thank two anonymous reviewers for thorough reviews of an earlier version of the manuscript that led to significant improvements to its final form. In addition we would like to thank John Zachara for helpful comments that improved the manuscript. This work was supported by the Hanford Science and Technology program managed by the Groundwater/Vadose Zone Integration Project.
- Received June 16, 2003.