Langmuir 1998, 14, 5959-5967
5959
Molecular Computer Simulations of the Swelling Properties and Interlayer Structure of Cesium Montmorillonite David E. Smith Department of Chemistry and Biochemistry, MSC 3C, New Mexico State University, P.O. Box 30001, Las Cruces, New Mexico 88003-8001 Received January 5, 1998. In Final Form: July 1, 1998 The crystalline swelling properties and interlayer structure of a cesium montmorillonite clay were investigated using molecular computer simulations. Two classes of dry-clay structures, proposed previously to explain X-ray diffraction and NMR experiments, were identified using Monte Carlo annealing calculations. Hydrated clays with water contents ranging from 0.044 to 0.440 gH2O/gclay were investigated using molecular dynamics simulations. Layer spacings calculated as a function of water content were found to be similar to experimental swelling curves, showing a distinct plateau at the monolayer-hydrate spacing. Hydration energies were calculated as a function of water content and expressed in three complementary forms. The immersion energy form was found to be most useful, revealing an apparently global minimum in the swelling-coordinate energy that corresponds to the monolayer hydrate. This is in agreement with experimental measurements and may help clarify the energetic origins of discrete, crystalline swelling processes in clay minerals. For two- and four-layer hydrates, cesium ions readily formed two different types of inner-sphere complexes with the clay surface. Ions associated with negatively charged tetrahedral substitution sites formed exclusively inner-sphere complexes and occupied hexagonal cavities adjacent to the substitutions. Other cesium ions occupied both inner-sphere and outer-sphere configurations with roughly equal probability. The ease with which cesium associates with the clay surface may be responsible for the formation of monolayer hydrates in cesium-substituted clays and for selective binding of cesium to many clay minerals.
1. Introduction Significant releases of radioactive materials into the environment have occurred over the past fifty years through weapons testing, accidents such as the one at Chernobyl, and the intentional or inadvertent discharge of nuclear waste into soils and groundwater.1 Cesium137 is a major contributor to the overall radioactivity of released waste and an important target element in environmental remediation efforts.2 The transport and biological uptake of radioactive cesium in the environment is governed by many factors including adsorption onto ion-exchanging soil constituents.3,4 Montmorillonites and other 2:1 clay minerals have among the largest cationexchange capacities of common soil components,5 and they play a correspondingly significant role in determining the fate of radioactive cesium in the environment. Molecular theories that treat ion binding to clay minerals and that account for groundwater and clay composition could assist in predicting the transport properties of radioactive ions in the environment. The binding strengths of cesium to a variety of clay minerals have been determined experimentally.6-15 Without ex(1) U.S. Dept. of Energy. Closing the circle on the Splitting of the Atom; U.S. Government Printing Office: Washington, DC, 1995. (2) Cornell, R. M. J. Radioanal. Nucl. Chem. Artic. 1993, 171, 483. (3) Coughtrey, P. J.; Thorne, M. C. Radionuclide Distribution and Transport in Terrestrial and Aquatic Ecosystems; A. A. Balkema: Rotterdam, 1983; Vol. 1. (4) Staunton, S. Eur. J. Soil Sci. 1994, 45, 409. (5) Sposito, G. The Chemistry of Soils; Oxford: New York, 1989. (6) Eliason, J. R. Am. Mineral. 1966, 51, 324. (7) Gast, R. G. Soil Sci. Soc. Am. Proc. 1969, 33, 37. (8) Sawhney, B. L. Clays Clay Miner. 1972, 20, 93. (9) Gast, R. G. Soil Sci. Soc. Am. Proc. 1972, 36, 14. (10) Maes, A.; Cremers, A. J. Chem. Soc., Faraday Trans. 1978, 74, 1234. (11) Xu, S.; Harsh, J. B. Soil Sci. Soc. Am. J. 1990, 54, 1596.
ception, cesium has been found to bind to clays most selectively among the alkali-metal cations. Cesium binding strengths are also known to increase with clay layer charge10 and to depend significantly upon layer-charge location.11,12 Existing theories of ion binding, reviewed recently,16,17 are based either on simple electrostatic descriptions of ion-water and ion-clay interactions or upon empirical hard and soft, acid and base concepts. Continued development of molecular theories of ion binding requires knowledge of clay layer-charge distributions and of the structure of the ion-clay complexes. This information is, in general, unavailable experimentally but may be investigated conveniently using molecular simulation methods18 and clay models that have been developed recently.19-25 One objective of our research is to use simulation methods to correlate ion binding strengths with (12) Xu, S.; Harsh, J. B. Clays Clay Miner. 1992, 40, 567. (13) Oscarson, D. W.; Hume, H. B.; King, F. Clays Clay Miner. 1994, 42, 731. (14) Adeleye, S. A.; Clay, P. G.; Oladipo, M. O. A. J. Mater. Sci. 1994, 29, 954. (15) Ohnuki, T.; Kozai, N. Radiochim. Acta 1994, 66/67, 327. (16) Laudelout, H. In Chemistry of Clays and Clay Minerals; Newman, A., Ed.; Wiley: New York, 1987. (17) Xu, S.; Harsh, J. B. Soil Sci. Soc. Am. J. 1990, 54, 357. (18) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K., 1987. (19) Bleam, W. F. Rev. Geophys. 1993, 31, 51. (20) Skipper, N. T.; Refson, K.; McConnell, J. D. C. Clay Miner. 1989, 24, 411. (21) Skipper, N. T.; Refson, K.; McConnell, J. D. C. J. Chem. Phys. 1991, 94, 7434. (22) Skipper, N. T.; Chang, F.-R. C.; Sposito, G. Clays Clay Miner. 1995, 43, 285. (23) Delville, A. Langmuir 1991, 7, 547. (24) Delville, A. J. Phys. Chem. 1995, 99, 2033. (25) Teppen, B. J.; Rasmussen, K.; Bertsch, P. M.; Miller, D. M.; Scha¨fer, L. J. Phys. Chem. B. 1997, 101, 1579.
S0743-7463(98)00015-8 CCC: $15.00 © 1998 American Chemical Society Published on Web 09/12/1998
5960 Langmuir, Vol. 14, No. 20, 1998
Smith
the structure (including clay layer-charge distribution) of the ion-clay complexes. Ion binding to clay minerals is influenced by changes in interlayer hydration that occur upon ion exchange.10 Cesium adsorption on hydrated montmorillonites is known to lead to dehydration of the interlayer region and formation of a monolayer hydrate structure.8,26-28 This is presumably related to the weak hydration energy of cesium relative to other alkali- and alkaline-earth-metal ions. Understanding the molecular origins of crystalline swelling (including the formation of discrete-hydrate structures) and the relationship of swelling to ion exchange in clay minerals is a second objective of our research. In this paper, we report results of Monte Carlo (MC) and molecular dynamics (MD) computer simulations of the crystalline swelling properties and interlayer structure of a cesium-montmorillonite clay. The relationship of these properties to calculated ion-exchange thermodynamics is presently being investigated and will be described elsewhere. The swelling investigations have provided clear evidence for an energetic origin of discrete, crystalline swelling in clay minerals and for the formation of mixed-layer structures in partially hydrated Csmontmorillonite. Structural investigations have revealed a tendency for cesium ions to associate with the clay surface. This surface affinity may be correlated with the preferential formation of monolayer hydrates for cesiumsubstituted clays and with the selective binding of cesium to many clay minerals. 2. Computational Methods 2.1. Interaction Potentials: Water and Cesium. Water is represented in these simulations by the “extended” simple point charge (SPC/E) model of Berendsen and co-workers.29 While conceptually and computationally simple, the SPC/E model effectively reproduces the properties of both bulk water and simple ionic solutions.30,31 Of popular water models, it is one of the few that includes the self-polarization energy correction that is appropriate for models parametrized with a static polarization. (It is this correction that distinguishes SPC/E water from its SPC parent.32) It was also shown recently31 that SPC/E water behaves similarly to the polarizable water model developed by Dang and coworkers.33,34 This similarity will facilitate future comparison simulations addressing the role of water polarizability in hydrated clay behavior. SPC/E water is a rigid, three-site model with an O-H bond length of 1.0 Å and a tetrahedral H-O-H bond angle of 109.47°. There are fixed charges located on each atomic site, yielding a permanent dipole moment of 2.35 D and a Lennard-Jones (LJ) interaction
σ ULJ ) 4 r
12
σ r
6
[( ) ( ) ]
Table 1. Solution Simulation Parameters element O H Cs
σ, Å
, kJ‚mol-1
q, e
3.166 0.000 3.830
0.650 0.000 0.418
-0.848 0.424 1.000
centered on the oxygen atom. The values of the water parameters are displayed in Table 1. Cesium ions are also represented by a fixed charge plus LJ interaction potential. The LJ parameters, given in Table 1, were adjusted to give agreement with experimental hydration data. A 150 ps simulation of one cesium chloride ion pair solvated by 214 water molecules was performed using methods described elsewhere.31 The calculated solvation energy of -624 ( 18 kJ‚mol-1 for the ion pair is in reasonable agreement with the experimental enthalpy of -655 kJ‚mol-1.35 The calculated coordination number (8.4) and Cs+-O peak positions (3.09 Å) also agree well with the experimental values of 8 and 2.95-3.15 Å, respectively.36,37 The potential energy of the interlayer solution due to water-water and ion-water interactions is therefore given by
U)
1
∑∑ 2 i j
{ [( ) ( ) ] }
′ 4
σij
12
-
rij
σij rij
6
+
qi q j rij
(2)
where the j summation is taken over all atoms except those in the same molecule as atom i. The ion-oxygen cross interaction terms in the LJ potential are generated using the usual combining rules:
1 σij ) (σi + σj) 2
(3)
ij ) (ij)1/2
(4)
2.2. Interaction Potentials: Clay. Montmorillonites are moderately charged, swelling clay minerals composed of an inner sheet, in which aluminum atoms are octahedrally coordinated with oxygen, sandwiched by two outer sheets, in which silicon atoms are tetrahedrally coordinated with oxygen. The resulting “2:1” sheets combine together via van der Waals and electrostatic interactions to form layered materials. The charge on the clay results most commonly from isomorphic substitution of magnesium for aluminum in the octahedral layer, or of aluminum for silicon in the tetrahedral layers. Exchangeable positive ions reside in the clay interlayer region and balance the negative layer charge. The Cs-substituted, Wyoming-type montmorillonite modeled here has a unit-cell formula of
Cs0.75[Si7.75Al0.25](Al3.5Mg0.5)O20(OH)4 (1)
(26) Mooney, R. W.; Keenan, A. G.; Wood, L. A. J. Am. Chem. Soc. 1952, 74, 1371. (27) Calvet, R. Ann. Agron. 1973, 24, 77. (28) Be´rend, I.; Cases, J.-M.; Francois, M.; Uriot, J.-P.; Michot, L.; Masion, A.; Thomas, F. Clays Clay Miner. 1995, 43, 324. (29) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (30) Watanabe, K.; Klein, M. L. Chem. Phys. 1989, 131, 157. (31) Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 100, 3757. (32) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981. (33) Dang, L. X. J. Chem. Phys. 1992, 97, 2659. (34) Caldwall, J.; Dang, L. X.; Kollman, P. A. J. Am. Chem. Soc. 1990, 112, 9144.
This clay has a net layer charge of -0.75 e per unit cell, with two-thirds of the charge located in the octahedral layer, and interlayer cesium ions to balance the layer charge. This particular clay was chosen to facilitate comparison with previous simulation studies of alkalimetal-substituted montmorillonites38-41 and because of (35) Friedman, H. L.; Krischnan, C. V. In Water, A Comprehensive Treatise; Franks, F., Ed.; Plenum: New York, 1973; Vol. 6. (36) Ohtomo, N.; Arakawa, K. Bull. Chem. Soc. Jpn. 1979, 52, 2755. (37) Pa´linka´s, G.; Radnai, T.; Hajdu, H. Z. Naturforsch. A 1980, 35, 107. (38) Skipper, N. T.; Sposito, G.; Chang, F.-R. C. Clays Clay Miner. 1995, 43, 294. (39) Boek, E. S.; Coveney, P. V.; Skipper, N. T. J. Am. Chem. Soc. 1995, 117, 12608.
Structures of Cesium Montmorillonite
Langmuir, Vol. 14, No. 20, 1998 5961
Table 2. Clay Simulation Parameters layer tetrahedral apical octahedral
element O Si Al O O H Al Mg
σ, Å
, kJ‚mol-1
q, e
3.166 1.840 1.840 3.166 3.166 0.000 0.000 0.000
0.650 13.18 13.18 0.650 0.650 0.000 0.000 0.000
-0.800 1.200 0.200 -1.000 -1.424 0.424 3.000 2.000
the abundance and importance of montmorillonite clays in the environment. The clay simulation model used here is derived directly from empirical models introduced by Skipper and coworkers.20-22 In these models, the clay is represented as a rigid lattice with atom coordinates taken from X-ray diffraction studies of talc, an uncharged 2:1 clay mineral.42 Atomic charges and van der Waals parameters for the oxygen atoms were assigned on the basis of a proposed similarity between water oxygen and clay basal oxygen behavior. These models have been used to investigate Li-, Na-, and K-montmorillonites38-41 and Na-vermiculite38 and have been generally successful at representing known structural and thermodynamic properties of these systems. The success of these previous studies lends credibility to the use of rigid lattice models and the current parametrization. Nevertheless, future investigations to identify sensitivity of the simulations to clay relaxation (particularly for hydroxyl groups) and charge parametrization are desirable, and some such studies are currently in progress. The existing clay models have been modified in this work to achieve consistency with the SPC/E water model and to provide uniform LJ-type parametrization of the van der Waals interactions. The structural OH groups in the octahedral layer are represented with the SPC/E water parameters. The LJ parameters on the SiO4 oxygens have also been changed, although the oxygen charges have not been modified. Parameters for the tetrahedral-layer Si and Al atoms were adjusted to yield Si-O and Al-O interactions similar to those reported by Skipper et al.22 Both LJ and charge parameters for all clay atoms are displayed in Table 2. Atom positions are identical to those given by Skipper et al.22 and are not reproduced here. 2.3. Simulation Methods. The simulation system consists of a periodically replicated, rectangular box with two montmorillonite layers each containing eight clay unit cells and having x- and y-axis lengths of 21.12 Å by 18.28 Å, respectively. The two clay layers are mirror images. This allows for simultaneous overlap of hexagonal (technically ditrigonal) cavities in both interlayer regions and is required to correctly model dry-clay structures. Registry shifts of the two clay sheets are sampled using a MC algorithm described below. The z-axis is perpendicular to the clay sheets and varies in length with water content and applied pressure. While the clay layer spacing is thus variable, the two interlayer spacings within the simulation cell are constrained to be equal. A snapshot of the simulation cell for a monolayer hydrate of Cs-montmorillonite is displayed in Figure 1. The hydrated clay simulations were performed using the SOLVENT MD code described elsewhere.43 The (40) Chang, F.-R. C.; Skipper, N. T.; Sposito, G. Langmuir 1995, 11, 2734. (41) Chang, F.-R. C.; Skipper, N. T.; Sposito, G. Langmuir 1997, 13, 2074. (42) Brindley, G. W.; Brown, G. Crystal Structures of Clay Minerals and Their X-ray Identification; Mineralogical Society: London, 1980. (43) Smith, D. E.; Haymet, A. D. J. J. Chem. Phys. 1992, 96, 8450.
Figure 1. Image of the periodically replicated simulation cell for a monolayer hydrate of Cs-montmorillonite.
equations of motion were integrated using a 1.5 fs time step and the velocity-Verlet algorithm.44 Ewald summation methods were used to treat Coulombic interactions,18,45 and a minimum image convention was used for all short-ranged terms. Water bond lengths were constrained using the RATTLE46 variant of the SHAKE47 algorithm. All simulations were performed at a constant temperature of 298 K using a “weak-coupling” algorithm and a coupling time constant of 500 fs.48 The clay sheets were rigid and immovable except for swelling and registry degrees of freedom. The z-axis length was allowed to fluctuate through a “weak-coupling” pressure algorithm and a target pressure of 1.0 bar.48 Clay registry shifts were accomplished with a MC algorithm in which a small, random registry move is attempted and accepted or rejected using the usual Metropolis criterion.18 Three separate starting configurations were generated for a hydrated clay system containing 192 water molecules, or a water content of 0.264 gH2O/gclay. The first configuration was created by inserting two 96-molecule, bulk water slabs into the interlayer regions of a clay simulation cell with an initially wide interlayer spacing. The cesium ions were initially placed at random points on the clay surface. Following equilibration of this system, a second starting configuration was generated by applying a biasing potential to the clay-ion interaction and forcing ions to the center of the interlayer regions. The bias was then removed and the new system equilibrated. A third starting configuration was generated with an initially random ion and water distribution. All three starting configurations generated statistically equivalent layer spacings, energies, and structures following equilibration of ∼100 ps. Configurations for the systems with smaller water contents were generated by removing water from an equilibrated, 192-molecule configuration, followed by equilibration for between 30 and 90 ps depending upon the number of water molecules removed. A configuration (44) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. J. Chem. Phys. 1982, 76, 637. (45) de Leeuw, S. W.; Perram, J. W.; Smith, E. R. Proc. R. Soc. London, A 1980, 373, 27. (46) Andersen, H. C. J. Comput. Phys. 1983, 52, 24. (47) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (48) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684.
5962 Langmuir, Vol. 14, No. 20, 1998
Smith
Table 3. Energiesa and Clay Layer Spacingsb for Symmetric (SYM) and Asymmetric (ASYM) Configurations, As Determined from Dry-Clay MC Simulations at Three Annealing Rates anneal rate no. of SYM no. of ASYM 〈U〉SYM 〈U〉ASYM 〈z〉SYM 〈z〉ASYM
slow
mod
rapid
10 0 0.0
10 0 3.3
10.64
10.68
5 5 4.7 15.9 10.72 11.57
a The energy, 〈U〉, is reported as the relative average energy per cesium ion in units of kJ‚mol-1. b The clay layer spacing, 〈z〉, is the average spacing in Å obtained for the given class of configuration.
for a 320-molecule, 0.440 gH2O/gclay system was generated with an initially random water and ion distribution and equilibrated for 80 ps. In all cases, no significant drift in layer spacing or energy was observed following equilibration. Averages were collected on each system for a total of 300 ps. Uncertainties were calculated using a block averaging technique in which each final average was broken into 10 30 ps blocks. The reported uncertainties correspond to 95% confidence limits in the mean calculated using standard statistical techniques. This analysis assumes, without independent confirmation, that the 30 ps data blocks are uncorrelated. The determination of equilibrium structure and energy for the dry clay at small layer spacing is complicated by the presence of deep ion-clay potential energy minima. Ambient temperature MD simulations do not effectively sample possible ion-clay configurations. In addition, MD algorithms are not well-suited for this system since each ion interacts nearly independently with the clay lattice, and kinetic energy transfer from ion to ion is limited. To address these issues, we have implemented a MC-based simulated annealing method for sampling dry-clay structures. The method is a simple application of the Metropolis MC methodology.18 Each MC cycle consists of one attempted move for each of the 12 interlayer ions plus one attempted lattice registry shift and one attempted z-axis move. All moves are accepted or rejected on the basis of the usual temperature-dependent criteria. The range of movement is chosen such that the acceptance ratio near ambient temperature is between 0.1 and 0.5. 3. Results and Discussion 3.1. Dry Clay. The properties of the dry clay were investigated in order to provide energy and layer spacing reference points for the hydrated clay systems. Sets of simulated annealing calculations were performed at three different annealing rates. Each annealing run consisted of 1.0 × 104 MC equilibration cycles at a temperature of 2000 K followed by “rapid”, “moderate”, or “slow” cooling to a temperature of 298 K over a period of 1.7 × 103, 8.5 × 103, or 4.3 × 104 cycles, respectively. The 298 K systems were then sampled for an additional 5.0 × 103 cycles in order to collect averages of energy and clay layer spacing. Lengthy sampling is not required for the collection of averages, since only local, small motions are sampled at ambient temperature for the dry clay. Results are presented in Table 3. Two classes of dry-clay configurations were obtained. The first class is a symmetric (SYM) structure where the cesium ions are nestled within overlapping hexagonal cavities of neighboring clay sheets. Each ion is 12-fold
Figure 2. Illustration of the two classes of dry-clay structures. Only the basal oxygen atoms and one set of interlayer cesium ions are shown. The left and right figures show the symmetricclass (SYM) and asymmetric-class (ASYM) structures, respectively. The oxygen atoms in the closer basal plane are light in color, while those behind are darker.
coordinated with oxygen atoms. An example of this class of structure is illustrated in Figure 2. The second class of structure is asymmetric (ASYM), with the cesium ions occupying a hexagonal hole on one clay sheet and a triangular site above a silicon atom on the other sheet. Each ion is therefore 9-fold coordinated with oxygen atoms, as illustrated in Figure 2. Note that both structure classes contain sites that would yield 6-fold-coordinated cesium ions but that none of these sites is occupied in the dry clays. The existence of two types of energetically competitive dry structures for Cs-clays has been proposed on the basis of X-ray diffraction (XRD) studies of Cs-vermiculite49 and Cs-montmorillonite28 and on 133Cs NMR studies of Cs-hectorite.50 The XRD studies yielded layer spacings for the two structures of 10.6 and 11.5 Å for Cs-vermiculite and 10.7 and 11.5 Å for Cs-montmorillonite. These are in good agreement with the calculated SYM and ASYM spacings. The simulations indicate an energetic preference for the SYM structure, with ASYM configurations produced only in the rapid annealing runs. Nevertheless, the identification of two types of stable structures in the simulations provides qualitative support for the experimentally based proposal of competitive SYM and ASYM configurations. The calculated energies may be sensitive to the clay layer charge distribution and to details of the interaction potential model. Additional studies of the relative energy distributions of the two structure classes are therefore being undertaken. 3.2. Swelling Curves. Simulations were performed on 14 hydrated systems with between 32 and 320 water molecules, corresponding to water contents ranging from 0.044 to 0.440 gH2O/gclay. Equilibrium layer spacings, clay hydration energies, and ion and water structures were determined from 300 ps simulations for each system. The correlation of layer spacing and clay water content for Cs-montmorillonites has been determined from XRD experiments.26,27 The layer spacings calculated from the simulations are presented in Table 4 and plotted in Figure 3 along with the experimental data. The dry-clay spacing is for the minimum-energy SYM structure determined from the annealing runs. The experimental swelling curves displayed in Figure 3 both show a layer-spacing plateau of ∼12.5 Å, corresponding to a monolayer hydrate. The simulations reproduce this plateau, with slightly smaller spacing. At water contents of above 0.12 gH2O/ gclay the simulated swelling curve increases monotonically. In contrast, the experimental data of Mooney et al. show (49) Laperche, V.; Lambert, J. F.; Prost, R.; Fripiat, J. J. J. Phys. Chem. 1990, 94, 8821. (50) Weiss, C. A., Jr.; Kirkpatrick, R. J.; Altaner, S. P. Geochim. Cosmochim. Acta 1990, 54, 1655.
Structures of Cesium Montmorillonite
Langmuir, Vol. 14, No. 20, 1998 5963
Figure 3. Cs-montmorillonite swelling curves from the simulations (filled circles) and from experimental measurements of Calvet27 (diamonds) and Mooney et al.26 (open circles). Error bars give 95% confidence limits.
Figure 4. Cs-montmorillonite hydration energy, ∆U, as a function of clay water content. Error bars give 95% confidence limits. The solid line corresponds to the mean interaction energy of bulk SPC/E water.
Table 4. Simulated Layer Spacings and Hydration Energies for Cesium-Montmorillonite
useful information on individual-layer swelling properties for comparison with the bulk-clay behavior measured experimentally. The presence of a plateau in the swelling curves suggests that the expansion of individual Cs-montmorillonite layers to form stable monolayer hydrates occurs by a “pop and fill” mechanism. The first few adsorbed water molecules pop the layer open, and subsequent water molecules fill up the remaining interlayer space. The similarity of the shape of the curves in Figure 3 is evidence that the simulations reproduce well the general clay swelling behavior, although the small difference in plateau spacing may indicate a need for refinement of the interaction potentials. The similarity also suggests that mixed-layer hydrate formation is not required to account for the experimental XRD spacings. This should not, however, be used to conclude that mixed-layer hydrates do not form. 3.3. Hydration Thermodynamics. In this section, we describe the calculation of three thermodynamic quantitiessthe hydration energy, the isosteric heat of adsorption, and the immersion energysfrom average potential energies of the hydrated clay systems. Although all three of these quantities contain identical information, they may be related to different experimental measurements. Of the three, the immersion energy is found to yield the clearest representation of the clay hydration behavior, at least with regard to stepwise crystalline swelling and mixed-layer hydrate formation. The clay hydration energy may be defined as
H2O content, gH2O/gclay
spacing, Å
∆U, kJ‚mol-1
0.000 0.044 0.066 0.088 0.099 0.110 0.121 0.132 0.154 0.176 0.198 0.220 0.242 0.264 0.440
10.53 ( 0.01 12.05 ( 0.02 12.24 ( 0.02 12.31 ( 0.02 12.38 ( 0.02 12.50 ( 0.02 12.68 ( 0.06 13.08 ( 0.04 14.22 ( 0.06 14.80 ( 0.02 15.19 ( 0.04 15.74 ( 0.03 16.33 ( 0.03 17.11 ( 0.03 22.01 ( 0.02
-49.68 ( 0.25 -49.54 ( 0.29 -50.21 ( 0.15 -50.49 ( 0.15 -50.60 ( 0.17 -50.30 ( 0.30 -49.10 ( 0.15 -46.06 ( 0.14 -45.31 ( 0.09 -45.52 ( 0.08 -45.43 ( 0.06 -45.12 ( 0.09 -44.45 ( 0.08 -43.18 ( 0.05
a fixed layer spacing at the higher water contents. This is most likely due to adsorption of water on external clay surfaces. The comparison of experimental and simulated swelling curves is complicated by several factors related to simulation scope and constraints. For example, the simulations are designed to model only interlayer water, while experimentally measured water contents include contributions from external-surface and pore-space hydration. The latter possibly leads to capillary-condensation hysteresis in adsorption isotherms and swelling curves.51 In addition, it is known that clay swelling often proceeds through the formation of mixed-layer hydrates, with the XRD peak position giving an averaged layer spacing of the randomly interstratified structure.42 The experimental dry-clay spacing of 11.2 Å is possibly also an average of SYM and ASYM values, as discussed above. In contrast, the simulation constraint of uniform layer spacing and water content precludes the formation of mixed-layer hydrates. These issues suggest that a direct comparison of experimental and simulated swelling curves is unlikely to provide more than a qualitative validation of the simulation methods. Instead, the simulations provide (51) Adamson, A. W.; Gast, A. P. Physical Chemistry of Surfaces, 6th ed.; Wiley: New York, 1997.
∆U )
〈U(N)〉 - 〈U(0)〉 N
(5)
where N is the number of interlayer water molecules adsorbed on the clay and 〈U(N)〉 and 〈U(0)〉 are the average hydrated- and dry-clay potential energies, respectively. The hydration energy reflects not only the energy of adsorbed water but also changes in the clay-clay and ion-clay energies upon hydration. The calculated hydration energy for Cs-montmorillonite is displayed in Figure 4. The predominant features in the plot are an approach toward the average bulk water energy as the water content
5964 Langmuir, Vol. 14, No. 20, 1998
Smith
Figure 5. Isosteric heat of adsorption, qst, for Cs-montmorillonite. Error bars give 95% confidence limits. The solid line corresponds to the mean vaporization enthalpy of bulk SPC/E water.
is increased and two shallow minima near stable monolayer (∼0.1 gH2O/gclay) and two-layer (∼0.2 gH2O/gclay) hydrates. The hydration energy is an average or “integral” quantity. The related integral hydration enthalpy may be determined experimentally either by calorimetry or through integration of water adsorption isotherms at two different temperatures.51 In practice, the integration is difficult due to uncertainties in the low-pressure isotherm data. Calorimetry experiments are usually reported as immersion enthalpies rather than as integral hydration enthalpies. The enthalpic driving force for hydration of a clay at a given water content is not a function of the integral hydration enthalpy but rather depends on the partial molar or “differential” enthalpy. This quantity corresponds to the change in enthalpy of a clay of specified water content, per mole of water, upon adsorption of an additional infinitesimal amount of water. The negative of the differential hydration enthalpy is called the isosteric heat of adsorption and denoted qst. It may be calculated from the temperature dependence of the experimental water adsorption isotherm using the Clausius-Clapeyron equation,51
ln P (∂ ∂T )
qst ) RT2
Γ
(6)
where Γ is the water surface concentration that in turn is proportional to the water content. The value of qst may be determined from a simulation through a finitedifference derivative of hydrated-clay energies,
h) ) qst(N
〈U(N)〉 - 〈U(N′)〉 + RT N - N′
(7)
where N h is the average of N and N′. The factor of RT is needed to convert between energy and enthalpy for this vaporization process. The results, displayed in Figure 5, show oscillations about the bulk water vaporization enthalpy. A positive deviation from the bulk value is evidence of an enthalpic driving force for further hydration upon equilibration of the clay with bulk water. A negative deviation indicates, in contrast, a driving force for dehydration. The points where the heat of adsorption
Figure 6. Cs-montmorillonite immersion energy. Error bars give 95% confidence limits. The reference hydration state is defined as 0.12 gH2O/gclay, the immersion energy minimum.
crosses the bulk vaporization enthalpy with negative slope, at water contents of 0.12 and 0.23 gH2O/gclay, correspond to stable hydration minima. The enthalpy minima evident from the calculated isoteric heats of adsorption are more clearly revealed through a plot of the clay immersion energy, Q, defined as
Q ) 〈U(N)〉 - 〈U(N°)〉 - (N - N°)Ubulk
(8)
where 〈U(N°)〉 is the average energy of the clay at a reference hydration state N°, to be specified, and Ubulk ) -41.4 kJ‚mol-1 is the mean interaction energy of bulk SPC/E water. The immersion energy is related to qst through
∂Q (∂N )
qst ) -
T,V
+ RT - Ubulk
(9)
It is analogous to experimental immersion enthalpies, varying only in the small difference between enthalpy and energy and in the reference state, taken experimentally to be the fully immersed clay. If the reference state is chosen to be the dry clay, the immersion energy and hydration energy are related by
Qd ) N[∆U(N) - Ubulk]
(10)
where the d subscript indicates the dry-clay reference state. The primary difference between Qd and ∆U is the factor of N; the hydration energy is expressed per mole of water, while the immersion energy is expressed (conventionally) per gram of clay. The calculated immersion energy for Cs-montmorillonite is displayed in Figure 6. The reference hydration state is taken to be a water content of 0.12 gH2O/gclay, corresponding to the first immersion energy minimum. The curve shows two clear minima, a two-layer hydrate local minimum and a monolayer hydrate overall minimum. The data point at 0.45 gH2O/gclay provides evidence that the monolayer hydrate is a true global minimum. Proof of this assignment, however, would require additional calculations at high water contents. The calculated thermodynamic functions displayed in Figures 4-6 all imply that there is an energetic stabilization of the discrete monolayer and two-layer hydrates in
Structures of Cesium Montmorillonite
Cs-montmorillonite. This is most clearly evident in Figure 6, with the immersion energy curve interpreted as a potential energy surface for the layer-spacing coordinate of Cs-montmorillonite in contact with bulk water. Upon immersion, an initially dry clay would adsorb water to form a monolayer hydrate, while an initially expanded clay (prepared, for example, by exchange of Na-montmorillonite with Cs+) would dehydrate to the monolayer state. The dehydration would occur as long as a mechanism for surmounting the swelling-coordinate barrier between the two-layer and monolayer hydrates could be found. It is known experimentally that montmorillonites form integer-layer hydrates and that Cs-montmorillonites form monolayer hydrates in water.42 The immersion energy data are consistent with these observations and provide computational verification of discrete swelling process in clay minerals. Note, however, that entropic effects are not evaluated from the simulations and may perturb quantitative conclusions. Several previous simulations performed at constant water chemical potential have addressed clay swelling through measurements of disjoining pressures as a function of clay layer spacing.52-54 Oscillations in the pressure are a signature of a discrete swelling process. Delville observed pressure oscillations for mica when the potassium counterions were constrained to occupy clay surface sites but not when the ions were disssociated from the surface.53 For Na-montmorillonite, he observed a monotonic decrease in the swelling pressure with increasing layer spacing.52 In contrast, Karaborni et al. measured significant pressure oscillations for a similar Na-montmorillonite.54 The variability in these simulation results may be due to differences in interaction potentials or to computational difficulties associated with constant chemical potential simulation methods. Our work demonstrates that immersion energy calculations can provide a complementary and efficient method for investigations of discrete swelling. Although clay hydration thermodynamics may be determined from isotherm and calorimetry experiments, little data exist for Cs-montmorillonite systems. Even if data were available, comparison of experimental and simulated thermodynamics is complicated by the possible formation of mixed-layer hydrates and by external surface contributions to experimental measurements. As with swelling curves, the simulations provide information that is complementary to, rather than directly comparable with, experimental data. For example, the pronounced oscillations of qst displayed in Figure 5 are not reproduced in typical experimental measurements.55,56 The oscillations reveal single-layer behavior that is averaged out in bulk clays through layer-layer equilibration and mixed-layer hydrate formation. Berend et al.28 recently reported a Cs-montmorillonite immersion enthalpy as a function of pre-immersion water vapor pressure. Although not reported directly as a function of water content, their data are qualitatively consistent with typical immersion enthalpies.57-59 For a clay that forms monolayer hydrates, the value of Q would (52) Delville, A. Langmuir 1992, 8, 1796. (53) Delville, A. J. Phys. Chem. 1993, 97, 9703. (54) Karaborni, S.; Smit, B.; Heidug, W.; Urai, J.; van Oort, E. Science 1996, 271, 1102. (55) Mooney, R. W.; Keenan, A. G.; Wood, L. A. J. Am. Chem. Soc. 1952, 74, 1367. (56) Slabaugh, W. H. J. Phys. Chem. 1959, 63, 436. (57) Keren, R.; Shainberg, I. Clays Clay Miner. 1975, 23, 193. (58) Zhang, Z. Z.; Low, P. F. J. Colloid Interface Sci. 1989, 133, 461. (59) Fu, M. H.; Zhang, Z. Z.; Low, P. F. Clays Clay Miner. 1990, 38, 485.
Langmuir, Vol. 14, No. 20, 1998 5965
Figure 7. Z-axis probability distribution functions for cesium ions (first column) and water oxygen atoms (second column). Labels indicate system water content in gH2O/gclay.
be expected to decrease roughly exponentially to near zero at the monolayer water content, followed by small changes in Q that are associated with external-surface hydration. In contrast, the simulated immersion energy shows a slight negative curvature in the region between the dry clay and the monolayer hydrate, and oscillations at higher water contents. The latter region is sampled in the simulations through “forced hydration” of the interlayer spaces to unphysically high water contents and should not be compared with experimental measurements. For submonolayer water contents, the variation between the experimentally observed decay and the simulation results is likely indicative of the formation of mixed-layer hydrates. While the simulations do not allow for mixedlayer states (each periodically replicated interlayer region has the same water content), the formation of mixed-layer hydrates in experimental systems is known to occur,42 being driven by clay layer inhomogeneity. Assuming that the most favorably hydrated interlayers fill first leads to a prediction of exponential-like decay, as observed in the experimental immersion enthalpy. The experimental dry-clay immersion enthalpy determined by Berend et al.28 is ∼20 J/gclay less than the calculated dry-clay immersion energy. The difference may be due to inaccuracies in the interaction potential model. Alternatively, thermal desorption experiments indicated that the “dry” Cs-montmorillonite used in the immersion enthalpy measurements actually contained approximately one water molecule per interlayer ion.28 This prehydration of ∼0.033 gH2O/gclay would lower the experimental immersion enthalpy significantly and could account for much of the discrepancy. For comparison, suppose the simulated “dry” clay actually contained 0.033 gH2O/gclay. This would lead to a shift in the simulated immersion energy to the left by 0.033 units, and near quantitative agreement between the simulated and experimental dry-clay results. 3.4. Interlayer Distribution of Water and Ions. The structure of the interlayer solution has been determined through calculations of the z-axis probability distribution functions for both water oxygen atoms and cesium ions. The resulting functions are displayed in Figure 7 for systems with water contents of 0.121, 0.242, and 0.440 gH2O/gclay. The first two water contents correspond to monolayer and two-layer hydrates, determined
5966 Langmuir, Vol. 14, No. 20, 1998
Smith
Figure 8. Graphical definition of triangular (T) and hexagonal (H) regions of the clay surface. The lines connect basal oxygen atoms.
from the immersion energy minima discussed above. The last corresponds approximately to a four-layer hydrate. The ion distributions for all hydrates indicate favorable inner-sphere interactions between cesium and the clay surface. The ion distribution is single-peaked for the monolayer hydrate and splits into two peaks for the twolayer hydrate, with a small population in the central region. For the four-layer hydrate, approximately twothirds of the interlayer ions remain in contact configurations with the clay surface. The contact peaks in the ion distribution functions for the two-layer and four-layer hydrates show structure that is related to the surface heterogeneity of the clay, as discussed below. The water distribution in the monolayer hydrate shows a central peak that is slightly split, perhaps indicating that the layer spacing is somewhat broader than expected from optimizing water-clay interactions alone. The splitting may also be associated with hydration of ions in different ion-clay configurations. The water distribution for the two- and four-layer hydrates shows significant peaks at the clay surface and roughly uniform central distributions. Favorable ion binding to the clay surface probably plays a key role in swelling inhibition and the preferential formation of monolayer hydrates in Cs-clays.60 The calculated ion distributions displayed in Figure 7 show enhanced ion binding to the clay relative to simulated distributions for Li- and Na-montmorillonites.39-41 The distributions are similar to those obtained for a clay substituted with potassium,39 a known swelling inhibitor.61 3.5. Surface Binding Sites. An analysis of the surface binding sites has been performed through a decomposition of the cesium z-axis distribution functions into contributions from the various surface sites. The clay surface is split conceptually into regions of triangular “T-sites” and hexagonal “H-sites” defined by connecting the clay basal oxygen atoms, as shown in Figure 8. The surface binding sites are further distinguished on the basis of the layer charge locations. Nonsubstituted T-sites (Tn) are situated above a tetrahedral-layer silicon atom, while substituted T-sites (Tt) are above a tetrahedral-layer aluminum atom. Each tetrahedral layer in the simulation contains one aluminum substitution and therefore one Tt-site and 31 Tn-sites. The H-sites are considered tetrahedrally sub(60) Sposito, G. In Clay-Water Interface and its Rheological Implications; Guven, N., Pollastro, R., Eds.; The Clay Minerals Society: Boulder, CO, 1992. (61) Denis, J. H.; Keall, M. J.; Hall, P. L.; Meeten, G. H. Clay Miner. 1991, 26, 255.
Figure 9. Cesium ion probability distribution function for the clay with a water content of 0.242 gH2O/gclay. The total distribution function (thick line) is decomposed into contributions from different sites, as defined in the text.
stituted (Ht) when the hexagonal cavity is adjacent to a tetrahedral-layer aluminum atom. There are therefore three Ht-sites for each clay surface in the simulations. Other H-sites are either nonsubstituted (Hn) or octahedrally substituted (Ho), depending on the location of the octahedral-layer magnesium atoms. The decomposed interlayer density profile for cesium was calculated for the one-, two-, and four-layer hydrates. Results for the two-layer hydrate are displayed in Figure 9. The site-specific distribution functions were calculated relative to one clay layer only. The two-layer hydrate was chosen for display in order to illustrate the binding to one clay surface without significant competing effects from the second surface. The innermost portion of the main peak in the ion distribution function corresponds to ions in Ht-sites. There are six of these sites, three on the top of the clay layer and three on the bottom. One of the Ht-sites also contains an octahedral substitution, and cesium ions in this site account for the inner shoulder on the Ht peak. The Tt peak, while of smaller area, occupies the same spatial region as the Ht peak. These two peaks together integrate to ∼2 at a spacing of 6.5 Å, indicating that each tetrahedral substitution is on average associated with a single interlayer cesium ion. In previous simulations, it was observed that sodium ions form strong, inner-sphere complexes with tetrahedral substitution sites, continuously occupying what is denoted here a Tt site.40 In contrast, cesium ions preferentially occupy the Ht sites. The inner-sphere peak also contains features associated with other surface sites, most notably the nonsubstituted Tn-sites. This indicates that cesium has a general affinity for the clay surface irrespective of the location of the octahedral substitutions. In fact, and perhaps surprisingly, these ions show a clear preference for Tn-sites over substituted, Ho-sites. The behavior of the two-layer hydrate, including the tendency for inner-sphere ions to occupy Ht- and Tn-sites, was qualitatively reproduced in both the monolayer and four-layer hydrates. This suggests that cesium ions occupy two distinct inner-sphere environments, one associated with tetrahedral-layer charges and one with other surface regions, even for the monolayer hydrate. Recent Cs NMR studies of several smectite clays suggest that there may
Structures of Cesium Montmorillonite
be multiple binding sites for interlayer cesium.50,62,63 The experiments indicate that multiple sites occur for both octahedrally substituted hectorite and tetrahedrally substituted saponite. It is unlikely, therefore, that the substitution-dependent sites identified in the simulations can explain the NMR results, and the structural identity of the experimentally determined sites is still in question. 4. Conclusions In this paper, we have described a computational investigation of the crystalline swelling properties and interlayer structure of a Cs-montmorillonite clay. The simulations allowed for a direct determination of interlayer swelling thermodynamics without competition from edge or external surface effects. The results indicated that there is an energetic origin for discrete swelling in clay minerals. Although this conclusion is not surprising, it does provide a unique computational verification of a crystalline swelling mechanism in clay minerals. The development of microscopic theories of crystalline swelling should benefit from this and similar investigations, especially those that include calculations of swelling free energies.53,54 The structures of the cesium-clay surface complexes were also determined. The simulations suggest (62) Weiss, C. A., Jr.; Kirkpatrick, R. J.; Altaner, S. P. Am. Mineral. 1990, 75, 970. (63) Kim, Y.; Cygan, R. T.; Kirkpatrick, R. J. Geochim. Cosmochim. Acta 1996, 60, 1041.
Langmuir, Vol. 14, No. 20, 1998 5967
that cesium readily forms two types of inner-sphere complexes with the clay, one associated with tetrahedral substitution sites and a second with other surface regions. No strong spatial association with octahedral substitution sites is observed. The formation of stable inner-sphere complexes with the clay surface is likely a cause of the preferential formation of monolayer hydrates in Csmontmorillonite. Investigations of the relationship between ion-exchange thermodynamics and clay swelling and structural properties is a natural extension of the work presented here. Calculations of sodium-cesium ion exchange free energies on a montmorillonite clay are currently being performed and will be reported shortly. It is hoped that this research will assist in the development and testing of theories of ion exchange on clay minerals, and in building an understanding of binding and transport mechanisms for radioactive cesium and other pollutants in the environment. Acknowledgment. This research was supported by the U.S. Department of Energy Environmental Management Science Program and by the donors of the Petroleum Research Fund, administered by the American Chemical Society. LA980015Z
