MOLECULAR DYNAMICS SIMULATIONS OF ANION EXCLUSION IN CLAY INTERLAYER NANOPORES Open Access

How does water confinement in nanopores impact anion exclusion? The answer to this question has important implications for the migration of salts and groundwater through clay-rich rocks (Fritz, 1986; Leroy et al., 2006), the design of engineered nanofluidic systems (Schoch et al., 2008; Bocquet and Charlaix, 2010), the performance of desalination membranes (Richards et al., 2012), and the transport of ions through channels in cell membranes (Nikaido and Rosenberg, 1981; Nikaido, 2003). It also is of importance to high-level radioactive waste (HLRW) management, where the performance of geologic repositories is determined, in many cases, by the migration of anionic, long-lived fission products $(129I−,99TcO4−,36Cl−)$ through nanoporous clay barriers (Tournassat et al., 2007; Altmann, 2008; Claret et al., 2010; Glaus et al., 2011; Altmann et al., 2012).

At present, extant models of anion adsorption and diffusion in nanoporous clayey media postulate either that anions are completely excluded from clay interlayer nanopores (Sposito, 1992; Wersin et al., 2004; Bourg et al., 2008; Churakov and Gimmi, 2011; Tournassat and Appelo, 2011) or that anions are partially excluded from these nanopores in a manner that can be predicted accurately by assuming nanopore water is the same as bulk liquid water, albeit with a lower electrostatic potential imposed by negatively charged clay surfaces (Leroy et al., 2006; Birgersson and Karnland, 2009; Glaus et al., 2013; Tachi et al., 2014). Efforts to test these two competing hypotheses using anion exclusion measurements have been inconclusive because the results of such measurements depend not only on the mechanistic details of anion exclusion in individual nanopores, but also on the pore-size distribution in saturated, compacted clayey media (Tournassat and Appelo, 2011; Tournassat et al., 2016; Tinnacher et al., 2016), which is not yet fully understood (Holmboe et al., 2012).

Molecular dynamics (MD) simulations offer a way out of this conundrum by directly examining the details of anion exclusion in individual clay interlayer nanopores. However, three recent studies that quantified anion exclusion in clay interlayers by MD simulation have yielded contradictory results: one concluded that interlayer water can be treated as having the same aqueous chemistry as bulk liquid water (Hedström and Karnland, 2012), whereas two others reached the opposite conclusion, that anion exclusion in clay interlayers is much greater than would be expected based on the aqueous chemistry of bulk liquid water (Rotenberg et al., 2007b; Hsiao and Hedström, 2015). Consequently, the correct way in which to account for anion exclusion in macroscopic models of ion migration in clay formations and engineered clay barriers remains uncertain (Muurinen et al., 2007; Birgersson and Karnland, 2009; Appelo et al., 2010; Tournassat and Appelo, 2011). The present study aimed to resolve the differences among previous MD simulations. The MD analysis was carried out using montmorillonite as a prototypical smectite clay mineral. This clay mineral contributes predominantly to the specific surface area and surface charge characteristics of most engineered clay barriers and many soil and sedimentary formations. It is also the most widely examined clay mineral in studies of anion diffusion through engineered clay barriers (Bourg and Tournassat, 2015). While the present study focuses on clay interlayers, the results should be relevant to ion and salt migration in other charged nanopores with pore widths below 2 nm.

Equations 1 and 3 obviously have very similar forms. Because of this similarity, the terms “Donnan model,” or “Donnan volume,” are sometimes erroneously used to describe the MEP model despite its molecular-scale assumptions (Appelo and Wersin, 2007; Schoch et al., 2008; Birgersson and Karnland, 2009; Tertre et al., 2011; Tournassat et al., 2011; Glaus et al., 2013). The DE model, however, is based on chemical thermodynamics and, therefore, is fundamentally different from the MGC and MEP models, or any other model involving molecular-scale concepts.

Whether a Donnan equilibrium model or a model based on the Poisson-Boltzmann equation is best suited for predicting anion concentrations in clay mineral interlayers cannot be readily answered by comparing model predictions with macroscopic experimental results. Macroscopic data can be fitted to the DE model by adjusting the clay-phase activity coefficients in θ. Since neither the individual activity coefficient values nor the Donnan potential is accessible experimentally, θ cannot be constrained independently. Similarly, macroscopic anion-exclusion data can be fitted with the MGC or MEP model by adjusting the pore-size distribution or the distance of closest approach of ions to the clay surface (Muurinen et al., 2007; Tournassat and Appelo, 2011; Tournassat et al., 2016). Independent information must, therefore, be sought that can be used to test these modeling approaches. In the following, MD simulations are shown to be insightful in this respect.

Simulations of clay-water systems with three interlayer nanopores in contact with a larger mesopore were carried out with the program LAMMPS (Plimpton, 1995). Water molecules were assumed to be rigid according to the SHAKE algorithm (Ryckaert et al., 1977). The clay mineral structure was also assumed to be rigid, as in previous simulations of this type (Rotenberg et al., 2007b; Hedström and Karnland, 2012; Hsiao and Hedström, 2015), but clay H atoms were allowed to move. Recent MD simulations have shown that the equilibrium distributions of interlayer water and ions (but not their dynamics) are insensitive to the motions of clay atoms (Holmboe and Bourg, 2014). Interatomic interactions were described with the SPC/E model of liquid water (Berendsen et al., 1987), the CLAYFF model of mineral-water interactions (Cygan et al., 2004), and the ion-water interaction parameters of Smith and Dang (1994). Van der Waals interactions between unlike atomic species were derived using the Lorentz-Berthelot combining rules. The CLAYFF model includes Van der Waals and Coulomb interaction parameters for atoms in the bulk clay mineral structure and for certain functional groups on the clay edge surfaces (>Si–OH, >Al₂–OH), but not for other clay mineral edge groups (>AlOH, >Al–O–Si<). For the present study, CLAYFF-like interatomic parameters were derived for edge >AlOH and >Al–O–Si< groups as described in the Supplemental Materials section.

The cis-vacant structure of montmorillonite was taken from Tsipursky and Drits (1984). Individual clay mineral layers (hereafter referred to as TOT layers) were truncated at the (110) and ( $11¯0$⁠) structural planes in order to create edge surfaces, which were then healed by adding −OH or −H groups to under-coordinated edge atoms in order to obtain edge surfaces with zero net proton surface charge (Churakov, 2006). The resulting edge surface structures, which correspond to the AC chains of White and Zelazny (1988), are thought to be the most stable or one of the most stable smectite edge structures (White and Zelazny, 1988; Bickmore et al., 2003; Churakov, 2006; Liu et al., 2014; Newton and Sposito, 2015; Newton et al., 2016). Each TOT layer had a thickness of ~10 Å (ambiguities associated with the definition of the TOT layer thickness are discussed below), a finite length of 45.7 Å in the direction normal to the edge surfaces, and a width of 63.4 Å (effectively an infinite width, because of the periodic boundary conditions of the simulation cell). Each simulation cell contained three TOT layers, three interlayer nanopores, and a larger mesopore (Figure 1). The cell had an initial length of 90 Å in the direction normal to the clay mineral edges (before pressure equilibration), a width of 63.4 Å, and a height of 46.8 Å, 56.7 Å, or 76.5 Å in the direction normal to the clay basal surfaces for simulations of the two-, three-, and five-layer hydrates (referred to hereafter as 2WL, 3WL, and 5WL systems, respectively). Isomorphic substitutions of Mg for Al were randomly distributed on cis octahedral sites. In order to minimize uncertainties associated with the parameterization of edge surface sites, isomorphic substitutions were not allowed to occur near the edge surfaces. The mean number of isomorphic substitutions (108 substitutions randomly distributed between the three TOT layers) was selected in order to obtain a mean basal surface charge density near −0.1 C m⁻², a value typical of montmorillonite (Sposito, 1984).

Each simulation cell included 148 Na⁺ ions and 40 Cl⁻ ions as required to balance the negative surface charge of the TOT layers and establish a NaCl concentration in the mesopore of ~0.5 mol dm⁻³, ~0.4 mol dm⁻³, or ~0.3 mol dm⁻³ in the 2WL, 3WL, and 5WL systems, respectively. The NaCl concentration in the mesopore was chosen such that the MEP model calculations of Hedström and Karnland (2012) (based on Equation 2, with h_EDL equal to one-half the width of a clay interlayer nanopore and with a closest approach distance of a = 0) would predict one Cl⁻ ion in each interlayer nanopore, on average, during the entire simulation. Previous calculations have indicated that, at the ionic strengths examined in the present study, the region where the electrostatic potential spills over between the nanopore and mesopore extends ~1 nm from the clay mineral edges into the interlayer nanopores and into the mesopore (Kraepiel et al., 1999; Jardat et al., 2009). These simulated systems should, therefore, be sufficiently large to contain interlayer and mesopore regions that are not influenced by proximity to the clay edge surfaces.

As in the simulations of Hedström and Karnland (2012), a single Cl ion (hereafter denoted Cl_int) was initially placed in the center of each interlayer nanopore. The three Cl_int ions were held fixed during equilibration, then allowed to move freely during the simulation. Production runs were carried out in the NVT ensemble (constant composition, volume, and temperature) at T = 350 K with a 1 fs time step. Simulation runs were preceded by 0.5 ns of equilibration in the NVT ensemble, 0.5 ns of equilibration in an NP_yT ensemble (constant composition, stress in the direction normal to the edge surfaces, and temperature) at P_y = 1 bar and T = 350 K, and another 2 ns of equilibration in the NVT ensemble. Electrostatic and dispersion interactions beyond 15 Å were computed by the Ewald sum method during part of the simulation (22 ns, 17 ns, and 15 ns for the 2WL, 3WL, and 5WL systems, respectively), then by the similarly accurate particle-particle particle-mesh (PPPM) method (Eastwood et al., 1980; Isele-Holder et al., 2012) during the remaining part of a simulation (up to 53 ns, 40 ns, and 38 ns for the 2WL, 3WL, and 5WL systems, respectively).

The contributions of various electrostatic interactions to the potential energy of all Cl⁻ and Na⁺ ions were calculated using LAMMPS. The contributions of electrostatic interactions with water molecules, clay mineral atoms, and other ions (Na⁺, Cl⁻) to the total per-atom potential energy of Na⁺ and Cl⁻ ions were calculated by canceling the charges of all water atoms, then of all water and montmorillonite atoms during different reexaminations of a trajectory. For the last calculation, where only Na⁺ and Cl⁻ mutual electrostatic interactions were calculated, the system was not neutral (an excess of Na⁺ was present due to the TOT layer surface charge) and the computation of a relevant total electrostatic energy (including pairwise and k-space energy contributions) was not possible. Consequently, only pairwise energies were computed with a cutoff at 25 Å.

All Cl_int atoms left the 2WL and 3WL interlayer space within three nanoseconds of simulation. During the same period of time, Cl⁻ ions from the mesopore entered the 3WL nanopores but did not enter the 2WL nanopores. In the 5WL system, not all of Cl_int atoms left the interlayer space within the same period of time. According to the concentration maps that were computed after 3 ns of simulation, Cl⁻ has access to 3WL and 5WL interlayers, but not to 2WL interlayers (Figure 2). Note that in the 2WL system, the reported Cl⁻ concentration in the interlayer space is truly zero, not simply a low value that is unresolved by the color concentration scale.

Total Cl⁻ and Na⁺ concentration profiles and water density profiles (Figure 3) were computed with a spatial resolution of 0.5 Å on the z-axis in the y [65–85 Å] interval, i.e. in the middle of the interlayer nanopores. The calculated water-density profiles correspond to the concentration of O atoms in water molecules multiplied by the mass of one mole of water. The Na⁺ ions were concentrated in a single peak at the center of the 2WL interlayer space, while two main peaks were present in the 3WL and 5WL interlayer spaces (Figure 3) in agreement with previous MD and Monte Carlo (MC) simulations of Na-montmorillonite with structural charge arising from isomorphic substitutions in the octahedral sheet (Marry et al., 2002; Tambach et al., 2004; Holmboe and Bourg, 2014). The Na density maxima were located at a distance of ~4.5 Å from the siloxane surface O atoms of the TOT layer and separated from the TOT layer by a hydration layer characterized by a water O density peak. Peak positions were consistent with previous studies of clay mineral interlayer nanopores and larger pores (Marry et al., 2002; Tournassat et al., 2009; Bourg and Sposito, 2011). The height of the Na density peaks increased with the number of isomorphic substitutions in the nearby TOT layer, as expected from charge balance considerations. The height of the main water density peak was correlated, but was not strictly proportional, to that of the Na peak, indicating that the structure of the interfacial water layer was significantly influenced by both the clay surface and the adsorbed cations. Few Na⁺ ions were engaged in inner sphere surface complexes. In the case of Cl, density profiles showed significant anion exclusion even at the mid-plane of the nanopores (Figure 3). The distance of closest approach of Cl⁻ to the plane of O atoms on the basal surface was ~ 3.5 Å. This value coincides with the distance of closest approach of Na⁺ ions that were not bound in innersphere surface complexes.

At the scale of an individual interlayer nanopore, average concentrations calculated from the MD concentration profiles depend on an arbitrary choice of the boundary between the pore space and the solid phase. The DE model calculations presented in a previous section are based on a molality convention (mol per kilogram of water) while the MGC and MEP model calculations are based on a molarity convention (mol dm⁻³) with the underlying assumption of a water density of 1 kg dm⁻³ (note that the effect of temperature on water density was neglected). Consequently, the location of the clay mineral-water boundary was identified based on the choice that yields an average interlayer water density of 1 kg dm⁻³. The resulting effective TOT layer thickness equals 9.78 ± 0.45 Å for all three systems investigated in the present study (Table 1). Accordingly, average interlayer ion densities were calculated based on interlayer nanopore widths of 5.7, 9.1, and 15.6 Å in the 2, 3, and 5WL systems, respectively (Table 1).

For comparison with the mean interlayer concentrations predicted by MD simulation, interlayer concentrations of Na⁺ and Cl⁻ were calculated using the MGC and MEP models. The MGC model predictions were obtained by solving the PB equation (Equation 2) numerically using the values of σ and $cisoln$ from the MD simulation as input parameters (Table 2). The distance of closest approach of ions to the solid surface was identified based on the MD concentration profiles, using either the location of the ‘foot’ of the Na⁺ density profiles or the location of the Na⁺ density peak (Figure 3 and Table 1). In the case of the 2WL system, the second option was not available because only one Na⁺ density peak was present. The MEP model calculations assumed that the EDL filled the entire interlayer pore space and involved the same assumptions as the MGC model regarding the distance of closest approach of ions to the clay mineral surface (Table 2). Predicted Na⁺ and Cl⁻ density profiles were converted to average concentrations per mass of interlayer water in order to ensure consistency between different results.

The per-atom electrostatic energy corresponds to the coulomb energy attributed to each individual Cl⁻ (Na⁺) ion in the system. For the 2WL system (Figure 4), the computation of the per-atom energy for Cl⁻ in the interlayer space was not possible because chloride did not enter the interlayer space during the simulation. On the other hand, the total per-atom electrostatic energy was the same in the mesopore as in the interlayer in the 3WL and 5WL systems (Figures 5 and 6). The most remarkable outcome of this calculation is that, in the 3WL and 5WL systems, the electrostatic energy arising from Cl⁻ solvation was much higher (less favorable interaction) in interlayer nanopores than in the mesopore. The opposite effect was not observed in the case of Na⁺; electrostatic interaction of Na⁺ ions with water molecules was approximately equally favorable in the interlayer nanopores as in the mesopores. A similar asymmetry in Na⁺ and Cl⁻ solvation energies has been observed in model polar nanopores and ascribed to hindered orientations of confined water molecules hydrating Cl⁻ in nanopores (Crozier et al., 2001). As expected, interactions with other ions in the interlayers were attractive in the case of Cl⁻ and repulsive in the case of Na⁺ (lower row in Figures 4, 5, and 6).

In the case of the 2WL hydrate, no Cl⁻ ion entered the interlayer space during the course of the simulation, in agreement with the modeling results of Rotenberg et al. (2007b), but in disagreement with those of Hsiao and Hedström (2015). The Cl⁻ concentrations expected from the energy barriers predicted by Rotenberg et al. (2007b) or Hsiao and Hedström (2015) for bulk NaCl concentrations similar to those obtained in the present study are given in Table 3, together with the average number of Cl⁻ ions expected in the interlayer space during the course of an entire simulation. The simulations were carried out at the same temperature (350 K) as the simulations of Hsiao and Hedström (2015) and with similar simulation times (50 ns vs. 100–200 ns) and volumes (27×10⁴ ųvs. 15×10⁴ ų), thus ensuring roughly equally reliable output statistics. The fact that Cl⁻ ions did not enter the interlayer space cannot, therefore, be attributed to a lack of convergence in the present simulation, as Hsiao and Hedström have postulated to explain the difference between their results and those of Rotenberg et al. (2007b).

The mapping of the per-atom electrostatic energy contributions evinced the major role played by solvation in inhibiting the entry of Cl⁻ into the interlayer space. This effect was particularly pronounced in the 2WL system, perhaps because a significant fraction of the water molecules associated with the water density peaks nearest to the clay surfaces (i.e. most of the interlayer water in the 2WL hydrate) were engaged in solvating Na⁺ outer-sphere surface complexes while also donating a hydrogen bond to the siloxane surface O atoms (Marry et al., 2008; Bourg and Sposito, 2011). In short, water located in the first monolayer on the clay surfaces had a very different structure from bulk liquid water (Low, 1979), and this difference became acute in the case of the 2WL hydrate. These findings are consistent with the hypothesis formulated by Hsiao and Hedstrom (2015) that differences in solvation energy play an important role in inhibiting the entry of Cl⁻ in the interlayer space.

Molecular dynamics simulations are well known to be sensitive to the choice of parameters used to describe interatomic interactions. Consequently, differences in the extent of anion exclusion predicted by Rotenberg et al. (2007b) and Hsiao and Hedström (2015) may reflect differences in the parameterization of interatomic interactions in these two studies. The present work used the same interatomic potential models as Hsiao and Hedstrom (2015), namely, the CLAYFF and SPC/E models. To the best of the present authors’ knowledge, the only major methodological difference between the present simulations and those of Hsiao and Hedstrom (2015) is that these latter authors treated the negative structural charge of the clay mineral particles as uniformly distributed over all octahedrally coordinated atoms in the clay mineral structure. In the present simulations, the negative charge arising from each Al → Mg isomorphic substitution was assigned to the site of substitution and its six surrounding O atoms in accordance with the original CLAYFF model. The treatment of clay mineral negative structural charge used in the present study was derived from quantum mechanical calculations of the electrostatic potential near clay mineral surfaces (Cygan et al., 2004), which predicts interlayer water and ion distributions and dynamics in excellent agreement with experimental results (Bourg and Sposito, 2010; Ferrage et al., 2011; Marry et al., 2011). The present results do not allow an accurate quantification of the height of the energy barrier opposing Cl⁻ entry into the interlayer space of the 2WL hydrate, but strongly suggest that this barrier is greater than that predicted by Hsiao and Hedstrom (2015).

The numerical solution of the PB equation for the 3WL and 5WL hydrates was consistent with the average interlayer Na⁺ concentration and with the interlayer Na⁺ concentration profile predicted by the MD simulations (Figure 3). This close correspondence was particularly true if the distance of closest approach of ions to the clay mineral surface was identified with the foot of the Cl⁻ density profiles, which yields a = 2.0±0.8 Å. This predicted a value was in good agreement with the range of values regarded as a typical ion size in studies of ion adsorption (a =1.8 to 2 Å) (Sposito, 1992; Kim and Netz, 2006) and with the hydrated radii of Na⁺ and Cl⁻ (1.7 and 2.0 Å, respectively) (Jardat et al., 2009). The agreement between PB calculations and MD simulation predictions was somewhat worse in the case of the Cl⁻ concentration profiles than in the case of the Na⁺ profiles (Figure 3), perhaps reflecting the poorer statistics for interlayer Cl concentrations or the influence of short-range ion-ion interactions (and possibly ion-water interactions, as noted above) that are not accounted for in the PB equation. Nevertheless, reasonable quantitative agreement was found (Table 2). These calculations further showed that, for the systems examined, the average interlayer Na⁺ and Cl⁻ concentrations predicted with the MGC and MEP models were essentially identical, i.e. a full numerical resolution of the PB equation was not necessary to predict average ion concentrations in the nanopores. Instead, predicted interlayer anion concentrations were primarily sensitive to the distance of closest approach used in the MGC or MEP models, the best fit being obtained with a distance of closest approach intermediate between those obtained using the ‘closest approach’ and ‘Na peaks’ approximations. The similarity of the MGC and MEP models in the case of monovalent ions is consistent with recent predictions of the density profiles of Na⁺, Cl⁻, and Br⁻ (but not Ca²⁺) in a 31.5 Å-wide clay interlayer nanopore (Tinnacher et al., 2016).

According to Hsiao and Hedstrom (2015), solving the Poisson-Boltzmann equation in order to model the behavior of ions in clayey media is not necessary, because the use of the Donnan equilibrium (DE) model is warranted for compacted clay systems. A strict application of Henry’s Law (and, hence, of the DE model) to clay systems would require that the counterions be dissociated fully from the mineral surface, such that mean ionic concentrations in the hydrated clay phase be unambiguously determined (Babcock, 1963). The present MD simulations were carried out at a relatively high salinity and, therefore, the simulations cannot determine whether Henry’s Law is valid for counter-ions in a clay system. However, MGC model calculations (Sposito, 1992), X-ray reflectivity experiments (Schlegel et al., 2006; Lee et al., 2010, 2012), and MD simulation studies (Marry et al., 2002; Rotenberg et al., 2007a; Holmboe and Bourg, 2014) have provided strong evidence that cations remain condensed at clay mineral–water and mica–water interfaces, even under bulk solution conditions approaching infinite dilution. The results presented here indicate MGC model calculations and MD simulations are mutually consistent, even in highly compacted clays, and, therefore, suggest that the DE model cannot be valid, even in highly compacted clay.

The PB model does not rely on the validity of Henry’s Law. Consequently, the MEP model predictions are similar to those of the DE model and it remains valid within the limits of the underlying molecular-scale hypotheses and approximations. The interlayer pore scale of the 3WL system can be considered as the lower limit for water content in this model. Below this scale, the mean force potential that applies to Cl⁻ ions entering the interlayer region cannot be approximated using mean field theories, because Cl⁻ is significantly more weakly solvated by confined water than by bulk liquid water. Within the context of the MEP (or MGC) model, this effect can be modeled in a straightforward manner by accounting for the distance of closest approach of ions to the clay surface, a ~ 2 Å.

The present MD simulations demonstrated that anion exclusion from clay interlayer nanopores is significantly greater than expected from mean field theories, such as the modified Gouy-Chapman (MGC) or the simplified mean electrostatic potential (MEP) models. The calculations showed that this effect can be viewed as arising from the weaker solvation of Cl⁻ by confined water than by bulk liquid water and it can be modeled by accounting for the closest approach distance of ions to the clay surface, a ~ 2 Å. Calculations reported here further suggest the existence of a threshold value in nanopore width below which Cl⁻ cannot readily access clay interlayer nanopores.

Comparison of the MD simulation results with the underlying hypothesis of the Donnan model established that the Donnan model cannot be applied to clay-water mixtures because these systems do not satisfy the two basic assumptions made about the behavior of ions in a Donnan system. This problem is related to the geometry of the charged surfaces (infinite plane), which leads to counterion condensation in the vicinity of the surface (Sposito, 1992), thereby preventing the full dissociation of Na⁺ ions from the surface at infinite dilution required by the Donnan model. Although a mean electrostatic potential model derived from the Poisson-Boltzmann equation may be equivalent numerically to the Donnan model, it is fundamentally different in terms of the underlying hypotheses.

The ability of Cl⁻ to enter clay interlayer nanopores decreases sharply as pore size decreases from 9.1 Å (3WL hydrate) to 5.7 Å (2WL hydrate). This finding should be taken into account when predicting the membrane properties of compacted clay materials. The present results suggest that clay materials may behave as nearly ideal semi-permeable membranes under conditions where only 2WL interlayer spaces (or inter-particle pores of similar width) connect the diffusion paths between larger pores, in agreement with recent experimental results by Chagneau et al. (2015).

This work was supported by the Director, Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, of the U.S. Department of Energy under Contract DE-AC02-05CH11231. C. Tournassat acknowledges funding from L’Institut Carnot BRGM for his visit to the Lawrence Berkeley National Laboratory, and from the CNRS-NEEDS project TRANSREAC. G. Sposito acknowledges funding from the University of California at Berkeley under the auspices of a Chancellor’s Professorship. The research reported in this paper used resources of the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract DE-AC-02-05CH11231.