Langmuir 1991, 7, 547-555
547
Modeling the Clay-Water Interface Alfred Delvillet Laboratoire de chimie fine, biomim6tique et aux interfaces, Ecole Polytechnique, 91128 Palaiseau Cgdex, France Received February 22, 1990. I n Final Form: May 9, 1990 To model the clay-water interactions, we use a molecular description of the solvent and clay sheet. The clay-water and clay-cation potentials are obtained from quantum calculations. These potentials are used in Monte Carlo simulations of the solvation of two montmorillonite sheets and their interlamellar sodium counterions. We obtain an immersion enthalpy and mean water orientation in agreement with the experimental data. The results show the importance of the solvation of the interlamellar cations on the swelling properties of the clay.
I. Introduction Clays are lamellar aluminosilicates showing a large variety of physicochemical properties: swelling, adsorption, ionic exchange, surface acidity, etc. Their use are widespread: gelling, decoloring,l water pollutant elimination,2radioactive waste r e p ~ s i t o r y~, ~r a c k i n gand , ~ heterogeneous c a t a l y ~ i s . ~Experimental -~ investigations of the clay-substrate interactions are countless, performed by various techniques, such as thermal analysis8 or calorimetry,*" mechanical measurements,12-14or spectral measurements, such as neutron diffra~tion,'~ X-ray diffra~tion,'~J' IR,18J9 UV,20EPR,21 and NMR11p22-25spectroscopies. t Present address: CRMD-CNRS, 1B rue de la Ferollerie, 45071 Orleans Cedex 02,France. (1) Herjavec, S.; Maric, J.; Mlinaric, Z.; Krauthaker, V.; Razumovic, D. Jugosl. Vinograd. Vinar. 1987,21,19-22. (2)Ray, C.; Ramsey, R. Enuiron. Prog. 1987,6,145-154. (3) Baetsle, L. H.; Henrion, P.; Put, M. High-Level Nuclear Waste Disposal. Proceedings of thelnternational Topical Meeting; Burkholder, H. C., Ed.; Columbus OH, 1985;pp 313-323. (4)Ryland, L. B.; Tamele, M. W.; Wilson, J. N. Catalysis; Emmett, P. H., Ed.; Rheinhold New York, 1960;Vol. 7, Chapter 1. (5)Pinnavia, T. J. Science 1983,220,365-371. (6)Laszlo, P. Science 1987,235, 1473-1477. (7)Laszlo, P. Acc. Chem. Res. 1986, 19, 121-127. (8)Mackenzie,R. C. Aduanced Techniquesfor Clay Mineral Analysis; Fripiat, J. J., Ed.; Elsevier Science Publishers: Amsterdam, 1982;Chapter 2. (9)Mikhail, R. S.; Guindy, N. M.; Hanafi, S.J . Colloid Interface Sci. 1979,70,282-292. (10)Oliphant, J. L.; Low, P. F. J . Colloid Interface Sci. 1982,89,366nrn J 1J.
(11)Fripiat, J.; Cases, J.;FranFois, M.; Letellier, M. J . Colloid Interface Sci. 1982,89,378-400. (12)Callaghan, I. C.; Ottewill, R. H. Faraday Discuss. Chem. SOC. 1974,57, 110-118. (13)Israelachvili, J. N.: Adams, G. E. J . Chem. Soc., Faraday Trans. 1 1978,74,975-1001. (14)Israelachvili. J. N.:Pashlev. R. M. Nature 1983. 306. 249-250. (15)Cebula, D. J.'; Thomas, R. K.fWhite, J. W. J . Chem. Soc. Faraday Trans. 1 1980, 76, 314-321. (16)Pons, C. H.; Rousseaux, F.; Tchoubar, D. Clay Miner. 1981,16, 23-42. (17)Viani, B. E.; Low, P. F.; Roth, C. B. J . ColloidInterface Sci. 1983, 96,229-244. (18)Suquet, H.; Prost, R.; Pezerat, H. Clay Miner. 1977,12,113-126. (19)Sposito, G.; Prost, R. Chem. Reu. 1982,82,553-573. (20)Villemure, G.;Detellier, C.; Szabo, A. G. J . Am. Chem. SOC.1986, 108,4658-4659. (21)Boyd, S.A.; Mortland, M. M. Nature 1985,316,532-535. (22)Woessner, D. E.; Snowden, B. S. J. Chem. Phys. 1969,50,15161523. (23)Grandjean, J.;Laszlo, P. Clays Clay Miner. 1989, 37,403-408. (24)Bottero, J. Y.; Bruant, M.; Cases, J. M.; Canet, D.; Fiessinger, F. J . Colloid Interface Sci. 1988, 124,515-527. (25)Delville, A.; Grandjean, J.; Laszlo, P. J. Phys. Chem., in press.
0743-7463f 91f 2407-0547$02.50/0
Many authors have proposed models to reproduce the interactions of the solvent-ion-support However they neglect either the structure of the solvent, which is treated as a continuum, or the structure of the support, which is treated as a smooth surface. A complete and satisfactory modeling of all the observed data requires a molecular description of the interacting species. The capabilities of computers cannot afford a complete solution to this problem. This article is a first attempt in this direction. It results from a compromise between the complexity of the model, the statistical accuracy of the calculations, and the computing time. Here we limit our study to the clay-cation-water system. But there is no intrinsic restriction; this study can be potentially extended to systems including organic matter, for example. This article is divided in three parts. First, we perform quantum calculations of the clay-water and clay-cation interactions. Second, we parametrize these energies to obtain a simple analytical expression, easy to use for statistical simulations or molecular mechanic optimizations. Third, we use Monte Carlo simulations to treat the solvation of a montmorillonite clay. We calculate an enthalpy of immersiongJ and a mean orientation of the adsorbed water,23125 in agreement with experimental data. This treatment allows us to identify and quantify the driving forces responsible for clay swelling.
11. Results A. Quantum Approach. To calculate the clay-water and clay-cation energies, we use a quantum approach in the framework of the supermolecule. The calculations include only one water molecule or one counterion in the vicinity of a model of montmorillonite. To reproduce the (26)Sullivan, D. E.; Stell, G. J . Chem. Phys. 1978,69,5450-5457. (27)Snook, I.; van Megen, W. J. Chem. Phys. 1979, 70, 3099-3105. (28)Sullivan, D. E.;Levesque, D.; Weis, J. J. J . Chem. Phys. 1980,72, 1170-1174. (29)Eggebrecht, J. M.; Isbister, D. J.; Rasaiah, J. C. J . Chem. Phys. 1980, 73,3980-3986. (30)Lee, S. H.; Rasaiah, J. C.; Hubbard, J. B. J . Chem. Phys. 1986, 85,5232-5237. (31)Plischke, M.; Henderson, D. J . Chem. Phys. 1986,84,2846-2852. (32)Kjellander, R.;Marcelja, S. Chem. Phys. Lett. 1987,142,485-491. (33)Torrie, G.M.; Kusalik, P. G.; Patey, G. N. J . Chem. Phys. 1988, 89,3285-3294. (34)Schoen, M.; Cushman, J. H.; Diestler, D. J.; Rhykerd, C. L. J . Chem. Phys. 1988,88,1394-1405. (35)Attard. P.:Kiellander, R.: Mitchell, D. J.; Jonsson, B. J . Chem. Phys. 1988,89,1664-1680. (36)Delville, A.; Laszlo, P. Nouu. J . Chim. 1989,13, 481-491. Woodward, C.; Jonsson, B. J. Chem. Phys. 1989,91, (37)Akesson, T.; 2461-2469. (38)Aloisi, G.; Foresti, M. L.; Guidelli, R.; Barnes, P. J.Chem. Phys. 1989,91,5592-5596.
0 1991 American Chemical Society
548 Langmuir, Vol. 7, No. 3, 1991
Delville Table I. Energy of the Clay Fragments
8
A12Si&H18
I, eV EBV,eV Eev - EHO,eV AHcd, kJ/mol AHopt,kJ/mol
8.815 1.107 9.923 -13.3 -11.3
A13Si5024H18 Al4Si~,030H24 A15Si5030H24 3.342 4.073 7.415 -27.3 -26.7
8.825 1.360 10.126 -15.6 -15.6
3.638 3.560 7.198 -28.9 -28.9
Table 11. MIND0 Local Charge (au) of the Clay Atoms in the Unit Cell
Figure 1. Model of clay used in the MO calculations.
E
A12Si602,H,8
AI,Si,O,,H,,
AI,Si,O,,H,,
A15Si,0,,l12,
eV
===
Si Si Si Si
0 0 0 0 0 0 0 0 0 0 0 0 A1 AI A1 A1
Figure 2. Energy levels near the Fermi level of different clay
fragments.
symmetry of the cavity of the silicates, the clay fragment contains a hexagonal cavity of silicates plus two inner octahedral aluminia (see Figure 1). We ignore the second layer of silicates. This simplification modifies the energy of the clay fragment but has no effect on the energy of adsorption, since the contribution of the contact layer overcomes the contribution of this remote layer. We take into account the tetrahedral substitution of one SiTvby one AlIII, introducing a negative charge in the clay network. Some tested calculations include four octahedral aluminum. All dangling bonds are replaced by hydrogen atoms, which are not shown in Figure 1. The general formulas of the smaller clay fragments are AlzSis024H18 and A13Si&H18-. The formulas of the fragments with four octahedral aluminum are A14Si6030Hz4 and A15Si5030Hz4-. The coordinates of the clay fragments were obtained by crystallographic measurement^.^^ The calculation uses a MIND040 approximation, included in the MOPAC program of the QCPE. The orientation of the hydrogen in the center of the cavity is obtained by a minimization of the energy of the clay fragment. The OH bond is directed parallel to the clay sheet with the hydrogen pointing to the octahedral vacancy (see Figure l ) ,as expected for dioctahedral clays.lg Figure 2 shows the closest empty and occupied energy levels from the Fermi level. Adding two octahedral aluminums does not modify the position of the LUMO and HOMO orbitals (see Figure 2) but increases the splitting of the energy levels. Tetrahedral substitution greatly modifies the position of the energy levels, reducing the gap between the LUMO and HOMO (39) Maegdefrau, V. E.; Hofman, U. 2.Kristallgr., Kristallgeom., Kristallphys. Kristallchem. 1937, 9, 299-323. (40) Dewar, M. J. S.; Thiel, W. J.Am. Chem. SOC.1977,99,4899-4907.
0 0 0 0 0 0 H H H H H H H H
X
Y
0.0 0.0 2.59 2.59 1.295 3.885 3.885 1.295 0.0 2.59 0.0 2.59 2.59 0.0 5.29 0.0 1.727 1.727 4.317 4.317 0.863 3.453 0.863 3.453 0.863 3.453 2.127 -0.463 0.863 3.453 0.863 3.453 0.863 3.453
2.99 5.98 1.495 7.475 2.2425 2.2425 6.7275 6.7275 4.485 8.97 2.99 1.495 7.475 5.98 4.485 8.97 8.97 2.99 4.485 7.475 1.495 2.99 4.485 5.98 7.475 8.97 5.285 9.77 1.495 2.99 4.485 5.98 7.475 8.97
Z 0.0 0.0 0.0 0.0 0.59 0.59 0.59 0.59 0.59 0.59 -1.59 -1.59 -1.59 -1.59 -1.59 -1.59 -2.68 -2.68 -2.68 -2.68 -3.77 -3.77 -3.77 -3.77 -3.77 -3.77 -1.68 -1.68 -4.7 -4.7 -4.7 -4.7 -4.7 -4.7
4
1.7 1.7 1.7 1.7 -0.93 -0.93 -0.93 -0.93 -0.93 -0.642 -0.93 -0.93 -0.93 -0.93 -0.642 -0.642 1.1 1.1 1.1 1.1 -0.29 -0.642 -0.29 -0.642 -0.29 -0.29 0.21 0.14 0.14 0.21 0.21 0.21 0.21 0.21
from 10 to 7.3 eV (see Table I). This clearly shows the influence of the substitutions on the electrical properties of the metallic oxides. Table I1 shows the mean local charges (in atomic units) for the atoms of the unit cell of the montmorillonite. The A1 in a tetrahedral site bears a local charge of 0.82 au, while its two nearest neighboring octahedral A1 atoms bear a charge of 1.04 au instead of 1.10 au. We calculate the clay-water energy for different configurations of the water molecule (251 configurations), approaching different sites of the clay surface. The geometry of the water molecule and of the clay fragment is fixed. The water OH bond length is assumed to be 0.9572 A with a HOH angle equal to 104.52°.41 The lower energy is obtained when the water molecule points the two hydrogens to the clay surface (see Figure 3), in opposition to the repulsive energy when the oxygen is directed to the clay surface. A neutral fragment stabilizes the water molecule by -13 to -16 kJ/mol (see Table I). This stabilization energy is lower than the water-water stabilization in a bulk aqueous phase, as estimated by the vaporization energy (-40.8 kJ/mol). The tetrahedral substitution increases the water stabilization by a factor of 2; it becomes equal to -29 kJ/mol. Figure 4 shows the energy of the (Na+-Al3Si@24Hl8-) (41) Matsuoka, 0.;Clementi, E.; Yoshimine, M. J . Chem. Phys. 1976,
64, 1351-1361.
Modeling the Clay-Water Interface
Table 111. Lennard-Jones Parameters for the Clay-Water Interaction kJ A2/mol bio, kJ A6/mol b i H , kJ A6/mol cio, kJ A12/mol
ajo,
Si ALt
3.886 X 2.396 X 1.322 X -2.775 X -2.023 X
&et
0 H
Langmuir, Vol. 7, No. 3, 1991 549
loo
lo2 lo5 lo5 lo3
3.675 X -3.360 X 2.971 X 3.531 X -7.311 X
lo2 lo2 loo 10'
-8.200 X 2.185 X -9.474 x -1.464 X -2.238 X
10'
10'
lo5 103
lo3 10'
7.352 X 3.045 X 6.871 X 4.416 X 3.006 X
lo6 1O'O lo7 lo6 lo6
CIH,
kJ A12/mol
-2.390 X -9.389 X -4.624 X 2.512 X -8.584 X
lo4 lo8 lo6 lo4 lo3
B. Parametrization. 1. Clay-Water Interaction. The calculated energies of the (A12Si6024H18-H20) and the (A13Si5024H18-H20)- systems are fitted by LennardJones parameters and electrostatic contributions
E = E,,
+ E,,
(la)
with I
.....
\
'.
\
'
.
,
....."'
' .
.... c
I
and
/ '
~~
5
10
Figure 3. Clay-water MINDO energy (kJ/mol), with (- and - -) or without (- - - and tetrahedral substitution. The energy is attractive when the protons are directed to the clay and - -1 and repulsive when the oxygen is directed to the clay (- and - - -). The distance from the Si plane is measured in
-.
e)
(e.
.
A.
I
It reduces the number of parameters and ensures the
AH 1
\ \
.\
r
r
*
_ _ _ _ _ _ _ - - - - - - --
-*0°
-
t
,
-,-.-
/
I '\.-* ' *' .LOO
The i and j indices stand for the atoms of the support and the water molecule, respectively. Different sets of parameters (bio, ci0 and b i ~C, ~ H )are used for each type of atom of the clay (Si, octahedral Al, tetrahedral Al, 0, H). A supplementary term in l / r 2 is added to the classical Lennard-Jones potential in order to improve the fit. However, a restriction is used
/
/
/
/
/
'L d
L
20
10
Figure 4. Clay-counterion MINDO energy (kJ/mol) for the (Na-clay) (- - -) and (Ca-clay)+ (- -) systems. The distance from the Si plane is measured in A.
-
Table IV. Extrema of the Clay-Water Energies (kJ/mol) as a Function of the Distance (A) of a Water Molecule from the Si Plane no tetrahedral A1
one tetrahedral A1
D,A
Emin
Emax
Emin
Emax
10 5 4 3.7 3.5 3
-22.6 -42.6 -47.2 -46.4 -44.2 -23.5
-20.6 -38.6 -37.3 -28.8 -14.9 151
-25.8 -55.6 -63.2 -62.0 -59.6 -41.2
-21.9 -42.2 -40.9 -32.4 -18.4 148
and (Ca2+-A13Si5024H18-) systems. The interaction is strong, around 400 kJ/mol for sodium and a factor of 2 larger for calcium. It is mainly electrostatic, except at small distances, where the effect of the excluded volume increases.
convergence of the potential at large distance. The electrostatic energy is calculated from the MINDO local charges for each atom of the clay fragment and the water molecule. The electrostatic potential is multiplied by a factor a. The MINDO local charges are 0.159 au for the hydrogen and -0.318 au for the oxygen of the water molecule. These 26 parameters are adjusted by a simplex42 optimization procedure, fitting simultaneously all the data. The mean relative deviation is 5 9; for the neutral fragment and 3% for the substituted fragment. The Table I11 shows the Lennard-Jones parameters resulting from the fitting. The parameter a is adjusted to 1.741. The simulation of the solvation of the clay surface requires a larger model of clay. The model of clay used in the next section contains 24 unit cells and extends over 31.08 A in the X direction and over 35.88 8, in the Y direction. The charge density of the montmorillonite (6.277 X e l e ~ t r o n / A spreads ~ ) ~ ~ on all the octahedral aluminum, after a possible deduction of the tetrahedral charge. The Table IV shows the extrema of the clay-water interaction with and without tetrahedral substitution. The assumed water configuration is the more stable with the two hydrogens pointing to the clay surface. The water molecule remains parallel to the (OY,OZ) plane. The distances are given in A, between the atom of oxygen of the water molecule and the Si plane of the clay. Far from the clay surface, the energy is nearly constant and does not vary with tetrahedral substitution. For neutral clay, the water stabilization (-47 kJ/mol) has the same order of magnitude as the water-water interaction (42) Deming, S. N.; Morgan, S. L. Anal. Chem. 1973, 45A, 278-283. (43) LOW,P. F. Soil S C ~SOC. . Am. J. 1980, 44,667-676.
Delville
550 Langmuir, Vol. 7, No. 3, 1991 I
\
/
/
4.
. !
I.
I
I
I
I
* \
1
\
\
I .
I
I
I /
\
I
I
\
.
I I
0
\
'L
I
X
I X
/
\
Figure 5. Map of the clay-water potential. The water molecule is located 5 8( from the clay. The symbol ( 0 )represents a Si atom and the symbol (0) represents a tetraheral A1 atom: -55 kJ/mol; - - -, -50 kJ/mol; -, -45 kJ/mol.
-.
e,
Figure 8. Map of the clay-water potential 3 A from the clay: -.-,-41 kJ/mol; - - -, -20 kJ/mol; - -, 0 kJ/mol; -, 50 kJ/ mol; -, 100 kJ/mol.
-
-.
Table V. Lennard-Jones Parameters for the Clay-Sodium Interaction Si
I ' I
I 1
Alwt
ALt
'-
0 H
\ \
\
\ \*
-
I'
Figure 6. Map of the clay-water potential 4 8( from the clay: *, -63 kJ/mol; - - -, -60 kJ/mol; - -, -55 kJ/mol; - -, -50 kJ/mol; -, -45 kJ/mol.
--
b, kJ A6/mol
d, k J Ag/mol
-9.780 -1.632 -1.693 -1.202 -2.875
7.016 X lo2 -2.995 X lo4 -4.136 X lo6 -1.710 X lo4 -8.897 X lo5
lo3 lo6 lo5 x 104 x 104 X X X
c,
kJ A1*/mol
5.911 X 1.264 X 1.123 X 8.858 X 2.669 X
lo5 lo6 lo7 lo4 lo8
a saddle point for the translation from one cavity to the other. A t short distance, the oxygen atoms become strongly repulsive, due to the excluded volume effect. 2. Clay-Cation Interaction. Again, the energy of the (Na+-A1&024Hls-) system is fitted by an electrostatic and Lennard-Jones potential. A supplementary additive parameter C is necessary. The Coulomb interaction is of infinite range. It is impossible to obtain a configuration of reference
E = C + (YE,,+ E ,
(3)
In eq 3 the electrostatic potential is calculated from the local MIND0 charges
'L X
and the Lennard-Jones potential is given by 50
Figure 7. Map of the clay-water potential 3.5 A from the clay: .,-59 kJ/mol; - - -, -50 kJ/mol; - -, -40 kJ/mol; -, -30
..
- e .
kJ/mol.
(-41 kJ/mol). For a substituted clay, the water stabilization increases by -16 kJ/mol, while the total charge density of the clay is the same. In both cases, the larger stabilization is reached for a clay-water distance of 4 A. Close to the surface, the energy exhibits some large variations. The transverse diffusion is not free; the activation energy for the translation may reach an upper limit of 180-190 kJ/mol. Figures 5-8 show the isopotential curves of one water molecule as a function of the distance from the surface of the clay with a tetrahedralsubstitution. The configuration of the water molecule is the same as in Table V. A t large distance, the minimum occurs at the site of substitution. As the water molecule approaches, the minimum shifts gradually to the center of the cavity above the site of substitution. Secondary minima appear at the center of the other cavities. The site of substitution behaves like
Ew
bi
di
6
9
ci
=E--+-+ri, rii rii
12
(5)
i=l
The simplex procedure gives a mean relative deviation of 3% with the parameters shown in Table V. The coefficients a and C take the values 1.08and -7.3286 MJ/ mol, respectively. C. Monte Carlo Simulations. As stated previously, the model of clay used in the simulations contains 24 unit cells. Two clay sheets are facing each other, with a distance of 20 8, between the Si planes. Each clay sheet contains two tetrahedral aluminum atoms in order to reproduce the ratio of tetrahedral substitution of the montmorillonite.43 The interlamellar space counts 14 sodium counterions. We add the water-water and water-cation potentials to the clay-water and clay-cation potentials. The waterwater potential is obtained by Clementi et al.41 from ab initio calculations on the water dimer. The water-cation
Modeling the Clay-Water Interface
1 100
Langmuir, Vol. 7, No. 3, 1991 551 8
c'z'
gcr)
1
-
I
II
I I
I ,
11,
z Figure 9. Local concentration (M) of the interlamellar water. The distance is measured in A. 5
10
15
5
10
15
r
Figure 10. Mean radial distribution function of the water molecules around the interlamellar cations. potential is calculated within the same approximation; its expression comes from the ab initio calculations of Enghydration shell appears clearly, while the next shells do strom et al.44 not emerge from the statistical noise. The first hydration We use only lateral periodic conditions because of the shell of the clay extends to 4.8 f 0.1 A from the Si plane. disymmetry induced by the clay sheets. The potential is It contains 136 f 5 water molecules, for a total surface of cut to zero when the radial distance between two inter2230 A2. The mean surface covered by a water molecule acting centers becomes larger than half of the length of is 16.4 f 0.6 A2, in agreement with the molecular surface the ce1145 (i.e. 15.54 A). of water (15 A2 ll). The mean concentration of water Since the water density in the vicinity of the clay surface between 5 and 15 A (thus out of the first hydration shells) may differ from the bulk water density, we do not know is equal to 49 f 6 M, close to the bulk water concentration. really the total number of water molecules present in the Figure 10 shows the mean solvation of the Na+ couninterlamellar space. One must begin Monte Carlo46 simulations in the grand canonical e n ~ e m b l e ~ ~to? ~ ~ @terions. The radial distribution function exhibits a first narrow hydration pick, followed by a broad secondary pick. determine the mean number of water molecules. The No more structure raises from the noise, and the districhemical potential of free water is given by p / K T = bution function decreases gradually to 0.5. The first shell -11.1728.49 With the number of cations being constant, only their moves are taken into account. The maximum contains 5.68 f 0.02 water molecules per sodium cation. size of the moving is adjusted to maintain an acceptance To calculate the heat of immersion, we use a thermorate of 50%. To improve the convergence and avoid dynamic cycle (see Scheme I), splitting the immersion correlation, we use block averages50with a size of the block process in three steps equal to 20 times the number of particules. To avoid adding a water molecule in a forbidden position, AHi, = AHl4= AHl2 AHz3 AH34 (6) we introduce minimal distances of 2.4 A between the The first contribution, A H 1 2 , sums the enthalpy of water oxygen of the water molecule and the components of the stabilization in the bulk plus the enthalpy of swelling of clay, 2.2 A between two oxygens of water molecules, and a dry clay. The second term, A H 2 3 , is the heat of 1.6 A between the oxygen of water and the sodium cation. dissociation of the dry interface. This heat results from We start the simulation with dry clay and stop the grand the electrostatic interaction between all charged species: canonical simulations when the number of water molecules the two clay surfaces and the Na+ counterions. the last is constant. This gives an initial configuration of the water term, A H 3 4 , is the heat of formation of the same interface molecules in equilibrium with the temperature of the in the presence of water. The interlamellar distance system (298 K). For an interlamellar distance of 20 A we between the swollen clay sheets is set equal to 20 A. obtain 476 water molecules. We then calculate the average The Monte Carlo simulation gives the enthalpies A H 2 3 of the desired properties (energy, water distribution, etc.) (28.5 f 0.1 MJ/mol) and A H 3 4 (-63.7 f 0.2 MJ/mol). The in the framework of the canonical ensemble. This is the heat A H 1 2 contains two contributions: the stabilization of longest step of the calculation. The convergence of the results requires 1 500 000 configurations. It needs a water and the swelling of the dry clay. This last term calculation time of 10 days on an Alliant FX20 supercontains two new contributions: the long range electrominicomputer. static interaction of the interface and the short range van der Waals interaction between the clay surfaces. We We simulate the solvation of the clay surface and the interlamellar counterions. Figure 9 shows the distribution estimate the stabilization of the water molecule from the of the water molecules in the interlamellar space. A first enthalpy of vaporization of the 476 water molecules (19.42 MJ/mol). This approximation assumes that the steam behaves like a perfect gas. We have already calculated (44) Enstrom, S.; Jonsson, B.; Jonsson, B. J.Magn. Res. 1982,50,1-20. (45) Pangali, C.; Rao, M.; Berne, B. J. Mol. Phys. 1980,40,661-680. the electrostatic c ~ n t r i b u t i o nfor ~ ~an interlamellar space (46) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. filled by a continuum with a dielectric constant equal to H.; Teller, E. J. Chem. Phys. 1953,21, 1087-1092. (47) Adams, D. J. Mol. Phys. 1974,28, 1251-1252. 78.5. The electrostatic contribution is thus multiplied by (48) Adams, D. J. Mol. Phys. 1975,29, 307-311. the same dielectric constant (2 KJ/(mole site)36 X 14 X (49) Thermodynamique Chimique; Prigogine, I., Defay, R., Eds.; De78.5 = 2.2 MJ/mol). Finally the van der Waals interaction soer: LiBge, 1950. (50) Bishop, M.; Frinks, S. J. Chem. Phys. 1987,87, 3675-3676. between the two clay sheets is estimated by the QPEN
+
+
Delville
552 Langmuir, Vol. 7, No. 3, 1991
Table VI. Thermodynamic Data of Some Smectitesa
Scheme I
L
1
~1009 M2509 M7009 laponite" hectoritell bentonite66
I
II 1I
0
swollen h y d r a t e d c l a y
t
0 0 steam
I
w
I 2
t
dry clay
\
t
mode15l~~~ obtained from ab initio calculations in the watersilicate system. This model takes into account the loose and bound electronic pairs of the oxygen at the surface of the clay. We use the parameters of the water-silicate interaction to describe the silicate-silicate interaction. The contribution to A H 1 2 is small (1.7 MJ/mol). The sum of these three contributions gives A H 1 2 = 23.32 MJ/mol. The final heat of immersion is (-11.9 f 0.3) MJ/mol or (-0.89 f0.02) J / m 2since the clay fragments have a surface of 2230 A2. This corresponds to a heat of 712 J / g for a specific surface of 800 m2/g. This result is in agreement with the experimental data published on semctites (see Table VI). It is also possible to calculate the residual anisotropy of the orientation of the water molecules induced by their interaction with the clay surface. Indeed, deuterium NMR of diluted suspension of montmorillonite (Ecca gum BP, English China Clay, St. Austell, Cornwall) shows, in heavy water, a splitting of the signal of ~ a t e r . *This ~ , ~splitting ~ is due, on one hand, to the alignment of the clay sheets by the magnetic field and, on the other hand, to the orientation of the water molecules by the clay surface. It is possible to separate these two organizations since their characteristic time scales are completly different. The observed splitting is given by53 16 = 0.75xA(3 COS' BLD - 1)
S is the specific surface and
112.65 81.594 16.804 416 75 90 S i m
0.767 0.786 1.192 1.2 0.96
the immersion enthalpy.
in heavy water (215 kHz54)and A the residual anisotropy obtained from the mean orientation of the water molecules. This parameter varies between-1 and +l. OLD is the mean angle between the magnetic field and the normal to the clay surface. The anisotropy of the deuterium of heavy water is given by53
A = 0.05(3 COS' BDM - 1) - 0.41(sin2ODM COS 2 4 D M )
steam
t
146.9 103.8 14.1 345.1 79
(7)
where x is the quadrupolar coupling constant of deuterium (51) Sauer, J.; Morgeneyer, C.; Schroder, K. P. J . Phys. Chem. 1984, 88, 6375-6383.
(52) Skipper, N. T.; Refson, K.; McConnell, J. D. C. Clay Miner. 1989, 24, 411-425.
(53) Halle, B.; Wennerstrom, H. J . Chem. Phys. 1981, 75, 1928-1943.
(8)
where BDM and DM are the Euler angles defining the orientation of the water molecule with respect to the clay surface.53 The result is averaged over all the orientations of the water molecules. For a free water molecule, the thermal motion will average the anisotropy to zero. With 476 water molecules we obtain a splitting of (1.3 f 0.1) kHz, while the measured quantity is 10 HzZ3rz5 for a diluted suspension (0.02 g of clay for 1 mL of heavy water). Only water molecules in contact with the clay surface contribute to the observed splitting. The residual splitting is proportional to the fraction of water in contact with the clay surface. For a totally swollen clay (800 m2/ g) we obtain A6 =
(1.6 X loz1) X 476 X 1300 Hz = 13 (2.23 X 103)(3.347X 10")
HZ (9)
since 0.02 g of clay develops a surface of 1.6 X loz1A2, our model of clay has a surface equal to 2230 A2, and 1 mL of molecules. water contains 3.347 X This result must be multiplied by the angular coefficient (3 cos2ODL - 1)to take into account the orientation of the clay sheets in the magnetic field. We must also divide eq 9 by the mean number of clay platelets per tactoid, since the NMR is sensitive only to the adsorbed water molecules in fast exchange'3~~~ with bulk water. So it does not detect water molecules trapped in the tactoids.15 The mean number of clay sheets per tactoid is lower than 2 for sodium m ~ n t m o r i l l o n i t e . The ~ ~ , ~angular ~ coefficient may be equal to 2 or 1. The first value corresponds to a diamagnetic sheet, aligned parallel to the magnetic field, while the second value corresponds to a paramagnetic sheet, aligned perpendicular to the field. So the total correction varies between 0.5 and 2. The calculated value 13 f 1 Hz may be satisfactorily compared to the observed value of 10 Hz. 111. Discussion
The complexity of the problem precludes the use of an ab initio method for the MO calculations of the clay fragments. Ab initio calculation are feasible only for very small f r a g m e n t ~ . 5 l ,A~ semiempirical ~ approach must be used. For that reason we have selected the MIND0 level. The same calculations were performed with the extended (54) Waldstein, P.; Rabideau, S. H.; Jackson, J. A. J . Chem. Phys. 1964, 41, 2407-3411.
(55) Schramm, L.; Kwak, J. C. T. Clays Clay Miner. 1982,30,40-48. (56) Kassab, E.; Seiti, K.; Allavena, M. J . Phys. Chem. 1988,92,67056709.
Modeling the Clay-Water Interface
Huckel method.57 The number of heavy atoms in our fragments reaches the upper limit for the use of this program in reasonable delays (1 h of CPU time on the minisupercomputer SCS40). The MINDO local charges of the tetrahedral components of the clay (see Table 11: 1.7 au for Si and between -0.6 and -0.9 au for 0) are in agreement with the available results obtained from ab initio calculations on zeolite56 (1.6-2.3 au for Si and -0.8 to -1.0 au for 0). A statistical treatment of the clay-water interface requires simple analytical approximations of the calculated clay-water and clay-counterion interactions. This is the goal of the energy parametrization. The basic hypothesis of this approach is the pairwise additivity of the atomic potentials. It is the basis of any molecular mechanic58159 treatment. The two last lines of Table I gives some details on the quality of the fitting. The MINDO minimum energy is -13.3 kJ/mol for the neutral fragment and -27.3 kJ/mol for the charged fragment. The fitted values are -11.3 and -26.7 kJ/mol, respectively. To check the additivity of the potential, we introduce in the fitting of the energy two fragments containing four octahedral aluminum (A14Si6030H24and A15Si5030H24-). Only one configuration of the water molecule is used. We take the optimum configuration of the smaller fragments. The agreement between the MINDO and fitted energies is encouraging (see Table I). The localization of the water minima at the center of the hexagonal cavities has already been suggested from an analysis of IR data.lg Indeed the water molecule may form an hydrogen bond with the lone pair of the oxygen at the center of the cavity. This is possible for dioctahedral clays,lgsince the proton is directed parallel to the plane of the clay. For trioctahedral clays, the proton is perpendicular to the clay surface. The optimization of the water energy requires another orientation. The water molecule must point a lone pair of its oxygen in the direction of the proton located a t the center of the cavity.lg In this case, the approach of asingle water molecule induces a reorientation process. At large distance, the optimum configuration is the same: the two protons pointing toward the clay surface. The set of parameters allows the determination of the optimum position of a water molecule or a sodium counterion near to a clay surface. This work is equivallent to a molecular mechanic58p59 approach and has already been published for the interaction of water with a neutral clay.52 However, such studies are limitated, since they treat only one adsorbed molecule. They do not model the organization of a dense solvent near the surface of the support. Another limitation results from the principle of molecular mechanic. It searches for an absolute minimum of the energy instead of averaging the results over all the available configurations of the system. Only a statistical model will establish a junction between calculated energies for a diluted system (only two interacting particles) and measurements on dense phases. Our model of the clay-water interface is quite heterogeneous. It mixes the potentials obtained from different levels of approximation. This will have no drastic effects, since we need only to know the energy differences. The (57) Aronowitz, S.; Coyne, L.; Lawless, J.; Rishpon, J. Inorg. Chem. 1982,21, 3589-3593. (58) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S.; Weiner, P. J. A m . Chem. SOC.1984, 106,765784.
(59) Allinger, N. L.; Yuh, Y. H.; Lii, J. H. J. Am. Chem. SOC. 1989,111, 8551-8566.
Langmuir, Vol. 7, No. 3, 1991 553
absolute value of the energy is more sensitive to the level of approximation. Another problem in Monte Carlo simulations concerns the c ~ n v e r g e n c e . One ~ ~ needs to obtain results independent of the size of the system. The Lennard-Jones potential has a finite range and the use of the minimum image procedure with periodic boundary conditions60 suffices to avoid border artifact. However,the electrostatic potential, because of its long range, induces a divergence in the calculation of the energy and of the distribution of the interlamellar cations.6l A classical p r o ~ e d u r e ~ ~ s ~ l ~ ~ consists of adding an external potential in order to compensate for the long tail of the Coulomb potential lost in the cutoff. This procedure ensures the convergence of the electrostatic energy resulting from the dipoles38of the water molecules. It also ensures the convergence of the distribution of the counterions even when the interlamellar space contains only 20 ions.36 But the convergence of the electrostatic energy of the counterions requires at least 400 to 500 ions.36 Thus one must be very careful in the interpretation of results obtained by a simulation of only 14 ions. The absolute energy is certainly wrong, while the distributions of the counterions and a fortiori the water molecules are reliable. The simulation of the clay solvation involves a great number of particles, but this approach is powerful. As an example, it distinguishes between the condensed and the bond cations. Indeed, a condensed counterion holds an intact first solvation sphere. A bound counterion must lose a t least one water molecule from its first solvation sphere to coordinate to a site at the surface of the clay in a so-called specific interaction. The simulation of the solvation also avoids any problem related to the decrease of the local dielectric constant in the vicinity of the surfa~e.6~~65 We assume t, = 1 everywhere. The orientation of the water molecules, pointing their dipole parallel to the electric field, reduces the effective field. Figure 9 shows the organization of the solvent induced by the molecular surface. It is much lower than the organization obtained with a smooth This result is due to interference between different sites a t the surface of the clay, generating a more diffuse organization than a smooth surface.34 The mean number of water molecule in the first solvation sphere of the sodium counterion (5.68 i 0.02) deviates significantly from the hexaaqueous sodium.44It corresponds to the maximal bonding of four cations. This number is equal to the number of tetrahedral substitutions. The analysis of spectral measurements of the adsorbed water25 yields the same conclusion. The tetrahedral aluminates bear an excess of negative charge (-0.88 au). They are ideal candidates for site binding. However, to confirm this result, one needs to analyze the composition of the first hydration shell as a function of the position of the cation in the interlamellar space. However, to obtain results with statistical significance, one must introduce more ions in the simulation and thus increase the number of water molecules. The resulting comptuting time will become prohibitively long. Nevertheless, specific clay-sodium interaction may be (60) Wood, W. W.; Parker, F. R. J. Chem. Phys. 1957,27, 720-733. (61) Jonsson, B.; Wennerstrom, H.; Halle, B. J.Phys. Chem. 1980,84, 2179-2185. (62) Torrie, G. M.; Valleau, J. P. J. Chem. Phys. 1980, 73,5807-5816. (63) van Megen, W.; Snook, I. J . Chem. Phys. 1980, 73, 4656-4662. (64) Hunter, R. J. Zeta Potential in Colloid Science, Principles and Applications;Ottewill,R.H., Rowel1,R. L.,Eds.;Academic Press: London, 1981; Chapter 2. (65) Torrie, G.M.; Valleau, J. P.; Patey, G. N. J. Chem. Phys. 1982, 76, 4615-4622.
554 Langmuir, Vol. 7,No.3, 1991
Deluille
by Monte Carlo simulations. Their electrostatic contributions are exactly opposite and thus vanish. Many additional hypotheses are necessary to obtain the desired macroscopic information. First, we calculate the heat of immersion with a maximum interlamellar distance of the swollen clay equal to 20 A. This reduces the total number of water molecules filling the interlamellar space and is allowed since the immersion enthalpy decreases strongly with the water content of the clay. It is reduced by 1 order of magnitude for a clay interface containing three solvation layers.ll*ss Second, the estimation of the electrostatic contribution to clay swelling ignores the variation of the ionic condensation with the dielectric constant. This hypothesis is certainly well founded a t small interlamellar distances, since the ionic condensation disa~pears.3~ This is less justified a t large interlamellar distances. Finally, the small X range contribution of the silicate-silicate interaction to the swelling of the dry clay is evaluated from a watersilicate potential. Thus, the calculation is useful as long as the error does not exceed 10% of the final result. This is verified owing to the smallnessof these two contributions (2.2 and 1.7 MJ/mol, respectively). Finally, we use the Monte Carlo simulations to calculate the residual anisotropy of the water orientation near the clay surface. This calculation implicitly assumes that the mean orientations of H20 and D20 are the same. What is the origin of the clay swelling? We have now all the information necessary to understand this question and quantify its driving forces. A generally recognized h explanation of the long range swelling is the electrostatic Y repulsion between the clay sheets bearing the same charge. This view is partial since it neglects the contribution of the interlamellar counterions. One must treat the system formed by the charged surfaces plus the counterions as a I 2 whole. This system is neutral and attractiue. This is the Figure 11. First hydration shell of a sodium (*) near the result of our Monte Carlo simulation^^^ and also of the hexagonal cavity including a tetrahedral aluminum (a). observation. Without water, the clay sheets do not swell but remain collapsed. Ionic crystals have the same shown. As an example, Figure 11 illustrates one configproperty. Their stability is explained by the Madelung uration of the Monte Carlo simulation displaying only the constant resulting from their electrostatic interactions. first hydration shell of one Na+ counterion in the vicinity So the electrostatic interaction is not a driving force of the of the clay surface. The cation shown in Figure 11 is above clay swelling! a hexagonal cavity containing a tetrahedral aluminum. Let us subtract the electrostatic contribution from the Its first hydration shell contains five water molecules. One energy of the clay. The residual energy of the hydrated water molecule is pinched between the cation and the clay clay with an interlamellar distance of 20 A is -35.2 MJ/ surface. To optimize its energy, it directs one hydrogen mole. It contains water-water interaction (-9.7 MJ/mol), to the center of the cavity and one lone pair of the oxygen clay-water interaction (-8.2 MJ/mol), and water-sodium to the Na+ cation. This is a nice example of interference interaction (-17.3 MJ/mol). This result appears erroneous between the solvation of the clay sheet and the interlamelsince the water-water interaction (-9.7 MJ/mol) is much lar counterion. lower than the water stabilization in the bulk (-19.42 MJ/ The radial distribution function of the water molecules mol). On the other hand, the ion-water interaction (-17.3 around the interlamellar counterions (Figure 10) shows MJ/mol) is much larger than the heat of solvation of the two intriguing properties. First, the secondary pick is 14 sodium cations (-5.8 MJ/mol). This discrepancy comes broader than the first one, and second, the radial distrifrom an unfounded comparison between the ion-water bution function does not show a plateau at 1 but decreases enthalpy and the heat of solvation of the ions. Indeed, to 0.5. The broadening of the secondary pick may be due this last term contains not only the ion-water interaction to the average of results independently of the position of but also water-water interactions between two water the cation. The solvation of the condensed counterions molecules solvating the same cation. This last interaction may differ from the solvation of the counterions far from may be repulsive. Moreover, the solvation of the cation the clay surface. The decrease of the function is due to interferes with the structure of water. The final balance some obstruction from the clay surface, reducing the is a reduction of the enthalpy of solvation. number of water molecules organized around the counThe detail of the energy allows us to identify two driving terions. forces responsible of the clay swelling: the solvation of The calculation of the immersion enthalpy requires a the clay sheets which accounts for 32 % of the energy and thermodynamic cycle, because of the long range of the the solvation of the cation which accounts for 68% of the electrostatic interaction. The size of the simulated system is too small to ensure the convergence of the electrostatic (66) Zettlemoyer, A. C.; Young,G. J.; Chessick, J. J. J. Phys. Chem. contribution. The two heats A H 2 3 and A H 3 4 are calculated 1955,59,962-966.
yl
w
Modeling the Clay-Water Interface
Langmuir, Vol. 7, No. 3, 1991 555
Table VII. Influence of the Nature of the Cations on the Size of the Tactoids
K+
no. of sheets per tactoid 1 1.5 2
Rb+ cs+
2.6 3
cation
Li+ Na+
A&+,
kJ/mol -530 -415
-320 -300 -217
energy. This explains why the nature of the interlamellar cation modifies the swelling properties of the clay (see Table VII). From lithium to cesium, the mean size of the t a ~ t o i d s ' increases ~,~~ by a factor of 3 while the heat of solvation of the cation6' decreases by a factor of 2. To confirm the influence of the counterion on the swelling property of the clay, we have performed an additional simulation for a neutral clay without any interlamellar counterion. The resulting A H 1 2 equals 21.12 MJ/mol, A H 2 3 equals 0, and A H 3 4 equals -17.27 f 0.1 MJ/mol. The immersion is endothermic (Mi,,,= +3.85 f 0.1 MJ/mol). This explains why talc, which does not possess interlamellar counterions, interacts poorly with water. Of course this result is punctual. It derives from simulations performed for one model of clay, in one particular configuration, with only one type of ion. From a qualitative point of view, the conclusions will probably be the same for any smectites and perhaps other colloidal systems. The electrostatic energy is repulsive; the driving forces of the swelling are the solvation of the counterions and of the clay surface. Of course, the relative contribution of these different terms will vary with the nature of the clay surface and of the interlamellar cations. The next question is: what is the driving force of the clay swelling at large interlamellar distance? These solvation interactions are only short range interactions. Their contributions are negligible when the clay interface contains more than three hydration shells. Adding water to such system does not modify the hydration shells of the clay surface and of the interlamellar counterions. The factor responsible of swelling is not enthalpic but entropic! Let us imagine a suspension of monoionic clay (67)Rashin, A. A,; Honig, R. J . Phys. Chem. 1985,89,5588-5593.
separated from pure bulk water by a semipermeable membrane. The second law of thermodynamic wants to equilibrate the chemical potential or the activity of any exchangeable species. The activity of the compensating cation is zero in the bulk water and positive in the clay suspension. The cations must move from the clay suspension to the bulk water to equalize their activity. However, this is impossible in the absence of an external electric field, because the clay suspension must remain neutral. So the water molecules of the bulk phase will pass through the membrane to dilute the clay suspension and reduce the activity of the counterions. The van't Hoff law68 describes this entropic effect and provides a good description of the swelling property of dilute clay susp e n s i o n ~ .This ~ ~ ~entropic ~~ contribution does not exist for a dry clay, since its counterions are not solvated exchangeable species. They must be considered as structural parts of the clay interface .
IV. Conclusions This work is a first attempt a t a molecular description of the solvation of the clay surface and the interlamellar counterions. It gives a good description of the immersion enthalpy and the mean water organization of sodium montmorillonite. The results of the simulations help to understand and quantify the driving forces responsible of clay swelling. This approach must be generalized to other systems, incorporating organic molecules. A molecular model of heterogeneous catalysis will also provide a better understanding of the physicochemical processes involved: physisorption, surface diffusion, chemisorption, and desorption. Acknowledgment. We thank Professor P. Laszlo for his interest during the progress of this work and Professor J. Y. Lallemand (Laboratoire de SynthBse Organique, Ecole Polytechnique, Palaiseau) for access to the Alliant FX20 minisupercomputer used for the Monte Carlo simulations. We thank the SCS Co., Paris, for the loan of the SCS40 minisupercomputer used for the MIND0 calculations. (68)Langmuir, I. J . Chem. Phys. 1938,6,873-896. (69)Delville, A.;Laszlo, P. Langmuir 1990,6,1289-1294.
