2734
Langmuir 1995,11, 2734-2741
Computer Simulation of Interlayer Molecular Structure in Sodium Montmorillonite Hydrates Fang-Ru Chou Chang Department of Environmental Science, Policy, and Management, 108 Hilgard Hall, University of California, Berkeley, California 94720-3110
N. T.Skipper Department of Physics and Astronomy, University College, Gower Street, London W C l E 6BT, United Kingdom
Garrison Sposito* Department of Environmental Science, Policy, and Management, 108 Hilgard Hall, University of California, Berkeley, California 94720-3110 Received November 4, 1994. I n Final Form: April 14, 1995@ Molecular structure in the interlayers of sodiumWyoming montmorillonite with one, two, or three water monolayers was investigated by Monte Carlo and molecular dynamics simulation, based on tested waterwater, water-counterion, clay-water, and clay-counterion potential functions. Calculated layer spacings and thermodynamic properties, as well as interlayer water configurations and interlayer-species selfdiffusion coefficients, were in good agreement with available experimental data. Inner-sphere surface complexesof Na+ with tetrahedral charge sites persisted in all hydrates, whereas the outer-sphere surface complexes with octahedral charge sites found in the one-layerhydrate dissociatedpartially into an incipient diffuse layer in the two- and three-layer hydrates. Even in the three-layer hydrate, interlayer water molecules did not replicate their structure or mobility in bulk aqueous solution. Improved diffraction and spectroscopic data on interlayer water are needed to guide future simulation studies.
Introduction The molecular structure of water and the distribution of counterions in the interlayer region of swelling 2:l layer type clay minerals such as smectites are the subject of intense current research.lr2 Perhaps the most controversial issue is the extent to which counterions solvate and thereby contribute to the swelling process in smectitewater systems. Claims have been made for both a dominant contribution3 and a negligible influence4 of counterion solvation, particularly in the case of the prototypical swelling smectite, Na-montmorillonite. Cases et ale5have reviewed the recent experimental literature directed to this issue and noted its problematic nature. Their own experiments, using a variety of classical surface chemical techniques, indicated that natural Na-montmorillonite forms one-, two-, and three-layer hydrates after adsorbing water vapor a t relative humidity up to 90%. Solvation of Na+ counterions was concluded by Cases et al.5 to occur a t relative humidities in the range 16-28%) followedby water vapor adsorption on the siloxane surface to complete a monolayer a t 50% relative humidity. The two-layer hydrate was only about three-fourths completed before a three-layer hydrate had begun to stabilize. Yamada et a1.6recently confirmed these relative humidity ranges for stability of the one- and two-layer hydrates of Abstract published in Advance A C S Abstracts, June 1, 1995. (1) Spdsito, G.; Prost, R. Chem. Rev. 1982,82,553. (2)Giiven, N. Molecular aspects of clay/water interactions. In ClayWater Interface and its Rheological Implications; Giiven, N., Pollastro, R. M., Eds.; The Clay Minerals Society: Boulder, CO, 1992;p 1. (3)Sposito, G. The diffuse-ion swarm near smectite particles suspended in 1:l electrolyte solutions: Modified Gouy-Chapman theory and quasicrystal formation. In Clay. Water Interface and its Rheological Implications; Giiven, N., Pollastro, R. M., Eds.; The Clay Minerals Society: Boulder, CO, 1992;p 128. (4)Low, P. F. Langmuir 1987,3, 18. (5) Cases, J. M.; BBrend, I.; Besson, G.; Francois, M.; Uriot, J. P.; Thomas, F.; Poirier, J. E. Langmuir 1992,8, 2730. @
a synthetic, macrocrystalline Na-montmorillonite they prepared with only octahedral layer charge. However, a three-layer hydrate of this Na-montmorillonite sample was not observed, possibly because of its very uniform layer charge distribution.6 At 72%relative humidity and above, the swelling process observed by Cases et ale5was isoenthalpic, suggesting to them that adsorbed water in higher hydrates is thermodynamically like water in aqueous solutions, a finding that is in contrast with other studies.l>' Research on molecular behavior in bulk liquid water and aqueous solutions has long benefited from insights provided by Monte Carlo (MC) and molecular dynamics (MD) simulation^.^^^ The underlying philosophy of these simulations is to construct a mathematical description of interactions represented as potential functions and then sample the configuration of a manageable system of molecules in order to ascertain its properties.1° Comparison of the simulation results with experiment not only tests the success with which molecular interactions have been described by model potential functions but also points to new experiments needed to close gaps in understanding. Bleamll has reviewed the few published molecular simu(6) Yamada, H.; Nakazawa, H.; Hashizume, H.; Shimomura, S.; Watanabe, T. Clays Clay Miner. 1994,42,77. (7)Fu, M. H.; Zhang, Z. Z.; Low, P. F. Clays Clay Miner. 1990,38, 485. (8)Beveridge, D.L.; Mezei, M.; Mehroortra, P. K.; Marchese, F. T.; Ravi-Shanker, G.; Vasu, T.; Swaminathan, S. Monte Carlo computer simulation studies of the equilibrium properties and structure of liquid water. In Molecular-Based Study of Fluids; Haile, J. M., Mansoori, G. A., Eds.; American Chemical Society: Washington, DC, 1983;p 297. (9)Bopp, P. Molecular dynamics simulations of aqueous ionic solutions: In The Physics and Chemistry of Aqueous Ionic Solutions; Bellissent-Funel, M.-C., Neilson, G. W., Eds.; D. Reidel: Boston, MA, 1987;p 217. (10)Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987;pp 1, 205. (11)Bleam, W. F. Rev. Geophys. 1993,31,51.
0743-746319512411-2734$09.00/0 0 1995 American Chemical Society
Na-Montmorillonite Double-Layer Structure
Langmuir, Vol. 11, No. 7, 1995 2735
'p'
Si atom in top tetrahedral sheet
Y
isomorphous substitution of Si by A1 atom in top tetrahedral sheet
0
A1 atom in octahedral sheet
0
isomorphous substitution of A1 by Mg atom in octahedral sheet
'p'
Si atom in bottom tetrahedral sheet
Y
isomorphous substitution of Si by A1 atom in bottom tetrahedral sheet
Figure 1. Schematic view of the molecular structure of a Wyoming-type montmorillonite. The unit cell is denoted by a shaded
area. The simulation cell is denoted by a solid gray line enclosing eight unit cells separated by dashed lines.14
lations of smectite interlayer water structure and noted their slow evolution from rather primitive concepts of water-water and clay-water interactions to fairly realistic modeling of both these interactions and the equally important counterion-water and counterion-clay mineral interactions. Skipper et al.12J3and Refson et al.14 have performed three-dimensional MC and MD simulations of high-charge, trioctahedral, or dioctahedral smectites bearing either one or two monolayers of adsorbed water and either Mg2+or Na+ counterions. Their results for a model Na-smectitecontaining only octahedral layer charge and two monolayersof adsorbed water indicated that Na+ counterions formed both inner-sphere and outer-sphere surface complexes. Interlayer water molecules, although they exhibited a much lower self-diffusion coefficient than in bulk liquid water, were not well organized around the Na+ cations, partly because of competing interactions through hydrogen bonding among themselves and with the clay mineral. The same qualitative picture emerged from MC simulations of a model Na-smectite bearing only tetrahedral layer charge and either one or two monolayers of water. These results may reflect the fact2 that Na+ solvation in aqueous solution falls between that of the strongly-solvating Li+ and the innocuous K+. Skipper et al.15J6recently performed molecular MC simulations of monolayer hydrates of dioctahedral Nasmectites containingboth tetrahedral and octahedral layer charge (e.g., Wyoming-type montmorillonite). Their methodological study15 indicated that potential functions for the interlayer species could be based reliably on the MCY model17of water-water interactions and that a simulation domain covering about eight unit cells of a smectite layer could be used successfully with periodic boundary conditions to approximate the behavior of macroscopic clay layers. For the one-layer hydrate, the MC simulations16 revealed the tendency of Na+ counterions to form inner(12)Skipper, N.T.;Refson, K.; McConnell, J. D. C. J. Chem. Phys. 1991,94, 7434.
(13)Skipper, N.T.;Refson, K.; McConnell, J. D. C. In Geochemistry
of Clay-Pore Fluid Interactions; Manning, D. C., Hall, P. L., Hughs, C.
R., Eds.; Chapman and Hall: London, 1993;p 40. (14)Refson, IC;Skipper, N. T.; McConnell, J. D. C. In Geochemistry of Clay-Pore Fluid Interactions; Manning, D. C., Hall, P. L., Hughs, C. R., Eds.; Chapman and Hall: London, 1993;p 1. (15)Skipper, N.T.;Chang, F.-R. C.; Sposito, G. Clays Clay Miner., in press. (16)Skipper, N. T.;Sposito, G.; Chang, F.-R. C. Clays Clay Miner., in press. (17)Matsouka, 0.; Clementi, E.;Yoshimine, M. J.Chem. Phys. 1976, 64, 1351.
sphere surface complexes as the tetrahedral contribution to layer charge increased. Good semiquantitative to quantitative agreement was found between simulation results and experimental diffraction or thermodynamic data. In this paper, we report results of the first concurrent MC and MD simulations of the one-, two-, and three-layer hydrates of Na-Wyomingmontmorillonite,building on the methodologies of Skipper et al.15 and Refson et al.14 The principal objectives of our study were to examine the perturbation of water structure by a model clay mineral surface that is of great importance in swelling clay-water ~ y s t e m s l -and ~ to investigate the partitioning of Na+ counterions into surface complexes and diffuse doublelayer species.
Simulation Methods The general principles of MC and MD simulations for liquids are described in detail in the monograph by Allen and Tildesley.lo We performed MC simulations of hydrated Na-Wyomingmontmorilloniteto obtain estimates of the equilibrium potential energy, layer spacing, interlayer density profiles, and radial distribution functions in an isothermal-isobaric ensemble, with the external pressure applied normal to the clay 1 a ~ e r s . lThen ~ we performed MD simulations of the same system to investigate the self-diffusion coefficient and trajectories of both water molecules and counterions. The model clay mineral studied was a Wyoming-type montmorillonite with the unit-cell formula: Nao.72[Si,.75Alo.25]A23.5Mg0.5020(OH)4. Monte Carlo Simulations. p e simulation cell comprised opposing 21.12 x 18.28 A2 clay mineral patches (about eight unit cells, Figure 1)interacting with six Na+ cations and varying numbers of water molecules. The one-, two-, and three-layer hydrates of Na-montmorillonite were investigated, corresponding to 32, 64, or 96 water molecules in the simulation cell and to the representative water contents of 0.1,0.2, or 0.3 kg of water per kg of clay, r e s p e c t i ~ e l y . ~Simulations ?~J~ of the clay-water systems were conducted at 300 K and lo5 Pa normal stress using the program MONTE, as describedby Skipper et al.15The simulation cell is repeated in all directions via threedimensional periodic boundary conditions. Details of the development and testing of the water-water, waterclay, water-counterion, clay-clay, clay-counterion, and (18) Mooney, R. W.; Keenan, A. G.; Wood, L. A. J. Am. Chem. SOC. 1952, 74, 1371.
Chang et al.
2736 Langmuir, Vol. 11, No. 7, 1995 counterion-counterion potentials used are discussed by Skipper et al.12J3J5 In a simulation of dehydrated Na-montmorillonite, all Na+ were arranged initially on the midplane of a simulation cell with a lo-A layer spacing. Water molecules in the simulations of the hydrated clay mineral were initially arranged randomly and the Na+counterions were asgigned positions over the clay charge sites a t either 5 or 6 A from the centers of both clay layers. The water molecules and Na+ counterions can be selected randomly to move arbitrary distances in any direction; the water molecules also are allowed to rotate. The layer spacing (z dimension of the simulation cell) is dependent on water content and so was allowed to change during simulation by movement of the upper clay layer. For$he monolayer hydrate, the simulation started with a 13-A layer spacing. Because the initial random configuration of the water molecules can deviate significantly from the equilibrium configuration, we allowed only the water molecules to move for the first 10 000 attempted steps. Then another 20 000 steps were attempted in which the upper clay layer also was allowed to move vertically. In this way, the water molecules and upper clay layer could adjust to optimum positions before the Na+ shifted away from their initial positions, thereby avoiding having the counterions flee from the unequilibrated water structure and become entrapped unphysically a t sites directly on the clay mineral surface. Thereafter, all components-upper clay layer, water molecules, and Na+-were allowed to move during 60 000 attempted steps. Equilibrium was conceded when the interlayer potential energy and layer spacing reached plateaus about which they fluctuated normally. After equilibrium was achieved, the simulations were allowed to proceed for a n additional 300 000 attempted steps to collect data on the interlayer atomic density profiles, radial distribution functions, and orientations of interlayer water molecules. For the two- and three-layer hydrates, the number of steps in the first two stages was increased to give each water molecule about the same number of attempted moves before the Na+ were shifted. The initial layer spacings for the two- and three-layer hydrates were 16 and 19 A. The total number of attempted moves required to equilibrate and collect data on dehydrated Na-montmorillonite was 87 500; whereas for the two- and three-layer hydrates 1200 000 moves were required. Interlayer atomic density profiles, radial distribution functions, and orientations of interlayer water molecules were calculated from data collected after every 500 attempted steps for all simulations. Molecular Dynamics Simulations. Equilibrium molecular configurations, with fixed layer spacings as found in the MC simulations of the one-, two-, and threelayer hydrates, were used as the initial configurations in the MD simulations. Simulations using the program MOLDY14involved a n elapsed time of 205 ps for the one-, two-, and three-layer hydrates, using 0.0005 ps per step. Interlayer molecular configurations were stored every 100 steps, then used for mean-square-displacement calculations and trajectory p10ts.l~Self-diffusion coefficients, D , for the motions of either water molecules or Na+ along the mineral surface were estimated from the two-dimensional Einstein relation
D
I d
= - - (7-2)
4 dt
where (1.2) is the (two-dimensional) mean-square displacement and t is the simulation elapsed time. Equation 1is to be preferred over the alternative, integration of the
Table 1. Thermodynamic Properties of Na-Montmorillonite Hydrates Calculated by MC Simulation hydrate AUa (kJ mol-I) d b (A) ec (Mg m+) one-layer 53.5 f 1.8 12.08 f 0.07 1.10 15.28f 0.06 0.91 two-layer 33.5 f 0.9 three-layer 32.3 f 0.7 18.77 f 0.06 0.83 Potential energy per mole of interlayer water and standard deviation. Layer spacing (c-axis) and standard deviation. Calculated interlayer water density, the ratio of water mass per simulation cell to interlayer volume (simulation cell area x interlayer width along c-axis). displacement autocorrelation k c t i o n , whenever modestly long simulation times are involved, as in the present study.1° The maximum mean-square displacement of most water molecules and Na+ exceeded the lateral dimensions of the simulation cell, except in the one-layer hydrate, where all species were very immobile. Interlayer molecular configurations were used to construct average density profiles of Na+ and water molecules to compare with the MC results.
Results Thermodynamic Properties. Monte Carlo simulation results for the differential potential energy of interlayer water, the layer spacing under 100 kPa applied pressure, and the interlayer water density a t 300 Kin the three Na-montmorillonite hydrates are listed in Table I. The differential potential energy of interlayer water was calculated as the difference in total system potential energy between the n and n - 1layer hydrates ( n = 1 , 2 , or 3), divided by the difference in the number of moles of interlayer water between the two. It may be compared to the internal energy of bulk MCY water,8 36.4 f 0.4 k J mol-l. The differential potential energy in the one-layer hydrate is significantly larger than this value and is also in good agreement with measured isosteric heats of adsorption or heats of immersion (after the latter is added to the heat of condensation for liquid water) for monolayer hydrates of Na-Wyoming montmorillonite (48-52 k J mol-').15 By contrast, the differential potential energies for the two- and three-layer hydrates are essentially the same as the internal energy of MCY water, in agreement with the conclusionsof Cases et aL5 The MC layer spacings in Table 1 also are in very good agreement with experimental X-ray diffraction data (12.1 A, n = 1;15.6 A,n = 2; 18.1 A, n = The calculated interlayer water densities decrease from slightly above to slightly below the density of bulk liquid water as the water content increases, a characteristic of the MCY potential function. l9 CounterionDistributions. In the one-layer hydrate of Na-Wyoming montmorillonite, Skipper et al.16 found two species of adsorbed Na+ in their MC simulation: an inner-sphere surface complex near sites of isomorphic substitution in the tetrahedral sheet and an outer-sphere surface complex over sites of isomorphic substitution in the octahedral sheet. In the present simulation of twoand three-layer hydrates, the inner-sphere surface complexes persisted near tetrahedral substitution sites. The Na+associated with octahedral sites, however, partitioned into two species: outer-sphere complexes bound with the ditrigonal cavities of the siloxane surface and highly solvated species residing near the midplane of the interlayer region. These latter species can be regarded as part of an incipient diffuse-ion swarm.3 The partitioning of 3)587316J8
(19) Mulla, D. J. In GeochemicalProcessesat Mineral Surfaces;Davis, J . A., Hayes, K. F., Eds.;American Chemical Society, Washington, DC, 1986; p 20.
Na-Montmorillonite Double-Layer Structure
Langmuir, Vol. 11, No. 7, 1995 2737
l0.W
I0.M 9.00
9.00
~+-+++-+
8.W s.-
1.w
7.M
6.m
6.m
5.W 4.w
5.W 4.m
3.W 2.w n
3
*
3.w
n
1.w
5,
0.w
*
.I.W -2.w
-3.00 4.W
2.00
1.w 0.m .l.W .2.W .3.W
4.w
1
I
I
I
I
0.w
5.w
1o.m
15.00
2o.w
!
I
I
I
0.0
3.m
1o.m
15.w
4.00 .7.m
4.W
-,.-
-9.w
-8.0 -9.m
.IO.W
.10.00
.lO.W
4.m
o.m
3.W
l0.W
'
10.w 9.w
".-
5.w 4.w
3.00 2.m
5 N
.8.W -9.w
1.w 0.W .1.m .2.W
.a.w
.IO.W
.lam
-5.m
0.m
5.m
I
I
I0.W
20.w
Figure 2. Trajectory plots for a Na* counterion near a
Figure 3. Trajectory plots for a Na+ counterion near an octahedralcharge site in the three-layer hydrate: (a)XYplane; (b)XZ plane.
the three adsorbed species was approximately equal in the two multilayer hydrates. These trends are consistent with available spectroscopic and X-ray diffraction data on the speciation of adsorbed Na+ in hydrated Wyomingtype m o n t m ~ r i l l o n i t e . ~ ~ - ~ ~ The contrasting behavior of adsorbed Na+ near the two types of charge substitution site is dramatically illustrated in MD trajectory plots for the three-layer hydrate. Counterions bound near sites of tetrahedral charge substitution hover close to their attracting locus (Figure 2), constrained there to trace out triangular cycles over the three surface oxygen atoms on which the layer charge
deficit is principally manifest.24 Counterions bound near sites of octahedral charge substitution roam much more freely over the siloxane surface (Figure 3),and they stratify their z dimension position either a t the interlayer midplane (two-layer hydrate) or a few angstroms above and below it (three-layer hydrates) as indicated in Figure 4. Interlayer Water Structure. As reported also by Skipper et a1.,16most of the interlayer water molecules in the MC simulation of the monolayer hydrate of Namontmorillonite were found to reside a t the midplane (Figure 5). In the two-layer hydrate, water molecules were distributed roughly equally around two planes offset from the midplane, and in the three-layer hydrate, besides the two planes of water molecules now offset even further from the midplane, a second group of water molecules was clustered loosely near the midplane into two minor
tetrahedral charge . site in the three-layer hydrate: (a)XYplane; (b)XZ plane.
(20)Pezerat, H.; MBring, J. C.R.Acad. Sci., Ser. D 1967,265, 529. (21)Sposito, G.; Prost, R.; Gaultier, J.-P. Clays Clay Miner. 1983,31, 9.
(22)Ben Brahim, J.;Besson, G.; Tchoubar, C. J.Appl. Crystallogr. 1984,17, 179. (23) Ben Brahim, J.; Annagan, N.; Besson, G.; Tchoubar, C. Clay Miner. 1986,21, 111.
(24)Sposito, G. The Surface Chemistry of Soils;Oxford University Press: New York, 1984;p 14.
Chang et al.
2738 Langmuir, Vol. 11, No. 7, 1995 4.00
1
1
r -a-
I
one-layer
...D.- two-layer three-layer
3.00
b
. I
B
2.00
n 1S O
I Y I
3.W
ILj
2.m
I I I t
1.w
\
I
I
\ I 1
I
\
V I
1.00
\
y/
0.50
0.00 -5.00
0.00
5.00
z (A> Figure 6. Interlayer density profiles for water molecules (strictly, 0 atoms) in Na-montmorillonite hydrates, with the midplane taken as the origin in the z direction.
12.M
;
I1.W
-
:::: 1o.m
I
,
I
I
I
I
I
om
Lm
I
I
.
I
9w 100 1w
6W 5W
rw 3m 2W
IW OW 4 W
.l
w
400
U h . I
m
8
n
0.00
5.00
z Figure 4. Interlayer density profiles for Na+ counterions in
(a)one-layerhydrate, (b)two-layerhydrate, and (c) three-layer hydrate. Results are from MD simulations, plotted relative to the midplane as origin. layers (Figure 5). The solvation shell water of surfacecomplexed Na+ usually contributed to the major layers
that were offset from the midplane, while water molecules associated with the diffuse-ion swarm were found in the midplane region. The small peaks in the wings of the distributions in Figure 5 represent water molecules located in ditrigonal cavities of the siloxane surface (not necessarily associated with outer-sphere surface complexes or octahedral charge sites), as detected spectroscopically by P r o ~ t .The ~ ~ MD trajectory plots suggested that the solvation-shell, ditrigonal-cavity, and network water molecules all undergo exchange in the two- and threelayer hydrates on the time scale of the MD simulations. Water molecules can enter and leave ditrigonal cavities (as observed also in MD simulations by Refson et al.14for the two-layer hydrate of Mg-smectite), or can associate with and dissociate from counterion Na+from time to time, in agreement with neutron scattering data.l Orientations of the interlayer water molecules in the MC simulations are summarized in Figure 6 . In the onelayer hydrate, most water molecules orient with their dipole moment vectors parallel to the siloxane surface (position 4 in Figure 6), but some may orient with their OH bonds almost parallel to the siloxane surface (close to position 2). In the two-layer hydrate, most water molecules rotate about 38@clockwise from position 4 so that their OH bonds orient perpendicularly to the siloxane surface (position 3). In the three-layer hydrate, water molecules on the plane nearest to the siloxane surface have an orientation similar to that in the two-layer hydrate (positions 3 and 51, whereas water molecules near the midplane can orient among different positions (Figure 6b). The 0-0 and H-0 radial distribution function (RDFs) in the Na-montmorillonite hydrates and in bulk water16 are compare4 in Figures 7 and 8. The first 0-0 peak is a t 2.7 o r 2.8 A for bulk water or the Na-montmorillonite hydrates, respectively. The second 0-0 peak in the onelayer hydrate, near 4.5 A, is much better defined than in bulk water.16 The H-0 RDFs indicate that hydrogen bonds in the one- and two-layqr hydratgs are longer than that in bulk water (i.e., 2.0 A vs 1.8 A). These trends (25) Prost, R. In Proceedings of the International Clay Conference, MBxico City, 1975;Bailey, S. W., Ed.; Applied Publishing Ltd.: Wilmette, IL, 1976; p 351.
Langmuir, Vol. 11, No. 7, 1995 2739
Na-Montmorillonite Double-Layer Structure 3.00
-a-
2.80
bulk water
2.60 I
I
I
I
0
0
w
90
38 52
52 38
90 0
2.40 2.20 2.00 1.80
I .60 I .40 1.20
0
w
I
I
I
128
142
180
I .00
90
0.80
38
52
0.60
0.40
10-3
0.20 0.00 1.50
2.00
2.50
3.00
3.50
4.00
4.50
5.00
Distance from oxygen atom (A)
Figure 7. Graphs of the 0-0 RDF for water molecules in the
bulk liquid16and interlayer water. 2.00
-m-
1.90
bulk water
1.80
I .70
20.00
-8
Ii
".- I 0.00
1.60 1.50
. I
I
I
I
50.00
100.00
lSO.00
1.40
3
1.30
g
1.20
.n
).1
a
0
s 3
1.10 1.00
4
0.90
I
;Ei
11~3
m.00 540.00
520.00
0.70
E
0.60 0.50 0.40
500.00
s
0.80
E
0.30
480.00 460.00 440.00
0.20
L
420.00
p
=
6 h c1 .Y
'3
'')
=
380.00 360.00 340.00
1.50
2.00
2.50
3.00
3.50
4.00
4.50
5.00
Distance from hydrogen atom (A)
320.00
300.00
Figure 8. Graphs of the H-0 RDF for water 0 atoms in the
280.00
bulk liquid16and interlayer water.
260.00
240.00 200.00
180.00 160.00 140.00 20.00
40.00
60.00
80.00
w Figure 6. Orientation of interlayer water molecules relative
to a normal on the nearest siloxane surface (arrow) in Namontmorillonite hydrates: (a) sample orientations of dipole axis and H-H axis relative to the arrow; (b)plots of the dipole axis distribution P( in Na-montmorillonite hydrates, where 8 is the angle relative to the arrow in (a); (c) plots of the H-H axis distribution P(V)16in Na-montmorillonite hydrates, where V is the angle the H-H axis makes with the arrow in (a)when both are in the same place. imply that the interlayer water is more structured than bulk water, although the differences between interlayer
and bulk water are less pronounced in two- and threelayer hydrates. The Na-0 RDFs in the Na-montmorillonite hydrates and for Na+ in bulk wat?rl6 are shown in Figure 9. The first peak occurs near 2.3 Ain the hydrates, which is the same as for Na+ in bulk water. However, the second Na-0 peak in bulk water is better defined than that in the Na-montmorillonite hydrates. Interlayer Water and Counterion Mobility. Selfdiffusion coefficients for water molecules and Na+ counterions, estimated from linear regression of the calculated mean-square displacements on elapsed time in the MD simulations, are listed in Table 2 [denotedby (MD)]along with relevant experimental data.26t28J9The MD simula(26)Cebula, D.J.;Thomas, R. K.;White, J. W . Clays Clay Miner. 1981,29,241. (27)Sposito, G. J. Chem. Phys. 1981,74, 6943.
Chang et al.
2740 Langmuir, Vol. 11, No. 7, 1995
calculations and experimental self-diffusion coefficients is quite good for the three-layer hydrate. The qualitative behavior of both the MD and experimental self-diffusion coefficients as functions of water content are in accord: a very rapid rise in magnitude in going from one to two water layers, with self-diffusion being significantly more restricted in the one-layer hydrate than in bulk water (last row in Table 2). The MD simulations indicated also that, in the two- and three-layer hydrates, Na+ in outersphere surface complexes diffused much more rapidly than Na+ in inner-sphere complexes.
8.00 7.50
7.00 6.50 6.00
5.50 5.00 4.50 4.00 3.50
3.00 2.50
2.00 1.50
I .00 0.50
0.00 l,50
2.00
2.50
3.W
3.50
4.00
4.50
5.00
Distance from sodium atom (A)
Figure 9. Graphs of t h e Na-0 RDF for Na+ in bulk solution16 and interlayer water. Table 2. Self-DiffusionCoefficients of Interlayer Water Molecules and Sodium Counterions at 300 K
one-layer two-layer three-layer bulk water
2.5 x 1.4 x
4x 1x
1.3 x
1x 1.5
10-9
3.9 x
1x
5.1 x 1.8 x
1x 2x 9 x 10-10
D, for one- and three-layer Li-montmorillonite hydrates;26D, for the two-layer hydrate of N a - m o n t m o r i l l ~ n i t e ;bulk ~ ~ water estimate is 2/fl, a t 298 RZ7 D Nfor ~ Na-Wyomingmontmorillonite hydrates, and 2/3 of the value of the infinite-dilution-limit D N in ~ aqueous solution at 298 K.2s a
tions involved a total elapsed time of 205 ps for the one-, two-, or three-layer hydrates, which is long compared either to the residence time of a water molecule in the bulk liquid (1ps)27or to that of a water molecule solvating Na+ in solution (10 P S ) . ~ O It is still possible, however, that a n even longer total elapsed time would permit more accurate estimates of D, and D N ~ . The MD estimate of D, in the one-layer hydrate is of the same order as D, for water molecules in the one-layer hydrate of Li-montmorillonite, as obtained from a jump/ rotational-diffusion model analysis of quasielastic neutron scattering data by Cebula et a1.26 The differencethat exists in the two D, may reflect either inadequate neutron scattering data analysis or insufficient MD sampling time, since a water molecule residence time of 43 ps was inferred from the neutron scattering data Hall et al.29 analyzed quasielastic neutron scatteringdata for the twolayer hydrate of Na-montmorillonite using a model that assumed no diffusive motion (on the neutron-scattering time scale) of water molecules solvating Na+ counterions and prescribed translational jump diffusion for all other water molecules-thus a rather different model analysis from that of Cebula et a1.26-to deduce D, = 9.4 x 10-lo m2 s-l, in good agreement with the MD estimate in Table 2 (second row of column 1). Agreement between the MD (28) Nye, P. H.Adu. Agron. 1979,31, 225. (29) Hall,P.L.;Ross,D. K.;Tuck, J . J.;Hayes,M. H.B.InProceedings of the ZAEA Symposium “Neutron Inelastic Scattering“, Vienna, 1977; International Atomic Energy Agency: Vienna, 1978; Vol. 1, p 617. (30)Enderby, J. E.; Nielson, G. W. Rep. Prog. Phys. 1981,44, 593.
Discussion The Na+ density profiles in Figure 4 show that innersphere surface complexes coexist with outer-sphere surface complexes in all three hydrates. Some of the outer-sphere surface complexes dissociate from the siloxane surface to become diffuse-layer species in the multilayer hydrates. Thus, the partitioning of adsorbed Na+ involves all three possible surface species,24with the fraction of inner-sphere surface complexes determined by the fraction of the layer charge arising from isomorphic substitutions in the tetrahedral sheet. In natural Na-montmorillonite, layer charge is not uniformly distributed among the clay particles, as assumed in the present study. It is possible that some clay particles may bear solely inner-sphere surface complexes, others a mixture of adsorbed species, and so on, with only the average partitioning approximating the results in Figure 4. In the multilayer Na-montmorillonite hydrates, the solvated counterions and water molecules are relatively mobile (Figure 3) and can undergo exchange among peak positions in the z dimension. However, the mobilities of Na+ and the water molecules do not always increase with increasing water content for these two hydrates (Table 2). If the smaller Na+ mobility in the three-layer hydrate is not a computational artifact, it may reflect a greater electrostatic gradient which makes it more difficult for Na+ cations to depart from the siloxane surface. This hypothesis could also explain why a more populated diffuse-ion swarm was observed in the two-layer hydrate (Figure 4). In the one-layer Na-montmorillonite hydrate, the water molecules reside mainly near the midplane (Figure 5) because there is only one water network, whereas in the two-layer hydrate, their number density is higher near the siloxane surface than near the midplane becaus? there are effectively two water networks a t about 6.0 A from the center of the octahedral sheets. In the three-layer hydrate, the midplane water is less well-structured than the near-surface water networks (Figure 5). This may reflect a greater rotational freedom of the midplane water molecules, as noted also in a Fourier analysis of neutron diffraction data for Na-Wyoming m~ntmorillonite~l over the same range of water contents a s in the present study (cf. Figure 6). Several issues arose in the course of the simulations that should be kept in mind when interpreting the results of the present study. For example, we found that similar interlayer potential energy levels among differing cation configurations are separatedby “energy barriers”, making it difficult to determine the cation configurations precisely within a n affordable simulation time. These “energy barriers” can prevent Na+ cations from passing through the water network and thereby make sampling all configurations impossible within a finite time. This problem of a system becoming trapped in a too-smallregion of phase space means that MC simulations started with (31) Hawkins, R.K.;Egelstaff, P. A. CZays Clay Miner. 1980,28,9.
Nu-Montmorillonite Double-Layer Structure different initial configurations andlor model parameters can result spuriously in different “equilibrium”interlayer molecular structures. Therefore, we cannot warrant that an MC simulation started from an arbitrary initial configuration would certainly lead to a unique equilibrium interlayer molecular configuration. However, we did find that simulations started from different initial configurations were consistent in their thermodynamic properties and in the positions of Na+ in inner-sphere complexes. The mean interlayer potential energy varied by less than 3%or 4% in the two- or three-layer hydrate, and the layer spacing varied by less than 5%in the two-layer hydrate. In the three-layer hydrate, however, drift of the layer spacing was found if the simulations were started with Na+ placed solely near octahedral surface charge sites. The increase in layer spacing during the course of a MC simulation was usually accompanied by a partitioning of the midplane water molecules into two networks. The MCY potential function used in the present study is known to have the problem of underestimating the mass density of bulk water.15 This should not cause difficulty in oneand two-layer hydrates, where screened Coulomb forces hold the system together. In the three-layer hydrate, this force may no longer suffice if the layer charge is substantially neutralized by Na+cations very near the mineral surface. Thus, it may be necessary to place Na+ initially over octahedral charge sites in the middle of a simulation cell solely to stabilize the interlayer water structure, as well as the layer spacing, in the three-layer hydrate. These and other sampling-related issues are topics for continued simulation research.
Conclusions Monte Carlo and molecular dynamics simulations of interlayer molecular structure, performed using the
Langmuir, Vol. 11, No. 7, 1995 2741 methodologies of Skipper et al.15and Refson et al.,14were in agreement with available experimental data on the one- to three-layer hydrates of Na-Wyoming montmorillonite, the prototypical swelling clay mineral. Good predictions were obtained for layer spacings, interlayer water potential energies, and the significant deviation of the water structure from that in the bulk liquid. Expansion of the interlayers was accompanied by shifts in the partitioning of adsorbedNa+ species between outer-sphere surface complexes and the diffuse-ion swarm, but with persistence of inner-sphere surface complexes near sites of tetrahedral charge substitution. The interlayer water structure became increasingly less organized, primarily because of growing rotational and translational freedom of the molecules, but did not replicate that in bulk water or aqueous solution even in the three-layer hydrate. More accurate simulations should be accompanied by more precise diffraction, neutron-scattering, and spectroscopic experiments on molecular structure in the interlayer hydrates of Na-montmorillonite.
Acknowledgment. The research reported in this paper was supported in part by NSF Grant EAR-9206052. The authors thank the San Diego Supercomputer Center and the Pittsburgh Supercomputer Center for allocations of CPU time on Cray Y-MP8/864 and Y-MP/CSO supercomputers. The first author expresses gratitude to Bruce Berne for his hospitality during her visit at Columbia University. Gratitude also is expressed to Keith Refson, University of Oxford, for assistance with the MD simulations. Thanks to Joan Van Horn for excellent typing of the manuscript. LA940869T
