NUREG/CR-6757

Large-Scale Molecular 

Dynamics Simulations of 

Metal Sorption onto the 

Basal Surfaces of Clay 

Minerals 

Sandia National Laboratories 

,AG••R REG 

U.S. Nuclear Regulatory Commission 

Office of Nuclear Regulatory Research 

Washington, DC 20555-0001
AVAILABILITY OF REFERENCE MATERIALS 

IN NRC PUBLICATIONS

nnC neTerence Material

As of November 1999, you may electronically access NUREG-series publications and other NRC records at 

NRC's Public Electronic Reading Room at www. nrc.gov/NRC/ADAMS/index.html. 

Publicly released records include, to name a few, NUREG-series publications; Federal Register notices; applicant, licensee, and vendor documents and 

correspondence; NRC correspondence and internal 

memoranda; bulletins and information notices; inspection and investigative reports; licensee event 

reports; and Commission papers and their 

attachments. 

NRC publications in the NUREG series, NRC regulations, and Title 10, Energy, in the Code of 

Federal Regulations may also be purchased from one of these two sources. 

1. The Superintendent of Documents 

U.S. Government Printing Office 

Mail Stop SSOP 

Washington, DC 20402-0001 

Internet: bookstore.gpo.gov 

Telephone: 202-512-1800 

Fax: 202-512-2250 

2. The National Technical Information Service 

Springfield, VA 22161-0002 

www.ntis.gov 

1-800-553-6847 or, locally, 703-605-6000 

A single copy of each NRC draft report for comment is available free, to the extent of supply, upon written 

request as follows: 

Address: Office of the Chief Information Officer, 

Reproduction and Distribution 

Services Section 

U.S. Nuclear Regulatory Commission 

Washington, DC 20555-0001 

E-mail: DISTRIBUTION@nrc.gov 

Facsimile: 301-415-2289 

Some publications in the NUREG series that are posted at NRC's Web site address 

www.nrc.gov/NRC/NUREGS/indexnum.html 

are updated periodically and may differ from the last 

printed version. Although references to material found on a Web site bear the date the material was accessed, the material available on the date cited may subsequently be removed from the site.

Non-NRC Reference Material

Documents available from public and special technical 

libraries include all open literature items, such as books, journal articles, and transactions, Federal 

Register notices, Federal and State legislation, and 

congressional reports. Such documents as theses, dissertations, foreign reports and translations, and 

non-NRC conference proceedings may be purchased from their sponsoring organization. 

Copies of industry codes and standards used in a substantive manner in the NRC regulatory process are maintained at

The NRC Technical Library 

Two White Flint North 

11545 Rockville Pike 

Rockville, MD 20852-2738 

These standards are available in the library for 

reference use by the public. Codes and standards are usually copyrighted and may be purchased from the 

originating organization or, if they are American 

National Standards, from

American National Standards Institute 

11 West 4 2nd Street 

New York, NY 10036-8002 

www.ansi.org 

212-642-4900

Legally binding regulatory requirements are stated 

only in laws; NRC regulations; licenses, including 

technical specifications; or orders, not in 

NUREG-series publications. The views expressed in contractor-prepared publications in this series are not necessarily those of the NRC. 

The NUREG series comprises (1) technical and 

administrative reports and books prepared by the 

staff (NUREG-XXXX) or agency contractors 

(NUREG/CR-XXXX), (2) proceedings of 

conferences (NUREG/CP-XXXX), (3) reports resulting from international agreements (NUREG/IA-XXXX), (4) brochures 

(NUREG/BR-XXXX), and (5) compilations of legal decisions and orders of the Commission and 

Atomic and Safety Licensing Boards and of 

Directors' decisions under Section 2.206 of NRC's 

regulations (NUREG-0750).

DISCLAIMER: 

Neither 

This report was prepared as an account of work sponsored by an agency of the U.S. Government. the U.S. Government nor any agency thereof, nor any employee, makes any warranty, expressed or implied, or assumes any legal liability or responsibility for any third party's use, or the results of such use, of any information, apparatus, product, or process disclosed in this publication, or represents that its use by such third party would not infringe privately owned rights.
NUREG/CR-6757

Large-Scale Molecular Dynamics 

Simulations of Metal Sorption onto the 

Basal Surfaces of Clay Minerals 

Manuscript Completed: October 2001 

Date Published: February 2002 

Prepared by 

D. M. Teter, R. T. Cygan 

Sandia National Laboratories 

Albuquerque, NM 87185-0750 

E. O'Donnell, NRC Project Manager 

Prepared for 

Division of Systems Analysis and 

Regulatory Effectiveness 

Office of Nuclear Regulatory Research 

U.S. Nuclear Regulatory Commission 

Washington, DC 20555-0001 

NRC Job Code W6811
Abstract 

Large-scale molecular dynamics computer simulations have been completed to evaluate the sorption of Cs', Sr 2

+, and 

other related ions onto the (001) basal planes of kaolinite and montmorillonite. Simulations were undertaken as a 

function of metal ion concentration to evaluate sorption mechanisms and to determine distribution coefficients 

(Kd's). The simulations show that Cs' ions bind to both kaolinite and montmorillonite, while Sr2

' ions tend to re

main in solution. At high ionic strengths (greater than 0.5 M), the fundamental mechanism for Cs' sorption to the 

aluminol plane of kaolinite is through the formation of CsCI complexes where the CI ion is loosely bound to the 

outer hydroxyl groups. At lower concentrations (0.1 M and 0.01 M), the Cs' ions do not sorb to the aluminol plane, 

but prefer to bind adjacent to the center of the six-membered ring silica tetrahedral sheet as an inner-sphere complex. 

We have determined a Kd of 10 ± 4 ml/g for kaolinite to sorb CsW ions from a 0.01M CsCl solution, and a Kd of 0 for 

kaolinite to sorb Sr 2

+ ions from a 0.01M SrC12 solution. For the case of montmorillonite, we have determined a Kd of 

293 ± 30 ml/g to sorb Cs' ions from a 0.1M CsCl solution, and a Kd of 11± 4 ml/g to sorb Sr 2

+ ions from a 0.1M 

SrCI2 solution. The montmorillonite sorbs the cations more effectively than kaolinite because it has negatively

charged sorption sites that are created by isomorphic substitutions on both tetrahedral and octahedral sites. Calcula

tions show that the Sr2

÷ ion is solvated more strongly than Cs' and is not sorbed by the kaolinite or montmorillonite 

surfaces. Directions for further work are outlined and prioritized.

iii
Contents 

Abstract iii Forew ard ................... ........ ................................................... vii 

Acknowledgm ents ....................................................................... viii 

1. INTRODUCTION 1-1 

2. MATERIAL AND FINITE ELEMENT ANALYSIS MODELS .................................. 2-1 

2.1 Kaolinite Structure .............................................................. 2-1 

2.2 M ontmorillonite Structure ......................................................... 2-1 

3. MOLECULAR MODELING METHODS .................................................. 3-1 

3.1 Validation of Interatomic Potentials ..................................................... 3-2 

4. MOLECULAR MODELING RESULTS ................................................... 4-1 

4.1 M odeling the Kaolinite Unit Cell ................................................... 4-1 

4.2 The Solvation of Cs+ in W ater ..................................................... 4-1 

4.3 M odeling of the Cs+Cl- Complex ................................................... 4-2 

4.4 Modeling of Cesium and Strontium Sorption onto the Kaolinite (001) Surface ................. 4-3 

4.4.1 Preliminary Simulations ..................................................... 4-3 

4.4.2 Production Simulations ...................................................... 4-6 

5. CALCULATION OF DISTRIBUTION COEFFICIENTS .................................. 5-1 

5.1 K aolinite ...................................................................... 5-1 

5.2 M ontm orillonite ................................................................ 5-1 

5.3 Energetics and Kd Values ......................................................... 5-1 

6. CON CLU SION S .................................................................... 6-1 

7. SUGGESTIONS FORFURTHER WORK .................................................. 7-1 

8. REFERENCES ............................................................. 8-1 

v
Figures

Figure 1. A graphical representation of the CI unit cell of kaolinite. The top sheet is composed of A10 6 

octahedra (0) and the bottom sheet is composed of SiO4 tetrahedra (T) .................................................... 2-1 

Figure 2. A graphical representation of the unit cell of montmorillonite without interlayer cations and waters ....... 2-2 

Figure 3. A graphical representation of a Cs' ion hydrated by six waters in an octahedral geometry ....................... 4-2 

Figure 4. Representative configuration for the preliminary analysis of the equilibrated Cs' sorption onto the 

(00 1) basal surfaces of kaolinite .................................................................................................................. 4-4 

Figure 5. Density profile of the equilibrated 1 M CsCl solution interacting with the two basal surfaces of 

kaolinite based on the accumulation of 100 psec of trajectory analysis from the preliminary molecular 

dynam ics sim ulation .................................................................................................................................... 4-4 

Figure 6. Density maps of the equilibrated 1 M CsC1 solution interacting with the aluminol basal surface of 

kaolinite, based on the accumulation of 100 psec of trajectory analysis from the preliminary molecular 

dynam ics sim ulation .................................................................................................................................... 4-5 

Figure 7. Density maps of the equilibrated 1 M CsCl solution interacting with the siloxane basal surface of 

kaolinite based on the accumulation of 100 psec of trajectory analysis from the preliminary molecular 

d ynam ics sim u lation .................................................................................................................................... 4-6 

Figure 8. A graphical representation of a 0.5 M CsCl solution in steady-state equilibrium with a (001) kaolinite 

surface ......................................................................................................................................................... 4 -7 

Figure 9. A graphical representation of a 0.1 M CsC1 solution in steady-state equilibrium with a (001) kaolinite 

surface ......................................................................................................................................................... 4 -7 

Figure 10. A graphical representation of a 0.01 M CsCI solution in steady-state equilibrium with a (001) 

kao linite surface ........................................................................................................................................... 4 -8 

Figure 11. Details of the 0.01 M CsCI solution in steady-state equilibrium with a (001) kaolinite surface .............. 4-8 

Figure 12. A graphical representation of a 0.1 M Cs' solution in steady-state equilibrium with a (001) 

m ontm orillonite surface ............................................................................................................................... 4-9 

Figure 13. A graphical representation of a 0.1 M Sr 2

÷ solution in steady-state equilibrium with a (001) 

m ontm orillonite surface ............................................................................................................................... 4-9 

Figure 14. A graphical representation of 8 Cs÷ and 8 Sr 2

÷ ions in steady-state equilibrium with a (001) 

m ontm orillonite surface ............................................................................................................................... 4-9 

Tables 

Table 1. Forcefield parameters used for the Coulombic and Lennard-Jones terms of the nonbonded potential 

energy expression represented by Equation (1) ........................................................................................... 3-1 

Table 2. A comparison of the kaolinite unit cell parameters as determined by experiment, DFT-GGA 

calculations, and the forcefield parameters used in this study .................................................................... 4-1 

Table 3. Equilibrium Cs+-Owate distances as the water coordination number increases calculated using density 

functional theory and the forcefield described in Table 1 ............................................................................ 4-2 

Table 4. Calculation of the cohesive energy of the Cs+CI complex as a function of ionic separation using DFT 

and the forcefield parameters described in Table I ..................................................................................... 4-2 

Table 5. Comparison of energy for various ions based on water solvation or clay surface sorption environment.... 5-2

vi
Foreword

This contractor technical report was prepared by Sandia National Laboratories under their Interagency Work Order 

(JCN W681 1) with the U.S. Nuclear Regulatory Commission. 

This report summarizes the results of a task to apply molecular dynamics computer simulations for the sorption of monovalent and divalent metallic cations onto clay minerals. In particular, the report focuses on the sorption of Cs' and Sr". The results of these simulations give insights into the basic mechanisms controlling metal sorption onto 

common soil components such as clay minerals. The information on mechanisms controlling sorption reported here 

are being applied by Sandia to enhance the prediction of radionuclide transport through soil. 

This report is not a substitute for NRC regulations, and compliance is not required. The approaches and/or methods described in this NUREG/CR report are provided for information only. Publication of this report does not necessarily constitute NRC approval or agreement with the information contained herein. Use of product or trade names is for identification purposes only and does not constitute endorsement by the NRC or Sandia National Laboratories. 

Cheryl A. Trottier, Chief 

Radiation Protection, Environmental Risk, and Waste Management Branch 

Division of Systems Analysis and Regulatory Effectiveness 

Office of Nuclear Regulatory Research

vii
Acknowledgments

Special thanks are extended to Andrey Kalinichev for his assistance in performing the atom density maps in the 

analysis of the preliminary kaolinite simulations. The reviews provided by Louise Criscenti and Henry Westrich of 

an early draft of the manuscript are gratefully acknowledged. This work was supported by the U.S. Nuclear Regula

tory Commission, Office of Nuclear Regulatory Research. We are very grateful for the advice and support provided 

by our NRC program manager, Dr. Edward O'Donnell, during the funding period. Sandia National Laboratories is a 

multiprogram laboratory operated by the Sandia Corporation, a Lockheed Martin Company, for the United States 

Department of Energy under Contract DE-AC04-94AL85000.

Viii
1. INTRODUCTION

Radioactive 137Cs and 9°Sr have been dispersed globally 

by nuclear testing, the production of weapons grade 

uranium and plutonium, and more recently by the nu

clear accident at Chernobyl. The fate of 1

37

Cs and other 

radionuclides in the environment is controlled primarily 

by the ability of subsurface and colloidal minerals to 

naturally attenuate these contaminants through irre

versible sorption and precipitation processes. 137CS is 

also a concern because of the phenomena known as "cask weeping", whereupon initially sorbed 1

37

Cs de

sorbs from nuclear fuel transport casks during shipping 

due to changes in environmental chemistry (Teter et al., work in progress at SNL). Our ability to evaluate these 

processes is based on experimental methods such as 

surface charge titrations, sorption isotherm measure

ments, atomic force microscopy (AFM), nuclear mag

netic resonance (NMR) and X-ray adsorption spectro

scopies (EXAFS and XANES). However, despite the 

power of the available analytical methods, many of the 

details of the sorption process elude us because of the 

structural and chemical complexities of geomedia such 

as clays. The application of atomistic methods allows 

us to construct and test detailed models of chemical and 

physical processes such as sorption. When combined 

with experiments, these methods can provide a funda

mental understanding of sorption mechanisms that enables us to design more efficient methods for waste 

treatment and improve the immobilization and seques

tration of contaminants (Krumhansl et al., 2001). 

We are studying the sorption of Cs' and Sr 2

, on kaolin

ite surfaces because kaolinite provides a relatively sim

ple substrate for studying surface sorption due to its 

lack of interlayer metal sites. The sorption of Cs' on 

kaolinite has been experimentally verified by Kim et al. 

(1996a, 1996b, 1997, 1998), Komarneni (1978), Krum

hansl et al. (2001), Michaelian (1998a, 1998b), Saiers 

and Hornberger (1999), Weiss et al. (1990), and We

strich et al. (1996). The sorption isotherms suggest a 

pH-dependent affinity of Cs for the Al octahedral edge 

sites, with possibly weak sorption on the Si tetrahedral 

edge sites at high-pH values. Binding of Cs' to the 

basal Al octahedral sites appears to be independent of 

pH. These results are supported by titration experi

ments to determine the effect of pH on the surface 

charge of kaolinite and the subsequent surface com

plexation modeling (Brady et al. 1996, 1998). Nagy et 

al. (1995) obtained AFM images of individual kaolinite 

grains that provide evidence for significant edge surface 

area, possibly up to 50% of the surface area. These 

AFM images provide additional support for the strong 

influence of edge sites on Cs sorption.

1-1
2. MATERIAL AND FINITE ELEMENT ANALYSIS MODELS

2.1 Kaolinite Structure 

Kaolinite (Al2Si2O5(OH)4) is one of the simplest layered 

clays. Kaolinite is a layered aluminosilicate with a 1:1 

uncharged layered structure. Each layer consists of a 

sheet of A10 6 octahedra (0) building four-membered 

aluminate rings that are connected via common oxygen 

atoms to a sheet of SiO 4 tetrahedra (T) forming six

membered silicate rings. Only two of the three possible aluminum sites are occupied in the octahedral layer. The kaolinite structure contains two types of hydroxyl 

groups, both of which are contained in the Al octahedral 

sheet. The interlayer hydroxyl groups are bound to the 

hexa-coordinated aluminum atoms and cover one of the 

two basal surfaces of each layer. These hydroxyl 

groups form hydrogen bonds with the adjacent Si tetra

hedral sheet in the next TO layer repeated in the (001) 

direction. The inner hydroxyl groups are situated 

within the shared plane of the octahedral/tetrahedral 

sheets and are associated with the vacant octahedral Al 

sites. 

To examine the specific molecular controls on Cs sorp tion, we used an atomistic model of a kaolinite surface 

derived from the unit cell shown in Figure 1. This 

model is based upon the unit cell parameters recently 

determined by Neder et aL (1999) in a single-crystal 

0 

T

structure refinement using synchrotron radiation, with 

the internal parameters determined by Bish (1989, 1993) from Rietveld refinement of the neutron powder 

data. These unit cell parameters agree well with the 

previous studies by Young and Hewatt (1988) and Ad

ams (1983). 

2.2 Montmorillonite Structure 

Montmorillonite (approximately Na 3(Si 31Al)(Al1 4Mg2 ) 

Os0(OH) 16.nH 20) where this formula represents four 

formula units and is comprised of an octahedral sheet, is 

sandwiched between two tetrahedral sheets and is clas

sified as a 1:2 layer silicate mineral. In contrast to kao

linite, montmorillonite has a layered structure that pos

sesses a net negative charge due to isomorphic 

substitutions of metal ions in either the tetrahedral, oc

tahedral, or both sheets. Typically, aluminum substi

tutes for silicon atoms in the tetrahedral sheet and mag nesium for aluminum in the octahedral sheet (see formula). The charge deficiency allows for the intro

duction of interlayer cations, often hydrated by water 

molecules, that balance out the negative layer charge. The relatively small charge imbalance of the layer struc

ture and the strength of the hydrogen bonding associ

ated with the interlayer waters, allow montmorillonite 

clays to swell, like other smectite clays. The

Figure 1. A graphical representation of the C1 unit cell of kaolinite. The top sheet is composed of A10 6 octa hedra (0) and the bottom sheet is composed of SiO4 tetrahedra (T). The interlayer hydroxyls are visible above the 0 sheet and the inner hydroxyls are visible between the 0 and T sheets.

2-1

ý_ kk- (,., v
octahedral layer consists of a sheet of A10 6 that form 

four-membered aluminate rings. These rings are con

nected via common oxygen atoms to two sheets of SiO4 

tetrahedra, forming six-membered silicate rings. Like 

kaolinite, only two of the three possible aluminum sites 

are occupied in montmorillonite; both clay minerals are 

often referred to as dioctahedral clays. The hydroxyl 

groups of montmorillonite are situated within both of 

the shared planes and are associated with the vacant 

octahedral aluminum sites. Refined structures for 

montmorillonite and other smectite clays are lacking in 

0- -

AM

literature, due to the sub-micron grain size of these ma

terials, which prevents any quality single crystal re

finement, and the extensive stacking disorder that oc

curs in the z-direction which precludes the long range 

ordering needed to enhance diffraction signals. Fortu

nately, there have been several conceptual models for 

montmorillonite which provide a sound basis for the 

crystal structure. The montmorillonite model presented 

in Figure 2 is based on the electron diffraction determi

nation by Tsipursky and Drits (1984).

Mg

Al

Figure 2. A graphical representation of the unit cell of montmorillonite without interlayer cations and wa

ters. The top and bottom sheets are composed of SiO 4 tetrahedra (T), while the middle sheet is composed of 

AI0 6 octahedra (0). The darkened tetrahedra in the top right sheet represents an aluminum substitution, 

while the two light octahedra in the middle rear sheet represent magnesium substitutions (note arrows). 

These substitutions create negatively charged sites that are normally charged balanced by interlayer Na ions.

2-2
3. MOLECULAR MODELING METHODS

The computer simulation methods in this study de

termine the potential energy of the system using a set 

of analytical functions referred to as forcefields. This 

model calculates the energy and forces in the system 

based upon the atomic coordinates of the atoms. En

ergy contributions are derived from non-bonded in

teractions, such as the long-range Coulombic terms, 

as well as short-range repulsions and attractive Van 

de Waals interactions. Contributions to the potential 

energy also arise from bonded terms such as bond 

stretching and bond angle bending terms. 

The non-bonded interactions are the critical energy 

contributions that will control how a molecule or ion 

will approach and interact with a mineral surface. 

The analytical expression (Equation (1)) for the non

bonded energy expression includes the Coulombic 

and Lennard-Jones potential terms. The potential 

energy Uij for the interaction between ion i and ion j 

at an ionic separation of rij is given by: 

-2 (1) 

where qj and qj are the ionic charges. The Lennard

Jones parameters, gij and oij represent the potential 

well depth (in kJ/mol) and the interatomic distance 

(in A) at the energy minimum. The first term of 

Equation (1) represents the Coulombic energy and

can be either positive or negative, depending on the 

charges on the interacting ions and is the dominant 

non-bonded energy term. The next two terms repre

sent the short-range repulsive energy term and the 

attractive Van der Waals term as a function of ionic 

separation. The forcefield parameters used in this 

study are given in Table I with the short-range pa

rameters determined using the mixing rules given in 

Equations (2a) and (2b):

-16i. = ;.6. aij = - (2a,b) 2

The bonded interactions within the water molecule 

and hydroxyl groups are treated as harmonic oscilla

tors (Eq. 3):

U. (X) = L (3)

where x is an interatomic bond length or angle, x, is 

its corresponding equilibrium value and k, the force 

constant. For the flexible SPC water model (Telemen 

et al., 1987) used in this study, the bond stretching 

force constant for the water hydroxyl groups is 

keo M=4637 kJ mol"' A-2 and the equilibrium bond 

length is xeOH=1.0 A, while the bond angle force con

stant is 383 kJ'mol"1 rad2 , and the equilibrium H-O-H 

bond angle is 109.470.

Table 1. Forcefield parameters used for the Coulombic and Lennard-Jones terms of the nonbonded potential 

energy expression represented by Equation (1). 

Ion Charge e (kJ/mol) F (A) Reference 

Aloct +1.575 5.567x10 6 4.271 Cygan et aL. (2002) 

Mgoc, +1.360 0.902x10 6 5.265 Cygan et al. (2002) 

Silre +2.100 7.706x 10.6 3.302 Cygan et al. (2002) 

Hhvdroxv] +0.425 0.0000 0.000 Cygan et aL (2002) 

Ohvdroxyl -0.950 0.6502 3.166 Cygan et al. (2002) 

Obond -1.050 0.6502 3.166 Cygan et al. (2002) Cs÷ +1.000 0.0800 4.650 This study 

Na÷ +1.000 0.1000 2.750 This study 

Sr +2.000 0.0900 3.700 This study Cl -1.000 0.1600 4.650 This study 

Hwater +0.410 0.0000 0.000 Telemen et al. (1987) 

Owater -0.820 0.6502 3.166 Telemen et al. (1987)

3-1
3.1 Validation of Interatomic Po

tentials 

A critical stage of any model is validating it against 

the best available results before extending to predict 

new phenomena. In order to validate the forcefield

parameters to be used in this study, we compared 

their ability to reproduce experimental findings, as 

well as results obtained using quantum mechanical 

calculations. The following sections describe the 

ability of these forcefields to reproduce the kaolinite 

unit cell, the structure of Cs÷ water clusters and the 

structure and energetics of the CsCl complex.

3-2
4. MOLECULAR MODELING RESULTS

4.1 Modeling the Kaolinite Unit 

Cell 

The equilibrium unit cell parameters of kaolinite were 

calculated using an energy minimization scheme im

plementing the forcefield parameters. These results 

(Table 2) compare favorably with the recent experi

mental study by Neder et al. (1999) using synchrotron 

radiation. Furthermore, we used density functional 

theory (DFT) (Hohenberg and Kohn, 1964; Kohn and 

Sham, 1965) within the PW91 (Perdew and Wang, 

1992) generalized gradient approximation (GGA) for 

electronic exchange and correlation to calculate the 

zero-pressure kaolinite unit cell. These calculations 

used ultrasoft pseudopotentials (Vanderbilt, 1990) 

with a plane-wave basis set (Pickett, 1989). A kinetic 

energy cut-off of 500 eV and a [222] Monkhorst-Pack 

Brilloin zone (Monkhorst and Pack, 1976) was used 

for the calculation. The kaolinite unit cell parameters 

determined by the forcefield calculations compare 

favorably with those determined using DFT. These 

results suggest that the forcefield parameters used in 

this study are capable of modeling the kaolinite sur

face with a reasonable degree of accuracy. 

4.2 The Solvation of Cs' in Water 

In order to validate the ability of the Cs' forcefield 

(Smith and Dang, 1994b) used in this study, we un

dertook a series of simulations to determine the abil

ity of this forcefield to model Cs' solvation in an 

aqueous solution. 

The solvation of a Cs' ion in water was investigated 

using molecular dynamics simulations. The computer

simulation consisted of a 18 A x 18 A x18 A box 

containing 200 water molecules and a Cs' ion using 

the forcefield parameters given in Table 1. An NVT 

(constant number of ions, constant volume, and con

stant temperature) molecular dynamics simulation 

was performed at a temperature of 300K using a 

timestep of 0.0005 psec to capture the high frequency 

hydroxyl vibrations. At equilibrium (5 psec), the 

radial distribution function indicates a coordination of 

approximately eight water molecules in the first coor

dination sphere at a mean distance of approximately 

3.1 A from the Cs' ion. X-ray and neutron diffraction 

studies (Ohtaki and Radnai, 1993) suggest a coordi

nation of either six or eight, with Cs-O distances of 

3.12 A or 2.95 A based upon two separate studies. 

The agreement of the aqueous simulation results with 

the experimental data provides some validation in the 

use of the selected forcefield parameters. 

We then tested the ability of the Cs' potential to 

model discrete water clusters with the number of wa

ter molecules in the inner solvation sphere ranging 

from one to six. Figure 3 shows the Cs' ion being 

solvated by six water molecules in octahedral coordi

nation. These results were compared to those derived 

from quantum mechanical calculations and were per

formed using the molecular DMol 3 code (Delley, 

2000). These calculations were performed within the 

GGA-PW91 exchange and correlation functional, 

using a double numerical plus polarization (dnp) ba

sis set, with a fine numerical integration grid and an 

octahedral expansion of the multipole densities. The 

Cs+~Owater distances shown in Table 3 are fully con

verged with respect to the basis set and the numerical 

integration grid.

Table 2. A comparison of the kaolinite unit cell parameters as determined by experiment, DFT-GGA calcula

tions, and the forcefield parameters used in this study. 

Experiment DFT-GGA Forcefield 

(Neder et al., 1999) (this study) (this study) a (A) 5.154 5.199 5.169 

b (A) 8.942 9.018 8.923 c (A) 7.401 7.455 7.253 

a 91.690 91.730 91.680 _ 104.610 104.640 104.620 

y 89.820 89.850 90.460 

V (A 3

) 329.91 338.01 323.31

4-1
Figure 3. A graphical representation of a Cs' ion 

hydrated by six waters in an octahedral geometry. 

The sixth water molecule is hidden by the sphere 

representing the cesium ion.

16_ý

Table 3. Equilibrium Cs-Owater distances as the water coordination number increases calculated using den

sity functional theory and the forcefield described in Table 1. 

Cluster DFT-GGA Forcefield 

R(Cs-O) (A) R(Cs-O) (A) 

Cs'+ 1 H20 3.014 2.982 

Cs÷ + 2 H20 3.062 2.990 

Cs÷ + 3 H20 3.073 3.003 

Cs' + 4 H2 O 3.101 3.019 

Cs+ + 6 HO 3.180 3.059 

Table 4. Calculation of the cohesive energy of the Cs+CI complex as a function of ionic separation using DFT 

and the forcefield parameters described in Table 1. 

R(Cs+-CI") DFT-GGA Forcefield 

(A) (kcal/mol) (kcal/mol) 

2.5 -95.201 15.3228 

3.0 -112.160 -95.699 

3.5 -102.921 -93.084 

4.0 -90.436 -82.869 

4.5 -79.571 -73.879 

5.0 -71.018 -66.528

4-2

The Cs+-Owat distances calculated using the force

field compare favorably with those calculated using 

density functional theory. The CS-Owater distances 

increase with the water coordination number, be

cause, as more water molecules are packed around a 

Cs' ion, they will strive to have a larger Owater-Owater 

separation. This result further validates the Cs' po

tential as being suitable for modeling Cs' in an aque

ous solution. 

4.3 Modeling of the CsC1" Com

plex 

As a final validation of the Cs' and Cl forcefield pa

rameters, we calculated the equilibrium bond length 

for a Cs'C1 complex and the cohesive energy as a 

function of Cs'-C1 separation. These results were 

compared to those derived from quantum mechanical 

calculations using the Molecular D Mol 3 code.
These calculations were performed within the GGA

PW91 exchange and correlation functional using a 

dnp-basis set, with a fine numerical integration grid and an octahedral expansion of the multipole densi

ties. The R(Cs*-C[) distances shown in Table 4 are 

fully converged with respect to the basis set and the 

numerical integration grid. 

The equilibrium separation determined using the 

forcefield (3.153 A), is slightly larger than that ob

tained from DFT (2.933 A). The energetics derived 

from the forcefield show smaller binding energies and 

a far stronger repulsive nature at close distances. 

However, in the region of interest (3.0 A and greater), 

the forcefield reproduces the DFT results reasonably 

well and is applicable to the modeling of Cs÷Cr com

plexes in an aqueous solution. 

4.4 Modeling of Cesium and 

Strontium Sorption onto the 

Kaolinite (001) Surface 

4.4.1 Preliminary Simulations 

4.4.1.1 Kaolinite 

Large-scale molecular dynamics simulations to simu

late the sorption of Cs' and Sr2÷ ions on the (001) 

surfaces of kaolinite and montmorillonite, were per

formed using the forcefield parameters presented in 

Table I. Preliminary simulations were performed on 

a four-processor SGI workstation using the OFF mo

lecular dynamics code (Molecular Simulations Inc., 

1999). Production calculations were completed using 

an MPI-parallelized version of the GULP 1.3 code 

(Gale, 1997) on SGI Origin2000 supercomputers at 

Sandia National Laboratories. The GULP 1.3 code 

scales near-linearly with the number of processors 

and is able to take advantage of the high floating

point performance processors available on the SGI 

Origin2000. The different codes and platforms used 

in these simulations represent an effective use of 

computational facilities. The development of the clay 

forcefield and preliminary simulations were run on a 

Unix workstation with one to four processors. The 

preliminary simulations, including the one noted be

low, required anywhere from approximately three to 

four days to two weeks to complete, while the pro

duction simulations on the massively parallel super

computers were typically completed in a day. Unfortunately, most of the computational cost of the latter 

simulations was due to queuing and preprocessing 

requirements. 

Figure 4 provides the simulation cell used in the pre

liminary analysis of the sorption of a I M CsCl solu

tion onto two unique basal surfaces of kaolinite. This 

system was used as a basis for determining the suit

ability of the clay forcefield for the sorption simula

tions and for troubleshooting the approach prior to 

performing the production simulations. The periodic 

simulation cell was created as a 4 x 4 x 2 supercell of 

the kaolinite unit cell. The CsC1 solution was first 

equilibrated, separate from the kaolinite, using mo

lecular dynamics and an NVT ensemble for 100 psec. 

The solution and kaolinite substrate, comprised of 

approximately 2400 atoms, were then equilibrated 

together at 300 K for 100 psec and then for an addi

tional equilibrated period of 100 psec. 

The equilibrated trajectory from the preliminary CsCI 

simulations was analyzed by running a density and 

probability algorithm to extract the positions for each 

atom during the 100 psec period. The density profile 

normal to the kaolinite surfaces for the equilibrated 

part of the 1 M CsCl simulation with kaolinite, is pre

sented in Figure 5. The density plot indicates the 

relatively diffuse and weak structure of the solution 

components relative to the sharp density peaks for the 

crystalline components. The only exception to this is 

the CI peak that occurs at the aluminol surface and 

which represents the strongly sorbed inner-sphere 

complex. In general, the Cl- sorbs to the positively

charged aluminol surface while the Cs÷ sorbs to the 

negatively-charged siloxane surface. The density 

profile supports this statement; however, some Cs' 

also sorbs with the cr as an ion pair onto the alumi

nol surface, as indicated by less intense Cs peak at 

aluminol surface, and which is not unexpected for the 

fairly concentrated CsCI solution. Cesium also oc

curs at the siloxane surface as a small peak at 2 A 

(inner-sphere complex) and as a more diffuse peak at 

5 A (outer-sphere complex). An additional peak for 

CI" occurs at 12 A, indicating an adsorbed outer

sphere complex on the aluminol surface. Hydrogen 

and oxygen profiles are typically diffuse across the 

profile, except for the near interface region where 

surface-controlled hydrogen bonding occurs to pro

vide some structure. Water molecule orientations 

occur with the water dipole directing the hydrogens to 

the siloxane surface and oxygens to the aluminol sur

face.

4-3
.....

'I 

* -J

Mostly 

Cs+ 

Mostly 

CI0 

T 

0 

T

Figure 4. Representative configuration for the preliminary analysis of the equilibrated Cs÷ sorption onto the 

(001) basal surfaces of kaolinite. The smaller Cs' spheres sorb primarily to the tetrahedral siloxane surface, 

while the larger CI" spheres coordinate to the aluminol surface. Note that the simulation is periodic in three 

dimensions and that both basal surfaces of kaolinite are available to the CsCl solution. The very small 

spheres immediately above the kaolinite surface are the hydrogen atoms associated with the aluminol groups.

cc

* - Crystalline-> ISi

AI

0

CsuSLI:•olution Mr 0 

30 H 0 

20 CI 

101 

1 CS 

0 5 10 15 20

0 

ML i

25 30

z/ A

Figure 5. Density profile of the equilibrated 1 M CsCl solution interacting with the two basal surfaces of kao

linite based on the accumulation of 100 psec of trajectory analysis from the preliminary molecular dynamics 

simulation. The crystalline kaolinite is represented by the narrow and structured peaks on the right, and the 

solution species are denoted by the more diffuse peaks on the left. The strongly sorbing Cl" is indicated by the 

distinct peak at the aluminol interface.

4-4

70 

60

r, )

Cl) 

CO 

x 

-9

40

Q) 

CC

Ar,
A more detailed density plot of the equilibrated kao

linite-CsCI solution is provided in Figures 6 and 7 

where atomic density maps exhibit the contours for 

the occurrence of inner and outer-sphere complexes 

on the kaolinite surfaces. Figure 6 presents the re

sults for the aluminol surface where the Cr sorption sites clearly occur as inner-sphere complexes above 

the aluminum vacancies. Cesium also exists as an 

inner-sphere complex on this surface, but only when 

coordinated to the sorbed Cl'. Hydrogen atoms of the 

aluminol groups are displaced appropriately, based on 

their proximity to the Cs' and CF. The outer-sphere density map shows more diffuse contours for the CI 

complexes indicating the greater mobility of the 

sorbed species in this region above the kaolinite alu

minol surface. One of the inner-sphere CF" exhibits 

some transitory behavior into the outer-sphere region. Again, Cs' only occurs when associated with the 

sorbed Cl species as an ion pair complex. Sorption site densities derived from the density maps are 0.27 

Cs+/nm 2 and 1.0 Cl-/nm2 . 

The corresponding density maps for the siloxane sur

face of kaolinite are dominated by the sorption of 

Cs'; there is no occurrence of CI" sorption as either an 

inner or outer-sphere complex on this surface. Cesium sorption is observed as both inner and outer

sphere complexes (Figure 7). As an inner-sphere 

complex, Cst primarily binds at the center of the hex

agonal ring of the siloxane surface. Also, it is occa

sionally observed directly above the Sit44 tetrahedra 

of the ring. The outer-sphere density map exhibits 

significantly diffuse contours for Cs', indicating a less 

strongly bound complex. As observed for CF on the 

aluminol surface, Cs' exhibits dynamic motion and 

goes through transitions between inner and outer

sphere complexes. The sorption site density is 0.65 

Cst/nm2 on the siloxane surface of kaolinite. 

Distribution coefficient can be derived from the 

equilibrated simulation by a comparison of sorbed 

Cs' with the Cs' still remaining in solution. Ap proximately 65 % of the Cs' is sorbed as either an 

inner or outer-sphere complex. This corresponds to a 

Kd value of 0.0024, which if extrapolated to more 

applicable solution concentrations of 0.001 M CsC1, 

corresponds to a Kd of 24 ml/g. This value is in rea

sonable agreement with the experimental value of 12 

m1/g (Brady et al., 1998). More details on the 

method used to calculate K, values are provided be

low.

Inner Sphere Complex < 3A Outer Sphere Complex 3 to 6 A

10 

a / A 10 a/ A

Figure 6. Density maps of the equilibrated 1 M CsC1 solution interacting with the aluminol basal surface of kaolinite, based on the accumulation of 100 psec of trajectory analysis from the preliminary molecular dy namics simulation. The left figure indicates the inner-sphere region within 3 A of the aluminol surface and the outer-sphere region between 3 A to 6 A is on the right. The hexagonal structure of the underlying alumi

num and oxygen sublattices are provided in each density map.

15 20

4-5
Inner Sphere Complex < 3A 

35 r. v U 

o 1 e o 0 

0I 0 

30 

25. ~ ~ 6 

200 1 400 

a a 

1 0 0 6 0 

Cs* at siloxane cavity 

0 _ 

0 5 10' 15 20 

a / A

35

0ýx

25 

20

15

Outer Sphere Complex 3 to 6 A 

* oc a oo V 

0 6' 0 

eG *

0 0 0. 

90000 ' 

,6' 9 6' 0 

*00 **L0

e a & a 

50 

Cs* OS complexes dominate 

10 0 

0 

0 5 10 15 20 

a/ A

Figure 7. Density maps of the equilibrated 1 M CsCl solution interacting with the siloxane basal surface of 

kaolinite based on the accumulation of 100 psec of trajectory analysis from the preliminary molecular dynam

ics simulation. The left figure indicates the inner-sphere region within 3 A of the siloxane surface, and the 

outer-sphere (OS) region between 3 A to 6 A is on the right. The hexagonal structure of the underlying sili

con and oxygen sublattice is provided in each density map.

4.4.2 Production Simulations 

4.4.2.1 Kaolinite 

A two-layer slab model of the kaolinite (001) surface 

was created using an 8 x 4 x 2 supercell of the kaolin

ite unit cell. The interlayer distance was increased by 

32 A and populated with 1670 H20 molecules. This 

system was then equilibrated using an NVT molecular 

dynamics scheme for 20 psec with a time-step of 

0.001 psec. After reaching equilibrium, 15 Cs' and 

15 CY ions (0.5 M CsC1) were inserted halfway along 

the interlayer distance and allowed to relax using an 

energy minimization scheme. A series of ten con

stant-NVT molecular dynamics runs were completed 

using a time-step of 0.001 psec and a simulation 

length of 500 psec. In all of these simulations, the 

dominant mechanism for Cs' sorption onto the alumi

nol surface was through the formation of Cs+C1 

complexes and their ssubsequent binding to the outer 

hydroxyl groups (Figure 8). .A series of ten simula

tions of a -0.1M CsCl solution interacting with the 

(001) kaolinite surface were undertaken by reducing 

the number of Cs' and Cl ions from 15 to 4. With 

the ionic strength of the solution reduced, the forma

tion of Cs+C1- complexes was decreased. At these

lower ionic strengths, the Cs' ions did not tend to 

sorb to the aluminol surface, but instead bound adja

cent to the center of the six-membered ring silica tet

rahedral sheet in a transitory inner-sphere complex 

(Figure 9). The Cl ions tend to form outer-sphere 

complexes that were attracted to the (001) aluminol 

plane. 

In order to simulate a 0.01 M CsCl solution interact

ing with the (001) kaolinite surface, a two-layer slab 

model was constructed using a 24 x 12 x 2 supercell 

of the kaolinite unit cell. The interlayer distance was 

increased by 32 A and populated with 15,030 H20 

molecules. This large system (approximately 64,700 

atoms) was then equilibrated using molecular dynam

ics and a constant-NVT ensemble for 20 psec with a 

time step of 0.001 psec. After reaching equilibrium, 

4 Cs' and 4 Cl- ions (0.01 M CsC1) were inserted 

halfway along the interlayer distance and allowed to 

relax using an energy minimization scheme. A series 

of five NVT molecular dynamics runs were com

pleted using a timestep of 0.001 psec and a simulation 

length of 250 psec. The equilibrium minimum energy 

configuration was again determined to be with the Cs' 

ions bound adjacent to the center of the six

membered ring silica tetrahedral sheet (Figure

4-6

4•
- - N,

0

Figure 8. A graphical representation of a 0.5 M 

CsCI solution in steady-state equilibrium with a 

(001) kaolinite surface. At high ionic strengths, 

the formation of CsC1 complexes occurs which can 

lead to Cs' sorbing to the basal aluminol plane. 

10). This inner-sphere complex can be transitory with the Cs' ion desorbing from the surface and then 

either resorbing or going back into solution (Figure 11). 

4.4.2.2 Montmorillonite 

Similar molecular dynamics simulations were com

pleted for the montmorillonite basal surfaces interact

ing with either a Cs' or Sr2+ solution. Only one TOT 

sheet of the clay was used to represent the substrate 

surfaces. Due to the charged nature of the montmoril

lonite layers, only the metal cations were incorpo

rated into the solution. This ensured a neutral system where the negative charge of the tetrahedrally- and 

octahedrally-substituted montmorillonite, is balanced 

by the solvated metal cations. The large systems (ap proximately 65,000 atoms), were first equilibrated 

using molecular dynamics at constant-NVT ensemble 

for 20 psec with a time step of 0.001 psec. After 

reaching equilibrium, the metal cation solute was 

inserted halfway along the interlayer distance and 

allowed to relax using an energy minimization

CsCI 

complex

0

0 '~Ž~N ~ *N2'NN*. ~"~' 7.

A'M IAAA~ 

Figure 9. A graphical representation of a 0.1 M 

CsCI solution in steady-state equilibrium with a 

(001) kaolinite surface. This is the energetically

preferred configuration with the Cs÷ ions sorbed 

to the basal tetrahedral plane and the CI" ions 

loosely sorbed to the basal aluminol plane. 

scheme. A series of five NVT molecular dynamics 

runs were then completed, using a time-step of 0.001 

psec and a simulation length of 250 psec. 

Representative equilibrium results for the simulations 

of the Cs' and Sr2

÷ systems interacting with mont

morillonite are presented in Figures 12 and 13, re

spectively. The Cs -montmorillonite example exhib

its a highly sorbing environment, with most of the Cs' 

ions occurring as inner-sphere complexes and a few 

Cs' ions remain solvated and isolated in the solution. 

The strongest sorption of the Cs' occurs near mont

morillonite sites associated with tetrahedral (Al for 

Si) and ochahedral (Mg for Al) substitutions where 

there is an excess of negative charge. The transitory 

nature of the sorption, especially for outer-sphere Cs' 

complexes, may lead to the clay substrate and perhaps 

later, resorption of the metal. This behavior under

scores the need to apply statistically-sound methods 

for evaluating the sorption state of the system over 

the equilibrated simulation trajectory and for multiple simulations. The Sr 2 -montmorillonite example (Fig

ure 13) indicates an equilibrated system in marked

4-7

17 177-7 

ý.M MWIW - -, ___ - -_ - ý'ý _N ý 't. -- ,: '. k' ý I - -,- . -
Figure 10. A graphical representation of a 0.01 M CsCl solution in steady-state equilibrium with a (001) kao

linite surface. This is the energetically-preferred configuration with the Cs' ions sorbed to the basal tetrahe

dral plane and the Cl ions loosely sorbed to the basal aluminol plane. 

AN J 

0 

T 

SCs1 

- C S+ 

Figure 11. Details of the 0.01 M CsCl solution in steady-state equilibrium with a (001) kaolinite surface. Note 

that the Cs' ions on the left are sorbed as an inner-sphere complex onto the basal aluminol plane while the 

one on the right is in the process of transitioning from an inner-sphere complex to a loosely sorbed outer

sphere complex. This second Cs' is desorbed from the surface 10 psec after this dynamic configuration.

4-8
* <:2j..

- *7"

Figure 12. A graphical representation of a 0.1 M 

Cs' solution in steady-state equilibrium with a 

(001) montmorillonite surface. This is the ener

getically preferred configuration with the Cs' ions 

sorbed to the charged sites on the top and bottom 

tetrahedral planes. 

contrast to those observed for the Cs' simulations. 

The Sr 2

+ ion is rarely observed binding to the mont

morillonite basal surface. When sorption does occur, 

it is often transitory and usually exists as an outer

sphere complex. All of the Sr 2

+-montmorillonite 

simulations indicate a preference for the metal ion to 

be fully solvated by water molecules and remain in 

solution.

A AV: 

Figure 13. A graphical representation of a 0.1 M 

Sr 2

+ solution in steady-state equilibrium with a 

(001) montmorillonite surface. This is the ener

getically-preferred configuration with the Sr 2+ 

ions sorbed to the charged sites on the top and 

bottom tetrahedral planes. 

A set of molecular dynamics simulations was com

pleted to examine the competition presented in Figure 

14 where the Cs' occurs fully sorbed to the basal sur

face, as mostly inner-sphere complexes and the Sr 2

, 

remains in solution. Strontium ions often approach the basal surface during the simulation but never be

yond a brief existence (up to 10 psec) as an outer

sphere complex. Statistical fluctuations in energy and 

ion trajectories will impact the evolution of these sys

tems, especially when Cs' remains at a sorption site 

preventing the near approach of a Sr2

, ion.

t:...... . ""• > •:••N ' • " " ' ' . • -"• ' ... .. ..

. - • .. S.. 

Figure 14. A graphical representation of 8 Cs' and 8 Sr2

' ions in steady-state equilibrium with a (001) mont morillonite surface. This is the equilibrated configuration showing the competition between two different cations, with the Cs' ions are preferentially sorbed to the montmorillonite surface while the Sr 2

, ions remain 

solvated.

4-9

... .. ...

i

VVVT W V

il

411
5. CALCULATION OF DISTRIBUTION COEFFICIENTS

5.1 Kaolinite 

The distribution coefficient (Kd) is a measure of the 

ability of a material to sorb contaminants, and is de

fined as: 

Kd = 1 i C f _ 

where C, is the initial concentration, Cf is the final 

concentration, V1 is the volume of the liquid in ml and 

MA. is the mass of the sorbent in grams. The higher 

the Kd, the greater the ability of a material to remove 

a contaminant from solution. 

In order to evaluate the distribution coefficient of the 

kaolinite (001) surface to sorb Cs' ions using the re

sults of the computer simulations, the condition of 

sorption must first be defined. Based upon observing 

the sorption and desorption of Cs' in the molecular 

dynamics simulations, we defined sorption as the Cs' 

ion being within 5 A of the plane of the aluminol 

groups or the oxygens on the silicate tetrahedral 

plane. Using this condition, the trajectory files for 

the molecular dynamics simulations were evaluated 

and the Kd's were calculated. We determined a Kd of 

8 ± 3 mrl/g for the 0.1 M CsCI solution and 10 ± 4 

ml/g for the 0.01 M CsC1 solution. These results 

agree reasonably well with the experimental value of 

12 mug (Brady et al., 1998). It is important to note 

that the calculated Kd values assume that the entire 

contribution to Cs' sorption can be attributed to the 

basal surfaces of the kaolinite and montmorillonite. 

In part this is a good assumption for smectite clays such as montmorillonite, but may be less valid for 

kaolinite where metal sorption onto edge sites can be 

significant (Brady et al., 1996). Therefore, Kd values 

obtained from the kaolinite simulations should be 

considered as limiting cases, assuming the basal sur

faces are the only sorbing surfaces. More detailed 

simulations incorporating the exposed edge sites

with their variable protonation state and therefore a 

function of pH-would be required to assess the rela

tive role of edge versus basal sorption in kaolinite. 

5.2 Montmorillonite 

The significance of edge sites is less important for the 

Kd values derived from the montmorillonite simula

tions, where basal surfaces are the predominant sorption surface (Pabalan et al., 1998). For the case of 

our montmorillonite simulations, we have determined 

a Kd of 293 ± 30 ml/g to sorb Cs' ions from a 0.1 M 

CsCI solution and a Kd of 11 ± 4 ml/g to sorb Sr 2

1 

ions from a 0.1 M SrC12 solution. 

The competitive sorption simulation of Cs' and Sr2, 

onto montmorillonite indicates the strong dominance 

of Cs' over Sr 2

+ with Kd values that are in general 

agreement with the pure Cs' solution. 

5.3 Energetics and Kd Values 

A series of energy calculations was performed to ex

amine the energy of various cations sorbed on kaolin

ite, montmorillonite and pyrophyllite, relative to the 

energies of the same cations being fully solvated by 

water molecules in solution. Pyrophyllite 

(Si 8Al4020(OH) 4) represents a proxy for montmoril

lonite where there are no tetrahedral or octahedral 

substitutions; effectively this substrate represents 

weak montmorillonite sorption sites. The strong sorption sites on montmorillonite are basal sites asso

ciated with a tetrahedral substitution and nearby octa

hedral substitutions. Energy minimization calcula

tions were performed for each configuration. Results 

are presented for Cs', Na÷, and Sr 2

, ions in Table 5. 

The reported energy values represent the optimum 

energy differences for the sorption site or hydration 

state. 

The significance of these results is in comparing the 

energy values for each sorption site with the water 

hydration (solvation) energy for the chosen metal ion. 

For example, Cs' indicates a fairly small solvation 

energy of -60 kcal/mol whereas the kaolinite and 

strong montmorillonite sorption energies are some

what greater indicating a preference for Cs' to be 

sorbed onto the clay. Similarly, Na÷ exhibits a pref

erence for only the strong montmorillonite site com

pared to remaining in solution. In contrast, Sr 2

+ 

shows an extremely large hydration energy of -315 

kcal/mol that is related to the double ionic charge and 

relatively small ionic radius (Na' = 0.97 A; Sr 2 ÷ 

1.12 A; Cs÷ = 1.67 A). All of the clay substrates ex

hibit Sr 2

+ sorption energies that are similar to, or 

much less, than this large hydration energy. There

fore, it is not surprising to see limited sorption, if any 

at all, of Sr 2

÷ in any of the previous molecular dynam

ics simulations.

5-1
Table 5. Comparison of energy for various ions based on water solvation or clay surface sorption environ

ment. The strong montmorillonite site corresponds to ion sorption near the tetrahedral aluminum and octa

hedral magnesium substitutions. The weak montmorillonite site represents a pyrophyllite-like smectite hav

ing no tetrahedral or octahedral substitutions. 

Energy (kcal/mol) Cs+ Na' Sr2+ 

Water solvation -60 -110 -315 

Kaolinite -74 -105 -275 

Strong montmorillonite -106 -142 -333 

Weak montmorillonite -37 -50 -163

5-2
6. CONCLUSIONS

Large-scale molecular dynamics were completed to 

model the sorption of Cs÷ and Sr 2 ions onto the (001) 

planes of kaolinite and montmorillonite. Contrary to 

earlier simulations, (Cygan et al., 1998), it was found 

that at high ionic strength cesium tends only to sorb to 

the aluminol plane through the formation of CsC1 ion 

pairs. At lower concentrations, the Cs÷ ions tend to 

sorb adjacent to the center of the six-membered ring 

silicate tetrahedral sheet in a transitory inner-sphere 

complex. Cesium sorption onto montmorillonite is

preferred over Sr 2

÷ sorption due to the relatively large 

hydration energy of Sr 2

÷ and relatively small sorption 

energy of Sr2'. We are able to calculate distribution 

coefficients by defining the condition of sorption and 

by completing a series of large-scale and long-term 

molecular dynamics simulations in order to provide 

the necessary statistics. The results of this study indi

cate a strong link between the atomistic simulations 

and the macroscopic sorption properties.

6-1
7. SUGGESTIONS FOR FURTHER WORK

It is difficult to predict stable Kd's for Cs' sorbing to 

the (001) planes of kaolinite because of the formation 

of transitory inner- and outer-sphere complexes. In

stead of trying to model the distribution coefficient, it 

seems worthwhile to evaluate the selectivity of a sur

face for the desired contaminant ion, and determine 

how that surface could be chemically or physically 

modified to enhance irreversible contaminant sorp

tion. 

The results of this study show that Cs' preferentially 

binds to kaolinite and montmorillonite surfaces over 

Sr2

' because the Cs' ion is poorly solvated by the 

water and migrates to the surface, while the Sr2

÷ is 

more strongly solvated by water and tends to remain 

in solution. We find that Cs' and Sr 2

+ bind more ef

fectively to montmorillonite than kaolinite due to the 

presence of negatively-charged sites due to AI÷3 sub

stituting for Si4 and Mg÷2 substituting for A]÷3

. It is 

possible to examine the selectivity for a wide variety of cations to sorb to these and other layered materials. 

For the case of montmorillonite and kaolinite, our 

preliminary modeling studies suggest the following 

order of preference for major sorbing cations: 

NH4 I>Cs÷>K*>Na÷>Li4 . This interesting result sug

gests that NH4 I will exchange for Cs' on natural 

clays. This may be the mechanism that is responsible 

for the phenomenon known as "cask weeping", as 

atmospheric water contains significant amounts of 

dissolved NH4

4 ions. When a nuclear transport cask 

is covered with rain or dew, the Cs' ions that are 

bound to natural clays stuck to the surface of the cask

desorb and are replaced by NH.' ions. A possible 

treatment to minimize this phenomenon would be to 

wash the cask with an ammonium solution prior to 

transport. 

During the course of this study, we have developed 

more accurate forcefields for modeling the solvation 

of complex ions. We are now able to model the sorp

tion of such radioactive ions as 1, Se0 4

2 , Tc0 4 , Ra2

*, uranium-carbonate complexes, as well as RCRA 

metals such as CrO4

2 and HAsO 4

2

". Ultimately, our 

molecular simulation approach, when incorporating 

more efficient statistical mechanics algorithms such 

as thermodynamic integration, will allow the calcula

tion of free energies for these complex systems. 

Priorities for extensions of this study include the mo

lecular dynamics simulation of more complex ra

dionuclides such as the uranyl and uranyl-carbonate 

complexes. The latter are of significance due to their 

stability in solution at higher pH conditions leading to 

a reduction in the sorption of uranium. As noted 

above, the radioactive anions should be examined to 

derive relative rank ordering of sorption strengths. It 

would also be important to examine the influence of 

edge sites under various protonated states in order to 

assess the relative behavior of edge and basal sites for 

sorption. This last priority would be a difficult task 

due to the variation of acid pK values for the numer

ous cleavage surfaces and clay substitutions. How

ever, a limited set of simulations may be helpful in 

providing a basis for analyzing sorption trends.

7-1
8. REFERENCES

Adams, J.M. (1983) "Hydrogen atom positions in 

kaolinite by neutron profile refinement." Clays and 

Clay Minerals, 31, 352-356. 

Bish, D.L. (1993) "Rietveld refinement of the kaolin

ite structure at 1.5K". Clays and Clay Minerals, 41, 

738-744. 

Bish, D.L. and Von Dreele, R.B. (1989) "Rietveld 

refinement of non-hydrous atomic positions in kaolin

ite." Clays and Clay Minerals, 37, 289-296. 

Brady, P.V., Cygan, R.T., and Nagy, K.L. (1996) 

"Molecular controls on kaolinite surface charge." 

Journal of Colloid and Interface Science, 183, 356

364. 

Brady, P.V., Cygan, R.T., and Nagy, K.L. (1998) 

"Surface charge and metal adsorption to kaolinite." 

In "Adsorption of Metals by Geomedia " (E.A. Jenne, 

Ed.), Academic Press, San Diego. 

Cygan, R.T., Nagy, K.L., and Brady, P.V. (1998) 

"Molecular models of cesium sorption on kaolinite." 

In "Adsorption of Metals by Geomedia " (E.A. Jenne, 

Ed.), Academic Press, San Diego. 

Cygan, R.T., Liang, J.J., Kalinichev, A.G. (2002) 

"Molecular models of hydroxide, oxyhydroxide, and 

clay phases, and the development of a general force

field." Journal of Physical Chemistry, in press. 

Delley, B. (2000) "From molecules to solids using the 

DMol3 approach." Journal of Chemical Physics, 

113, 7756-7764. 

Gale, J.D. (1997) GULP: "A computer program for 

the symmetry adapted simulation of solids." Philoso

phical Magazine B, 73, 3-19. 

Giese, R.F. (1982) "Theoretical studies of the kaolin 

minerals: Electrostatic calculations." Bulletin of 

Mineralogy, 105, 417-424. 

Hess, A.C. and Saunders, V.R. (1992) "Periodic ab 

initio Hartree-Fock calculations of the low-symmetry 

mineral kaolinite." Journal of Physical Chemistry, 

96, 4367-4374. 

Hobbs, J.D., Cygan, R.T., Nagy, K.L., Schultz, P.A., 

and Sears, M.P. (1997) "Ab initio energy minimization of the kaolinite crystal structure." American 

Mineralogist, 82, 657-662. 

Hohenberg, P. and Kohn, W. (1964) "Inhomogeneous 

electron gas." Physical Review, 136, 864-861. 

Kim, Y., Cygan, R.T., and Kirkpatrick, R.J. (1996a) "'33Cs NMR and XPS investigation of cesium ad

sorbed onto clay minerals and related phases." Geo

chimica et Cosmochimica Acta, 60, 1041-1052. 

Kim, Y., Kirkpatrick, R.J., and Cygan, R.T. (1996b) "133Cs NMR study of cesium on the surfaces of kao

linite and illite." Geochimica et Cosmochimica Acta, 60, 4059-4074. 

Kim, Y. and Kirkpatrick, R.J. (1997) ,

23

Na and 1

33

Cs 

NMR study of cation adsorption on mineral surfaces: 

local environments, dynamics, and effects of mixed 

cations." Geochimica et Cosmochimica Acta, 61, 

5199-5208. 

Kim, Y., Kirkpatrick, R.J. (1998) "NMR T1 relaxa

tion study of 1

3 3

Cs and 23

Na adsorbed onto illite." 

American Mineralogist, 83, 661-665. 

Kohn, W., and Sham, L.J. (1965) "Self-consistent 

equations including exchange and correlation ef

fects." Physical Review A, 140, 1133-1138. 

Komarneni, S. (1978) "Cesium sorption and desorp

tion behavior of kaolinites." Soil Science Society of 

America Journal, 42, 531-532. 

Krumhansl, J.L., Brady, P.V., and Anderson, H.L. 

(2001) "Reactive barriers for 137Cs retention.." Jour

nal of Contaminant Hydrology, 47, 233-240. 

Michaelian, K.H., Lapides, I., Lahav, N., Yariv, S., and Brodsky, I. (1998a) "Infrared study of the inter

calation of kaolinite by caesium bromide and caesium 

iodide." Journal of Colloid and Interface Science, 

204, 389-393. 

Michaelian, K.H., Zhang, S.L., Yariv, S., and Lapi

des, I. (1998b) "Low-frequency raman spectra of 

kaolinite alkali halide complexes." Applied Clay 

Science, 13, 233-243. 

Molecular Simulations Inc. (1999) "Cerius-2 4.0 

User Guide, Forcefield-Based Simulations." Molecu

lar Simulations Inc., San Diego.

8-1
Monkhorst, H.J., and Pack, J.D. (1976) "Special 

points for Brillouin-zone integrations." Physical Re

view B, 13, 5188-5192. 

Nagy, K.L., Cygan, R.T., and Brady, P.V. (1995) 

"Kaolinite morphology: AFM observations, model 

predictions, and adsorption site density." Clay and 

Minerals Society Program Abstracts, 32, 95. 

Pabalan, R.T., Turner, D.R., Bertetti, P., and Prikryl, 

J. D. (1998) "Uraniumvl sorption onto selected min

eral surfaces." In Adsorption of Metals by Geomedia 

(E.A. Jenne, Ed.), Academic Press, San Diego. 

Neder, R.B., Burghammer, M., Grasl, T., Schulz, H., 

Brain, A.. and Fiedler, S. (1999) "Refinement of the 

kaolinite structure from single-crystal synchrotron 

data.". Clays and Clay Minerals, 47, 487-494. 

Ohtaki, H. and Radnai, T. (1993) "Structure and dy

namics of hydrated ions." Chemical Reviews, 93, 

1157-1204. 

Perdew, J.P. and Wang, Y. (1992) "Accurate and 

simple analytic expression of the electron-gas correla

tion energy." Physical Review B, 45, 13244-13249. 

Pickett, W.E. (1989) "Pseudopotential methods in 

condensed matter applications." Computational 

Physics Reports, 9, 115-197. 

Saiers, J.E. and Hornberger, G.M. (1999) "The influ

ence of ionic strength on the facilitated transport of 

cesium by kaolinite colloids." Water Resources Re

search, 35, 1713-1727. 

Smith, D.E. and Dang, L.X. (1994a) "Computer 

simulations of NaCl association in polarizable water." 

Journal of Chemical Physics, 100, 3757-3766.

Smith, D.E. and Dang, L.X. (1994b) "Computer 

simulations of cesium-water clusters: Do ion-water 

clusters form gas-phase clathrates?" Journal of 

Chemical Physics, 101, 7873-7881. 

Telemen, 0., Jonsson, B., and Engstrom, S. (1987) 

"A molecular dynamics simulation of a water model 

with intramolecular degrees of freedom." Molecular 

Physics, 60, 193-203. 

Tsipursky, S.I. and Drits, V.A. (1984) "The distribu

tion of octahedral cations in the 2:1 layers of dioc

tahdral smectites studied by electron diffraction." 

Clay Minerals, 19, 177-193 

Vanderbilt, D. (1990) "Soft self-consistent pseudopo

tentials in a generalized eigenvalue formalism." 

Physical Review B, 41, 1990. 

Weiss, C.A., Jr., Kirkpatrick, R.J., Altaner, S.P. 

(1990) "The structural environments of cations ad

sorbed onto clays: 133CS variable-temperature MAS 

NMR spectroscopy study of hectorite." Geochimica 

et Cosmochimica Acta, 54, 1655-1669. 

Westrich, H.R., Cygan, R.T., Brady, P.V., Nagy K.L., 

Anderson, H.L., Kim, Y., and Kirkpatrick, R.J. 

(1995) "The sorption behavior of Cs and Cd onto 

oxide and clay surfaces." In "Proceeding of the 

Waste Management Conference, WM'95. ",24-4. Fo

lio Infobase. 

Young, R.A. and Hewatt, A.W. (1988) "Verification 

of the triclinic crystal structure of kaolinite." Clays 

and Clay Minerals, 36, 225-232.

8-2
NRC FORM 335 U.S. NUCLEAR REGULATORY COMMISSION 1. REPORT NUMBER (2-89) 

(Assigned by NRC. Add Vol., Supp.. Rev., 

NRCM 1102, 

3201, 3202 BIBLIOGRAPHIC andAssignedub m Nubrs, idd any.)up, e. DATA SHEET and Addendum Numbers, if any.) 

(See instructions on the reverse) 

2. TITLE AND SUBTITLE 

NUREGICR-6757 

Large-Scale Molecular Dynamics Simulations of Metal Sorption onto the Basal Surfaces of Clay Minerals 

3. DATE REPORT PUBLISHED MONTH YEAR 

January 2002 

4. FIN OR GRANT NUMBER 

5. AUTHOR(S)W811 5. AUTOR(S)6. TYPE OF REPORT 

D.M. Teter and R.T. Cygan 

7. PERIOD COVERED (Inclusive Dates) 

December 1999 to November 2001 8. PERFORMING ORGANIZATION - NAME AND ADDRESS (if NRC, provide Division, Office or Region, U.S. Nuclear Regulatory Commission, and mailing address, if contractor, 

provide name and mailing address) 

Geochemistry Section 

Sandia National Laboratories 

Albuquerque, NM 87185-0750 

9. SPONSORING ORGANIZATION - NAME AND ADDRESS (If NRC, type ASame as above'^ if confraclto, provide NRC Division, Office orRegion, U.S. Nuclear Regulatory Commission, 

and mailing address.) 

Division of Systems Analysis and Regulatory Effectiveness 

Office of Nuclear Regulatory Research 

U.S. Nuclear Regulatory Commission 

Washington, DC 20555-0001 

10. SUPPLEMENTARY NOTES 

11. ABSTRACT (200 words or less) 

Large-scale 

related 

molecular dynamics computer simulations have been completed to evaluate the sorption of Cs+, Sr2+, and other 

concentration 

ions onto the (001) basal planes of kaolinite and montmorillonite. Simulations were undertaken as a function of metal ion 

ions bind to both 

to evaluate 

kaolinite sorption 

and montmorillonite, 

mechanisms and to determine distribution coefficients (Kd's). The simulations show that Cs+ 

M), 

while Sr2+ ions tend to remain in solution. At high ionic strengths (greater than 0.5 

where 

the 

the 

fundamental 

Clmechanism for Cs+ sorption to the aluminol plane of kaolinite is through the formation of CsCI complexes 

sorb 

ion is loosely bound to the outer hydroxyl groups. At lower concentrations (0.1 M and 0.01 M), the Cs+ ions do not to the aluminol plane, but prefer to bind adjacent to the center of the six-membered ring silica tetrahedral sheet as an inner-sphere 

Kd of 0 for kaolinite 

complex. We have determined a Kd of 10 ± 4 ml/g for kaolinite to sorb Cs+ ions from a 0.01M CsCI solution, and a 

293 ± 30 ml/g 

to sorb Sr2+ ions from a 0.01 M SrCl2 solution. For the case of montmorillonite, we have determined a Kd of 

The montmorillonite 

to sorb Cs+ ions from a 0.1M CsCI solution, and a Kd of 11± 4 mVg to sorb Sr2+ ions from a 0.1M SrCl2 solution. 

created by 

sorbs the cations more effectively than kaolinite because it has negatively-charged sorption sites that are 

strongly 

isomorphic substitutions on both tetrahedral and octahedral sites. Calculations show that the Sr2+ ion is solvated more than Cs+ and is not sorbed by the kaolinite or montmorillonite surfaces. Directions for further work are outlined and prioritized. 

12. KEY WORDS/DESCRIPTORS (List words or phrases that will assist researchers in locating the report.) 13. AVAILABILITY STATEMENT 

sorption, adsorption, desorption, Cs sorption, Sr sorption, sorption on clay, radionuclide transport unlimited 

14. SECURITY CLASSIFICATION 

(This Page) 

unclassified 

(This Report) 

unclassified 

15. NUMBER OF PAGES 

16. PRICE 

NRC FORM 335 (2-89)

I
Federal Recycling Program
NUREG/CR-6757 Large-Scale Molecular Dynamics Simulations of Metal Sorption onto the February 2002

Basal Surfaces of Clay Minerals

UNITED STATES 

NUCLEAR REGULATORY COMMISSION 

WASHINGTON, DC 20555-0001 

OFFICIAL BUSINESS 

PENALTY FOR PRIVATE USE, $300