J. Phys. Chem. B 2000, 104, 9163-9170
9163
Simulations of Clay Mineral Swelling and Hydration: Dependence upon Interlayer Ion Size and Charge Daniel A. Young and 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 11, 2000; In Final Form: July 19, 2000
The dependence of the crystalline swelling and hydration properties of clay minerals on interlayer ion size and charge was investigated using molecular dynamics computer simulations of Na-, Cs-, and Sr-substituted montmorillonites. For all clays studied, layer spacings measured as a function of water content exhibited plateaus at the one-layer hydrate spacing. Calculated immersion energies exhibited minima for integer-layer hydrates up through the three-layer hydrate, with apparent global minima identified with the one-layer hydrate for Cs-montmorillonite and with the two-layer hydrates for the Na- and Sr-montmorillonites. In addition, for Sr-montmorillonite, layer spacings jumped discontinuously between one-layer and two-layer separations and showed a second plateau at the two-layer hydrate spacing. The immersion energy curve for Sr-montmorillonite showed similar discontinuities. These results provide evidence for a constant water content swelling transition between one-layer and two-layer spacings. This transition was further characterized using pressure versus layer-spacing isotherms. The isotherms showed a loop structure with peaks that correlated with the coordination number of the interlayer ion. Integration of the Sr-clay isotherms yielded Gibbs energy curves with minima at both one-layer and two-layer spacings. The results indicate that the mechanism of swelling and hydration depends upon the interlayer ion charge. For monovalent ions studied, swelling and hydration were coincidental, with water occupying all available interlayer volume. In contrast, expansion of Sr-montmorillonite, driven by formation of a first and partial second hydration shell for Sr2+, occurred initially without the filling of the interlayer volume. Trapping in metastable expanded states was also observed for Sr-montmorillonite and identified as a possible mechanism for layer-spacing hysteresis.
1. Introduction Montmorillonites and other smectite clay minerals are abundant, layer silicates that are notable for their tendency to expand and hydrate upon exposure to water. The swelling process is driven in great part by hydration of exchangeable cations that reside between the negatively charged clay layers. Smectites play an important role in many environmental and engineered systems. For example, swelling and ion exchange properties of smectites impact the transport and bioavailability of both ionic nutrients and pollutants in the environment.1,2 These same properties make swelling clays useful as backfill materials in nuclear waste repositories.3,4 Smectite swelling also leads to unwanted consequences such as structural damage5 and the degradation of borehole stability in petroleum drilling operations.6 Two regimes of swelling have been observed in smectite clay minerals.7 The first is the short-range or crystalline regime where swelling occurs discretely through the formation of integer-layer or mixtures of integer-layer hydrates. Water contents in this regime range between zero and approximately four water layers. Under certain conditions, a second regime involving much more extensive expansion and hydration may be observed. The equilibrium hydration state of any smectite clay is known to depend upon the magnitude and location of the clay layer charge, the temperature, the applied pressure, and the water chemical potential.8-12 In addition, the size and charge of the interlayer ions significantly influence swelling and hydration properties.7,13-18 * Corresponding author. E-mail:
[email protected].
Although the extents of hydration and expansion have been determined for a wide range of smectites, fundamental issues regarding the structural and thermodynamic origins of swelling still remain unresolved. Two of these issues, related directly to the simulations reported here, involve hysteresis and hydration mechanisms. First, hysteresis in both water content and interlayer spacing is almost universally observed in swelling measurements.12,15,16,19 This makes the identification of equilibrium states difficult and therefore restricts the usefulness of equilibrium thermodynamic expressions related to smectite swelling.20 Capillary condensation is a likely cause of hysteresis in external surface hydration but cannot account for layer-spacing hysteresis, and its origins remain unclear.12,19,21 Second, there is evidence that interlayer expansion and hydration may not be coincidental. It is commonly observed that X-ray diffraction layer spacings, measured as a function of water vapor pressure, show plateaus corresponding to broad regions of nearly constant layer spacing. In contrast, gravimetric water contents, measured as a function of water vapor pressure, increase much more continuously.13-16 This has led to the suggestion that water adsorption occurs in two distinct steps involving the formation of the hydration shell of the cation followed by the filling of the remaining interlayer space.15-17 The first step leads to layer expansion and some water adsorption, while the second step involves an increase in water content without substantial swelling. Evidence for this proposed mechanism depends directly upon the determination of interlayer and external surface components of the gravimetric water content. The validity of this mechanism has been challenged recently by Laird22 using
10.1021/jp000146k CCC: $19.00 © 2000 American Chemical Society Published on Web 09/01/2000
9164 J. Phys. Chem. B, Vol. 104, No. 39, 2000
Young and Smith
TABLE 1: Solution Simulation Parameters element
σ/Å
O H Cs Na Sr
3.166 0.000 3.830 2.350 3.341
/kJ‚mol 0.650 0.000 0.418 0.544 0.418
-1
q/e -0.848 0.424 1.000 1.000 2.000
evidence that the interlayer water content in Mg-substituted clays correlates well with the available interlayer volume. He suggests that hydration and expansion occur in a concerted fashion, and that hydration of external surfaces accounts for the discrepancy between layer-spacing and water-content measurements. Molecular computer simulations can be a valuable tool in addressing unresolved issues related to smectite swelling. Numerous models for the simulation of expandable clay minerals have been developed in the past decade and applied with at least qualitative success to investigations of the structure and thermodynamics of clay swelling.23-28 Although simulation models are not perfect, computational methods do have several unique capabilities and thus provide a complementary tool to experimental methods. Direct access to structural information in simulations allows for the investigation of structure-function relationships. In addition, nonequilibrium (metastable) states associated with swelling hysteresis may be investigated systematically using computational methods. Finally, simulations allow for a study of complex clay systems under simplified conditions. For example, information on the behavior of individual interlayers can be obtained without complications from mixed-layer hydrate formation or interference from external surface effects.29 In this article, results from simulations of Na-, Cs-, and Srsubstituted montmorillonites in the crystalline swelling regime are presented. Cesium and strontium were chosen in part because of their abundance in high-level nuclear waste, and sodium was chosen as a reference since sodium smectites are widely studied. Comparison of results for the three ions allows for analysis, albeit incomplete, of the ion size and charge dependence of swelling and hydration thermodynamics. Results of both layerspacing and immersion-energy calculations were found to be consistent with the discrete nature of crystalline swelling. Trends with ion size and charge indicated correctly that cesium acts as a swelling inhibitor when compared with either sodium or strontium. The hydration mechanisms for monovalent and divalent ion-exchanged montmorillonites were found to be distinctly different. The Sr-clay alone showed a clear two-step filling process, suggesting that the hydration mechanism may be dependent upon the interlayer ion hydration energy. Twostep hydration processes such as that observed in Sr-montmorillonite imply the existence of metastable swelling states. Trapping in these states provides one reasonable mechanism for layer-spacing hysteresis. 2. Computational Methods 2.1. Interaction Potentials. The form of the solution and clay interaction potentials used here is equivalent to one described previously.29 Water was represented by the extended simple point charge (SPC/E) model of Berendsen and coworkers,30 while ions were modeled as charged Lennard-Jones (LJ) spheres. The LJ interaction parameters for water and ions are displayed in Table 1. The parameters for cesium were taken from a previous clay simulation,29 while those for sodium were from an ionic solution model.31 The strontium parameters were taken from a polarizable model study32 and were evaluated for
nonpolarizable systems using a 300 ps simulation of a SrCl2 solution. The solution contained one Sr2+ ion, two Cl- ions, and 429 water molecules. Parameters for nonpolarizable Clwere taken from a previous study.31 The solvation energy for SrCl2, calculated to be -2120 kJ‚mol-1, is in reasonable agreement with the experimental solvation enthalpy of -2201 kJ‚mol-1.33 The Sr-O peak position (2.65 Å) and coordination number (8.5) also agree well with the experimental values of 2.6 Å and 8, respectively.34 The potential energy of the interlayer solution is given by
U)
1
∑∑ 2 i j
{ [( ) ( ) ] }
′ 4ij
σij
12
-
rij
σij
6
+
qiqj
rij
rij
(1)
where qi is the charge on atom i, rij is the distance between atoms i and j, and the prime indicates the summations exclude intramolecular interactions. The Lennard-Jones cross-interaction terms, ij and σij, are generated from the standard combining relationships:
1 σij ) (σi + σj) 2
(2)
ij ) (ij)1/2
(3)
Montmorillonites, a subclass of smectites, are layered aluminosilicate minerals with a dioctahedral “2:1” structure.35 Each clay layer consists of a central sheet in which aluminum atoms are octahedrally coordinated with oxygens and two outer sheets of tetrahedrally coordinated silicon atoms. The sheets are connected through shared, apical oxygens to form layers approximately 10 Å thick and of typically micrometer lateral dimension. Isomorphic substitution of, for example, Mg for Al in the octahedral sheet or Al for Si in the tetrahedral sheets results in a net negative layer charge. This charge is balanced by interlayer cations. Clay minerals with this 2:1 structure, a total layer charge of between -0.5 and -1.2 e per unit cell, and the majority of the layer charge located in the octahedral sheet are classified as montmorillonites. The montmorillonite modeled here has a layer charge of -0.75 e and unit-cell formula
X0.75/n[Si8](Al3.25Mg0.75)O20(OH)4 where X is the interlayer ion (Na, Cs, or Sr) and n represents the charge on the ion. The clay is an idealized montmorillonite since all of the layer charge is located in the octahedral sheet. The clay-clay and clay-solution interaction potentials were developed from the models of Skipper and co-workers,23,25,36 modified to incorporate a Lennard-Jones (LJ) functional form for the van der Waals interactions.29 Tetrahedral substitution sites in this model have been associated with qualitatively incorrect ion exchange free energies37 and were therefore not included in this study. Each clay layer was treated as a rigid lattice with unit-cell coordinates given by Skipper and coworkers.25 Charge and LJ potential parameters for the clay are displayed in Table 2. Atoms in the octahedral sheet were not assigned LJ parameters since these primarily determine shortrange interactions and the interlayer water and ions never occupy the interior of the clay layer. 2.2. Simulation Methods. Simulations were performed using a periodically replicated system composed of two clay layers each containing eight unit cells and having x and y dimensions of 21.12 and 18.28 Å. Substitution sites in the octahedral sheet were chosen randomly, with the restriction that nearest-neighbor
Simulations of Clay Mineral Swelling and Hydration
J. Phys. Chem. B, Vol. 104, No. 39, 2000 9165
TABLE 2: Clay Simulation Parameters layer tet. apical oct.
element
σ/Å
/kJ‚mol-1
q/e
O Si O O H Al Mg
3.166 1.840 3.166 3.166 0.000 0.000 0.000
0.650 13.18 0.650 0.650 0.000 0.000 0.000
-0.800 1.200 -1.000 -1.424 0.424 3.000 2.000
sites could not both be substituted. Except for substitution sites, the two clay layers are mirror images, allowing for equivalent registry positions of the opposing siloxane surfaces in each interlayer region. The total charge of -6 e on each layer was balanced by interlayer Na+, Cs+, or Sr2+ ions. The z-axis, taken to be perpendicular to the clay layers, had a variable length depending upon the interlayer water content and applied pressure. Ions and water molecules were distributed evenly between the two interlayer regions in the simulation cell. The volumes of the two regions were also constrained to be equal as the z-axis length varied. All simulations were performed using the SOLVENT code described previously.38 Trajectories for the interlayer solution were calculated using the velocity-Verlet algorithm39 with a time step of 1.5 fs. Long-range, Coulombic interactions were treated using Ewald summation methods,40,41 whereas the minimum image convention was used for all short-range terms. Water bond lengths were constrained using the RATTLE42 variant of the SHAKE43 algorithm. Following equilibration, all calculations were performed at a constant temperature of 298 K using a weak-coupling algorithm and a coupling time constant of 500 fs.44 Individual calculations were performed at either constant layer spacing or constant pressure, where pressure refers to the component of the pressure tensor normal to the clay sheets. For constant-pressure runs, the z-axis length was allowed to fluctuate using a weak-coupling pressure algorithm and a target pressure of 1 bar.44 For all runs, the pressure was calculated from the volume dependence of the system potential energy using a finitedifference algorithm.28 For selected water contents, pressure versus layer spacing isotherms, P(z), were calculated as a function of the layer spacing z using a series of constant-volume simulations in which the average pressure was measured. Registry motions of the clay sheets were accomplished using a standard Monte Carlo (MC) algorithm. Random lateral moves of between -0.1 and 0.1 Å in the x and y directions were attempted for one clay sheet every five simulations steps, with acceptance of the moves based on the Metropolis criteria.40 Acceptance ratios ranged between 0.4 and 0.8 for the various systems studied. Twisting motions were not attempted as they disrupt the rectangular symmetry of the simulation box. Strong interactions of hydrated Sr2+ complexes with the clay surfaces were found to significantly limit registry movement, making complete sampling over an extended registry range intractable. To address this concern, the range of registry motion in some simulations was restricted by rejecting moves that would place the clay outside of a predefined registry region. This region was chosen to be triangular, as illustrated in Figure 1, and is the minimum area required to sample all registry configurations of the opposing siloxane surfaces. The right point of the triangle corresponds to complete overlap of the clay hexagonal cavities. In the upper left point of the registry area, the hexagonal cavities of each clay sheet are across from tetrahedral silicon sites on the opposing surface. The lower left corner corresponds to overlap of hexagonal cavities with siloxane oxygen atoms. A much larger range would be required to sample all configura-
Figure 1. Schematic diagram of the allowed registry configurations in restricted registry simulations. Configurations outside the triangle are forbidden. Vertices on the hexagonal lattice correspond to the silicon atoms on a portion of one siloxane surface. The filled circles represent a ring of silicon atoms on the opposing surface. Coordinates give the relative positions, in angstroms, of the two clay surfaces. Vertices of the triangle correspond to configurations in which the hexagonal cavities of one clay surface overlap with (i) the hexagonal cavities of the opposing clay surface (2.66,0), (ii) silicon tetrahedra of the opposing clay surface (0,1.53), and (iii) oxygen atoms of the opposing clay surface (0,0).
tions with respect to the clay hydroxyl groups and isomorphic substitution sites. However, comparison studies indicate little sensitivity to registry restrictions. The current approach lends itself well to graphical analysis of sampled registry configurations. Starting configurations for each interlayer solution were generated from existing cesium montmorillonite configurations containing either 96 or 160 waters in each layer.29 The initial configurations were all equilibrated for at least 60 ps and until no measurable drift in the system energy was observed. Configurations with less water were then generated by a “desorption” process in which randomly chosen water molecules were removed from the interlayer regions of either the 96 or 160 molecule systems, followed by equilibration for at least 30 ps. Both registry and ion motions in the strontium-substituted clay are slow, particularly at low water contents, relative to motions in the other clays. An annealing procedure that enhances sampling of these degrees of freedom was therefore used to generate some starting configurations for Sr-montmorillonite calculations. The procedure consisted of performing a 15 ps simulation at 500 K followed by cooling at a constant rate to 298 K over a 30 ps period. This was followed by equilibration for 15 ps. Runs which follow this procedure are labeled as annealed, while those with only 298 K equilibration are labeled nonannealed. Constant-pressure simulations were performed for Na-, Cs-, and Sr-montmorillonites with water contents ranging between N ) 24 and N ) 160 molecules. Averages were collected for 300 ps in each case, with uncertainties calculated using a blockaveraging technique and a 30 ps block size.29 Here and throughout N corresponds to the number of water molecules per interlayer region. For reference, N ) 44 in Cs-montmorillonite corresponds to a gravimetric water content of 0.121gH2O/ gclay. For larger values of N, the Cs- and Sr-montmorillonite water contents exceed the maximum values observed experimentally. These “unphysical” states were investigated in order to better understand why swelling is inhibited by some interlayer ions. Under some conditions, there may be multiple stable states for a system at a single pressure. The classic example of this
9166 J. Phys. Chem. B, Vol. 104, No. 39, 2000
Young and Smith
behavior is phase coexistence. Under such conditions, constantpressure simulations may produce nonunique results that depend upon the initial system configuration.28,45 As will be discussed below, this was observed for Sr-montmorillonite in the N ) 36-48 region. The origin of this behavior was investigated systematically by measuring pressure versus layer-spacing isotherms for Na- and Sr-montmorillonites at water contents of N ) 44, 48, and 56. For each water content, a series of simulations at constant layer spacing was performed, with spacings constrained to values ranging between 12.0 and 15.0 Å. Starting configurations at each spacing were annealed as in the constant-pressure runs and averages were then collected for 75 ps. No hysteresis was observed when starting configurations were generated sequentially either from smaller to larger or from larger to smaller spacings. 2.3. Swelling Thermodynamics. For each clay, the thermodynamics of hydration was characterized in terms of the immersion energy, Q, defined as
Q ) 〈U(N)〉 - 〈U(N°)〉 - (N - N°)Ubulk
(4)
where 〈U(N)〉 and 〈U(N°)〉 are the average internal energies of the clays containing N and N° water molecules per interlayer, respectively, and Ubulk ) -41.4 kJ‚mol-1 is the internal energy of bulk SPC/E water.31 All energies were adjusted by the selfpolarization energy implicit in the SPC/E model.30 The water content of the reference state, N°, was chosen such that the minimum value of Q is equal to zero. The immersion energy therefore corresponds to the energy released when a clay with interlayer water content N transfers water to or from bulk water to form the clay in its reference state. Oscillations in the immersion energy provide a signature for discrete, crystalline swelling processes and a means for identifying integer-layer hydrate states.29 The relative stabilities of multiple stable states in a clay system under conditions of constant T, N, and applied pressure, Pext, may be determined from the pressure versus layer-spacing isotherm, P(z). All stable and metastable states may be identified by points where P(z) ) Pext and where the slope of P(z) is negative. The Helmholtz energy for swelling at fixed N and T is then given by
∫
∆A(z) ) -A P(z) dz
(5)
where A is the cross-sectional area of the clay. Integration was performed using the trapezoid rule. For zero applied pressure, the Helmholtz and Gibbs energies are equal, and the minimum in ∆A(z) corresponds to the most stable state. For nonzero applied pressure, the Gibbs energy is given by ∆G(z) ) ∆A(z) + PextA∆z or, equivalently
∫
∆G(z) ) - A (P(z) - Pext) dz
(6)
The stable states correspond to minima in ∆G(z). Multiple stable states may occur as a result of oscillations in the calculated isotherm. This is analogous to the liquid-vapor phase coexistence that results from loops in the van der Waals equation of state. 3. Results and Discussion 3.1. Swelling Curves. Layer spacings for Na-, Cs-, and Sr-montmorillonites calculated as a function of water content are displayed in Figure 2. At high water contents, all of the swelling curves approach linear behavior with a slope that is consistent with the bulk water density. For clarity, therefore,
Figure 2. Simulated swelling curves for Cs-, Na-, and Srmontmorillonites.
data are displayed only out to N ) 120, a region encompassing the one- through three-layer hydrates. All three curves show a plateau in the one-layer hydrate region (z ) 12.0-12.5 Å), with spacings slightly larger for Cs than for Na and Sr. The transition from one-layer to larger spacings shows a significant dependence upon ion identity. Both Na- and Cs-montmorillonite exhibit a rapid yet continuous increase in layer spacing as water content is increased beyond the one-layer hydrate. The Sr-montmorillonite, in contrast, shows a discontinuous increase in spacing in the same region, with the location of the transition point dependent upon the preparation procedure. For Sr-clay systems prepared using the standard “desorption” process, with no annealing, the layer spacing remains a nearly constant 14.314.5 Å as the water content is reduced from 56 to 36, followed by a drop to a spacing of 12.2 Å at N ) 32. Annealing causes the transition point to shift upward to N ) 48. The results indicate that the two-layer spacing is favored at water contents intermediate between the one- and two-layer hydrates, and that the two-layer spacing may persist to very low water contents upon desorption. Comparison with experimental swelling curves is made difficult by uncertainties in the interlayer water content to use for comparison and by the formation of mixed-layer hydrates in the experimental systems.15 Nevertheless, experimental onelayer hydrate spacings are typically reported to be in the 12.012.5 Å range,14-16 in qualitative agreement with our results. The results for Na show somewhat smaller spacings than observed in previous simulations using both the TIP4P26 and MCY46 water models. This may be related to differences in clay-water interactions that were reformulated for use with the SPC/E model.29 The trend of increasing spacing with increasing ion size is opposite that reported by Boek and co-workers for Na- and K-montmorillonites.26 Although the origin of this difference is unclear, our result is reasonable given the greater ion volume and smaller electrostriction effect of cesium relative to sodium. 3.2. Immersion Energies. Calculated immersion energies for the three clays are displayed in Figure 3. All curves show oscillations indicative of crystalline swelling,29 with minima corresponding to energetically preferred integer-layer hydrates. Several trends with ion size and charge are evident. The magnitude of the immersion energies at the smallest water content (N ) 24) increases with increasing hydration energy (Cs+ < Na+ < Sr2+). This is in qualitative agreement with experimentally measured immersion enthalpies for dry clays.15,16 Also consistent with experimental measurements is the stability
Simulations of Clay Mineral Swelling and Hydration
Figure 3. Immersion energies for Cs-, Na-, and Sr-montmorillonites. Curves for Na and Sr are displaced for clarity, with dashed lines corresponding to the curve minima. Error bars indicate 95% confidence limits.
of the two-layer hydrate relative to the one-layer state. The twolayer hydrate is favored as hydration energy increases such that the one-layer state is preferred only for Cs-montmorillonite. Finally, the immersion energy minima shift outward in the same sequence. This suggests that the water content of a given integerlayer hydrate can depend significantly upon the interlayer ion identity. This behavior may be associated with a decrease in partial molar volume of the ions with increasing hydration energy,33 thus increasing the number of water molecules required to form an integer layer. Comparison of this trend with experimental measurements is made difficult by uncertainties in the interlayer water content in the experimental systems. Nevertheless, it appears that the integer-layer water contents predicted by the simulations are too high. In Sr-montmorillonite, for example, the water content at the two-layer minimum is equal to 96 molecules or 0.27gH2O/gclay. In contrast, experimental measurements yield a value of between 0.15 and 0.20 following subtraction of estimated external water contents.16 Although this discrepancy may reflect some limitation in the interaction model, we believe that it is due in great part to the neglect of entropic contributions in the immersion energy calculations. Inclusion of entropic effects results in a slight shift to smaller water contents in Cs-montmorillonite, with a swelling free energy reported recently28 showing a somewhat lower water content (N ≈ 41) for the one-layer hydrate than the immersion energy value (N ≈ 44) for the same clay, displayed in Figure 3. Recent calculations on a Mg-montmorillonite suggest that entropic effects are significantly larger in clays substituted with divalent cations.37 A grand canonical ensemble simulation on that clay yielded a two-layer hydrate spacing and water content in excellent agreement with experimental measurements.22 Work to determine swelling free energies for additional clays is in progress. As in the swelling curves, the immersion energy for Srmontmorillonite is unique in the one-layer hydrate region. The nonannealed Sr curve shows a discontinuous jump in Q at N ) 32, whereas the annealed results show a similar jump at N ) 48. The transition points correlate with the layer-spacing transitions in Figure 2. The results indicate that there are two stable branches of the immersion energy in the one-layer hydrate region. Hydration and dehydration processes therefore occur discontinuously, as in a phase transition, through jumps between the two branches. The observed dependence upon annealing, reminiscent of hysteresis in experimental systems, results from
J. Phys. Chem. B, Vol. 104, No. 39, 2000 9167
Figure 4. Pressure versus layer spacing isotherms for Sr-montmorillonite at various water contents. The dashed line corresponds to an ambient pressure of 1 bar.
Figure 5. Gibbs energy of swelling for Sr-montmorillonite at various water contents. Gibbs energy units correspond to kilojoules per mole of interlayer ions. The absolute value of the Gibbs energy is set arbitrarily to zero at the two-layer-spacing minimum for each curve.
trapping in a metastable portion of one or both branches around the equilibrium transition point. 3.3. Isotherms. The nature of the two stable branches of the immersion energy and swelling curves was investigated using pressure versus layer-spacing isotherms. Results for Srmontmorillonite with water contents of N ) 44, 48, and 56 are displayed in Figure 4. These water contents were chosen to bracket the observed transitions for the annealed curves shown in Figures 2 and 3. Mechanically stable states for these systems are identified by points where the isotherm has a negative slope and a pressure equal to the externally applied pressure (i.e., crossing points with the external pressure line shown in the figure). All three isotherms show oscillations that lead to multiple stable states for a range of external pressures. For N ) 44 and 48 and ambient external pressure, the isotherms predict stable states with spacings of approximately 12.5 and 14.5 Å, corresponding closely with the transition spacings in Figure 2. The Gibbs energy of swelling, ∆G(z), was calculated for each of the measured isotherms using eq 6. Results are shown in Figure 5. The stable states discussed above are now evident as Gibbs energy minima. The data predict a discontinuous swelling transition with increasing water content. At the lowest water content, the one-layer hydrate spacing is clearly favored, and the two-layer state is metastable. As the water content is increased to N ) 48, the Gibbs energies of the two minima
9168 J. Phys. Chem. B, Vol. 104, No. 39, 2000
Figure 6. Pressure versus layer-spacing isotherms for Sr-montmorillonite and Na-montmorillonite, both with a water content of N ) 48.
become nearly equal while their positions shift only slightly. A further increase in the water content yields a stable two-layer state and the eventual disappearance of the metastable one-layer state. Spacings intermediate between the one- and two-layer values are never identified with the most stable state, and swelling is thus predicted to proceed through a discontinuous jump between one- and two-layer spacings. Cases et al. analyzed swelling experiments for Sr-montmorillonite and concluded that the two-layer spacing forms once each interlayer Sr2+ is associated with 16 water molecules.16 This is in excellent agreement with the simulation value (N ≈ 48 or 16 molecules for each of three interlayer Sr2+ ions) for the equilibrium crossing point to the two-layer spacing. An N ) 48 isotherm was also calculated for Na-montmorillonite and is compared with the strontium isotherm in Figure 6. While both curves show oscillations, the sodium curve has only a single crossing point with the ambient pressure line. Multiple stable states are therefore only evident in Namontmorillonite at negative pressures. Na-montmorillonite isotherms at N ) 44 and 56 are not shown but are similar. The results thus predict that sodium swelling should be qualitatively different, exhibiting a continuous change in layer spacing as the water content is increased. This is consistent with the observed behavior in the immersion energy and swelling curves. 3.4. Interlayer Structure. The peaks in the calculated isotherms for both Na- and Sr-montmorillonite correlate clearly with the ion coordination number. This is illustrated for the N ) 48 Sr-montmorillonite isotherm in Figure 7. The ion coordination number was determined by measuring the average number of water molecules within 3.65 Å of Sr2+, a distance corresponding to the first minimum in the Sr-O radial distribution function in simulations of a dilute SrCl2 solution. (Cl- ions were included in the dilute solution to assure electroneutrality and did not occupy the first hydration sphere of the Sr2+ ion.) The peaks in the isotherm at 13.3 and 14.0 Å correspond to increases in the strontium coordination number to 7.0 and 8.0, respectively. Similar correlation between isotherm peaks and ion coordination number occurs for all water contents in both Na- and Sr-montmorillonite. For example, the N)56 Srmontmorillonite isotherm has peaks at 13.0 and 14.2 Å that correspond to the formation of seven- and eight-coordinate strontium ions, respectively. Average sodium coordination numbers are around 4.0 for spacings inside 13.5 Å and rise to a plateau of approximately 5.5 at a spacing of 14.5 Å. Registry motion in Sr-montmorillonite was severely restricted at low water contents relative to the other clays studied.
Young and Smith
Figure 7. Average coordination number (top graph) and pressure versus layer-spacing isotherm (bottom graph) for Sr-montmorillonite with a water content of N ) 48.
Figure 8. Registry map for Sr-montmorillonite at N ) 48 and z ) 12.7 Å. A snapshot of one interlayer ion with its first hydration shell is also displayed. Each point corresponds to a single registry configuration during the course of the simulation.
Figure 9. Registry map and snapshot as in Figure 8 except that z ) 14.3 Å.
This suggests that the hydrated strontium ions have a stronger interaction with the clay surfaces than the other ions. Srmontmorillonite trajectories were used to produce registry maps in which each registry position during the course of the simulation is represented as a point on a diagram such as that shown in Figure 1. Results for the stable one- and two-layer spacings (12.7 and 14.3 Å) are displayed in Figures 8 and 9, respectively, along with representative snapshots of the Sr2+ hydration structure. At one-layer spacings, the six-fold coordination structure around strontium typically includes four molecules
Simulations of Clay Mineral Swelling and Hydration
Figure 10. Snapshot of the Sr-montmorillonite interlayer region (N ) 48, z ) 14.3 Å) looking down along the z-axis. Four periodic images of the simulation cell are shown with the view taken along the z-axis. The clay lattice is depicted in white and the interlayer solution in black with Sr2+ ions shown as large spheres.
located in a plane horizontal to the clay layers, while the two remaining molecules extend into opposing hexagonal cavities of the siloxane surfaces. The structure requires that the clay maintain a registry configuration with approximately overlapping cavities, and the registry motions are restricted accordingly. For two-layer spacings, the eight-fold coordination around strontium usually consists of a twisted cubic structure with four molecules each above and below the strontium ion. The favored registry position corresponds to one in which the hexagonal cavity of one siloxane sheet overlaps with an oxygen atom on the opposing sheet. The registry motion here is less restricted than in the one-layer state. The water molecules in the eight-fold coordination structure are typically oriented such that one of the water hydrogens is directed toward the clay surface. We therefore speculate that the restricted registry motion enhances the energy of hydrogen-bonding interactions between these hydrogen atoms and the siloxane oxygen atoms without compromising the optimal ion-water coordination structure. A significant fraction of the interlayer volume is unoccupied in low water content, two-layer-spacing states for Sr-montmorillonite. A snapshot of the interlayer structure for N ) 48 and a spacing of 14.3 Å is displayed in Figure 10. For clarity, four periodic images are shown. Half of the interlayer water molecules at this water content are directly coordinated to strontium ions. The remaining water molecules are clustered loosely around the ions, usually as part of a second coordination shell. Unoccupied space is also clustered together and is associated with regions in which the clay has neutral charge, away from substitution sites in the the octahedral sheet. Naand Cs-montmorillonites with N ) 48 have equilibrium layer spacings that are smaller than 14.3 Å (Figure 2). Snapshots of these clays (not shown) reveal some waters coordinated to the interlayer ions with the rest of the molecules distributed evenly throughout the interlayer region. 3.5. Swelling and Hydration Mechanisms. The simulation results relate directly to the issue of whether clay swelling and interlayer hydration occur in a concerted fashion or as a twostep process. A necessary condition for the two-step mechanism is the existence of mechanically stable, partially filled expanded states in which the clay is propped open by hydrated cations. The calculated pressure-spacing isotherms indicate that these states exist, for water contents between one- and two-layer
J. Phys. Chem. B, Vol. 104, No. 39, 2000 9169 hydrates, but only for Sr-montmorillonite among the clays studied. The simulations therefore indicate that the hydration mechanism is dependent upon the hydration energy of the interlayer ion. Unfilled expanded states are unfavorable with respect to clay-clay van der Waals interactions. This energy penalty must be offset by the favorable energy of forming one or two hydration spheres for the ion. Since the hydration energy differs dramatically among ions, there should be a crossover between unstable and stable expanded states, and thus between the concerted and two-step mechanisms, as the interlayer ion hydration energy is increased. For the simulation system, the crossover between the concerted and two-step mechanisms occurs at an ion hydration energy somewhere between that for Na+ and Sr2+. This particular conclusion should be viewed cautiously as it is likely sensitive to details of the simulation model. For example, if the model were to overestimate the strength of the water-clay interactions, it would predict too great a stability for the the collapsed state and would thus underestimate the contribution of expanded states to swelling processes. Quantitative validation of clay simulation models is difficult given the inhomogeneity and disordered structure of experimental clay systems, interference from external surface effects, and the macroscopic nature of most experimental measurements. The existence of metastable states at a single interlayer water content suggests that one possible mechanism for layer-spacing hysteresis involves trapping in these metastable states. The current results suggest that trapping is most likely to occur upon desorption since the two-layer state maintains metastability to quite low water contents. The formation of expanded states is not, however, a necessary condition for layer-spacing hysteresis. For example, a hysteresis mechanism that involves coincidental changes in layer spacing and water content may be inferred from a Cs-montmorillonite swelling free energy reported recently.28 The loop structure of the pressure-spacing isotherms is reminiscent of similar loops associated with gas-liquid coexistence in, for example, a van der Waals fluid. For fluids, raising the temperature along the gas-liquid coexistence line leads eventually to a critical point and associated critical phenomena. It would be interesting to evaluate swelling clay systems to see whether there is evidence for similar critical behavior such as a sharp rise in the system compressibility. 4. Conclusions Molecular dynamics computer simulations have been used to characterize the ion size and charge dependence of swelling and hydration in montmorillonite clay minerals. The picture that has emerged is consistent with the discrete nature of crystalline swelling. For all clays studied, layer spacings measured as a function of water content exhibited plateaus at the one-layer hydrate spacing and calculated immersion energies showed minima for integer-layer hydrates up through the three-layer hydrate. Apparent global minima were identified with the onelayer hydrate for Cs-montmorillonite and with the two-layer hydrates for the Na- and Sr-montmorillonites, in qualitative agreement with experimental observations. The specific mechanism of swelling and hydration was found to be a function of the interlayer ion charge and presumably is related to the ion hydration energy. For Sr-montmorillonite only, layer spacings and immersion energies both jumped discontinuously between one-layer and two-layer values, with the specific transition points depending upon how the system was prepared. Pressure-spacing isotherms calculated for both Sr- and Na-montmorillonite showed a loop structure with
9170 J. Phys. Chem. B, Vol. 104, No. 39, 2000 peaks that clearly correlated with the coordination number of the interlayer ion. These peaks lead to stable expanded states at ambient pressure for the Sr-clay only. The results indicate that swelling and hydration in Sr-montmorillonite is a twostep process involving first expansion, driven by formation of a first and partial second hydration shell for Sr2+, followed by filling of the remaining interlayer volume. For monovalent ions studied, swelling and hydration were coincidental with water occupying all available interlayer volume. In clays that follow the two-step swelling and hydration mechanism, layer-spacing hysteresis may occur through trapping in metastable states. Many issues related to clay swelling and hydration remain unresolved. Several of these cannot be addressed using constantwater-content simulations such as those described in this article. For example, immersion energies provide insight into the energetic origins of swelling behavior but neglect entropic contributions. In addition, the environmental behavior of clay minerals and most experimental measurements of clay behavior are controlled by the water chemical potential and thus occur with variable water content. A grand canonical ensemble simulation method developed recently28 is therefore currently being used to explore additional aspects of clay swelling and hydration behavior. 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. References and Notes (1) Coughtrey, P. J.; Thorne, M. C. Radionuclide Distribution and Transport in Terrestrial and Aquatic Ecosystems; A.A. Balkema: Rotterdam, 1983; Vol. 1. (2) Staunton, S. Eur. J. Soil Sci. 1994, 45, 409. (3) Higgo, J. J. W. Prog. Nucl. Energy 1987, 19, 173. (4) Pusch, R. Clay Mineral. 1992, 27, 353. (5) Nelson, J. D.; Miller, D. J. ExpansiVe SoilssProblems and Practice in Foundation and PaVement Engineering; John Wiley & Sons, Inc.: New York, 1992. (6) Durand, C.; Forsans, T.; Ruffet, C.; Onaisi, A.; Audibert, A. ReV. Inst. Fr. Pet. 1995, 50, 187. (7) Norrish, K. Discuss. Faraday Soc. 1954, 18, 120. (8) Viani, B. E.; Low, P. F.; Roth, C. B. J. Colloid. Interface Sci. 1983, 96, 229. (9) Slade, P. G.; Quirk, J. P. J. Colloid Interface Sci. 1991, 144, 18. (10) Sato, T.; Watanabe, T.; Otsuka, R. Clays Clay Miner. 1992, 40, 103. (11) Zhang, F.; Zhang, Z. Z.; Low, P. F.; Roth, C. B. Clay Miner. 1993, 28, 25.
Young and Smith (12) Laird, D. A.; Shang, C.; Thompson, M. L. J. Colloid Interface Sci. 1995, 171, 240. (13) Mooney, R. W.; Keenan, A. G.; Wood, L. A. J. Am. Chem. Soc. 1952, 74, 1371. (14) Calvet, R. Ann. Agron. 1973, 24, 77. (15) Be´rend, I.; Cases, J. M.; Francois, M.; Uriot, J. P.; Michot, L.; Masion, A.; Thomas, F. Clays Clay Miner. 1995, 43, 324. (16) Cases, J. M.; Be´rend, I.; Francois, M.; Uriot, J. P.; Michot, L. J.; Thomas, F. Clays Clay Miner. 1997, 45, 8. (17) Chiou, C. T.; Rutherford, D. W. Clays Clay Miner. 1997, 45, 867. (18) Cancela, G. D.; Huertas, F. J.; Taboada, E. R.; Sa´nchez-Rasero, F.; Laguna, A. H. J. Colloid Interface Sci. 1997, 185, 343. (19) Fu, M. H.; Zhang, Z. Z.; Low, P. F. Clays Clay Miner. 1990, 38, 485. (20) Sposito, G.; Prost, R. Chem. ReV. 1982, 82, 553. (21) Everett, D. H. Adsorption Hysteresis. In The Solid-Gas Interface; Flood, E. A., Ed.; Marcel Dekker: New York, 1967; Vol. 2. (22) Laird, D. A. Clays Clay Miner. 1999, 47, 630. (23) Skipper, N. T.; Refson, K.; McConnell, J. D. C. J. Chem. Phys. 1991, 94, 7434. (24) Delville, A. Langmuir 1991, 7, 547. (25) Skipper, N. T.; Chang, F.-R. C.; Sposito, G. Clays Clay Miner. 1995, 43, 285. (26) Boek, E. S.; Coveney, P. V.; Skipper, N. T. J. Am. Chem. Soc. 1995, 117, 12608. (27) Teppen, B. J.; Rasmussen, K.; Bertsch, P. M.; Miller, D. M.; Scha¨fer, L. J. Phys. Chem. B. 1997, 101, 1579. (28) Shroll, R. M.; Smith, D. E. J. Chem. Phys. 1999, 111, 9025. (29) Smith, D. E. Langmuir 1998, 14, 5959. (30) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (31) Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 100, 3757. (32) Smith, D. E.; Dang, L. X. Chem. Phys. Lett. 1994, 230, 209. (33) Friedman, H. L.; Krischnan, C. V. In Water, A ComprehensiVe Treatise; Franks, F., Ed.; Plenum: New York, 1973; Vol. 6. (34) Albright, J. N. J. Chem. Phys. 1972, 56, 3783. (35) Sposito, G. The Chemistry of Soils; Oxford University Press: New York, 1989. (36) Skipper, N. T.; Refson, K.; McConnell, J. D. C. Clay Miner. 1989, 24, 411. (37) Smith, D. E. Unpublished. (38) Smith, D. E.; Haymet, A. D. J. J. Chem. Phys. 1992, 96, 8450. (39) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. J. Chem. Phys. 1982, 76, 637. (40) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, 1987. (41) de Leeuw, S. W.; Perram, J. W.; Smith, E. R. Proc. R. Soc. London, A 1980, 373, 27. (42) Andersen, H. C. J. Comput. Phys. 1983, 52, 24. (43) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (44) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (45) Delville, A.; Sokolowski, S. J. Phys. Chem. 1993, 97, 6261. (46) Chang, F.-R. C.; Skipper, N. T.; Sposito, G. Langmuir 1995, 11, 2734.
