Effect of layer bending on montmorillonite hydration and structure from molecular simulation

Abstract

Conceptual models of smectite hydration include planar (flat) clay layers that undergo stepwise expansion as successive monolayers of water molecules fill the interlayer regions. However, X-ray diffraction (XRD) studies indicate the presence of interstratified hydration states, suggesting non-uniform interlayer hydration in smectites. Additionally, recent theoretical studies have shown that clay layers can adopt bent configurations over nanometer-scale lateral dimensions with minimal effect on mechanical properties. Therefore, in this study we used molecular simulations to evaluate structural properties and water adsorption isotherms for montmorillonite models composed of bent clay layers in mixed hydration states. Results are compared with models consisting of planar clay layers with interstratified hydration states (e.g. 1W–2W). The small degree of bending in these models (up to 1.5 Å of vertical displacement over a 1.3 nm lateral dimension) had little or no effect on bond lengths and angle distributions within the clay layers. Except for models that included dry states, porosities and simulated water adsorption isotherms were nearly identical for bent or flat clay layers with the same averaged layer spacing. Similar agreement was seen with Na- and Ca-exchanged clays. While the small bent models did not retain their configurations during unconstrained molecular dynamics simulation with flexible clay layers, we show that bent structures are stable at much larger length scales by simulating a 41.6×7.1 nm² system that included dehydrated and hydrated regions in the same interlayer.

Introduction

Bentonite and argillaceous clay-bearing rocks are considered in the design concepts of several deep geological repositories for the disposal of high level nuclear waste packages as an engineered barrier (Sellin and Leupin, 2013) and as the surrounding host rock (Grambow, 2016), respectively. Varying humidity and temperature conditions due to chemical or geologic processes in the near-field could affect the structural properties of swelling clay phases in repository environments. Specifically, smectite clays can undergo intracrystalline expansion or contraction in response to changes in relative humidity (RH). These hydration/dehydration phenomena can lead to changes in clay swelling pressures exerting undue forces on surrounding barrier materials (e.g., waste canister surface) or induce shrinkage crack formation in bentonite material (Soe et al., 2009). This can result in changes in the bulk porosity and permeability of engineered barriers. Therefore, it is important to fully understand the hydration and dehydration processes of expansive clays at the molecular level.

Smectite clays such as montmorillonite (Mnt) are thought to exist in stable hydration states characterized by the number of water layers in between the aluminosilicate layers. As RH increases, smectites undergo stepwise intracrystalline swelling as the d-spacing abruptly increases with increasing number of water layers from zero (0W), one (1W), two (2W), and three (3W) in the interlayers (Bradley et al., 1937). More recent X-ray diffraction (XRD) studies have revealed some intermediate peaks rather than distinct peaks for 1W, 2W, or 3W hydration states (Moore and Hower, 1986; Yamada et al., 1994; Tamura et al., 2000). These intermediate peaks have been interpreted as a result of the random interstratification of two or more distinct hydration states (Sato et al., 1992; Sato et al., 1996; Tamura et al., 2000) and are included in XRD profile modeling (Ferrage et al., 2010; Ferrage, 2016; Chaaya et al., 2023).

A recent molecular dynamics (MD) simulation and experimental study of smectite clay dehydration revealed a two-stage dehydration process (Ho et al., 2022). Based on calorimetric data, the first stage includes a rapid mass loss and partial reduction in crystallographic d spacing. Accompanying MD simulations indicated that water molecules not coordinated to interlayer cations are lost during this first stage. In the second stage, both the rate of mass loss and reduction in d spacing are reduced as water molecules coordinating interlayer cations are removed. Variations in the strength of ion–water and ion–clay interactions have been linked to the hysteresis observed in smectite hydration–dehydration experiments (Fu et al., 1990; Laird et al., 1995), as an explanation for thermo-mechanical irreversibility in smectite swelling behavior (Boek and Coveney, 1995; Honorio et al., 2017; Brochard, 2021). Irregular water loss during smectite dehydration can also result in localized bending of clay layers (Honorio et al., 2018), particularly near particle edges (Ho et al., 2019).

Recent theoretical studies have included clay layer bending related to osmotic swelling and tactoid formation (Whittaker et al., 2023), as well as energetic and mechanical aspects of clay layer bending during transitions between thermodynamically stable hydration states (Honorio et al., 2017; Honorio et al., 2018). For example, the transitions from 0W to 1W and from 1W to 2W of Na-Mnt at a confining pressure of 10 MPa and a temperature of 300 K encounter energy barriers of 7 kT nm^–2 and 2 kT nm^–2, respectively, and the energy barriers decrease with increasing temperature (Honorio et al., 2017). These values are consistent with separate simulation studies at room temperature that revealed energy barriers between hydration states in Na-Mnt of ~8−17 kT nm^–2 for the 0W–1W transition, decreasing to ~2 kT nm^–2 for the 1W–2W transition and <1 kT nm^–2 for the 2W–3W transition (Brochard, 2021; Ho et al., 2019; Shen and Bourg, 2021). In the limit of large lateral dimensions for clay particles, Honorio et al. (2018) combined plate theory with molecular simulations to show that such hydration energy barriers are overcome by bending of clay layers (i.e. atomic displacement in c) over a relatively small lateral dimension (a or b) of ~1 nm. Thus it is expected that for clay particles with large aspect ratios, the hydration transition involves bending of clay layers in localized areas of the particle, rather than an instantaneous translational displacement of an entire layer.

The vast majority of molecular simulation studies of smectite clays use planar (flat) models of clay layers based on time-averaged structures obtained from diffraction studies of related minerals (Catti et al., 1994; Drits et al., 2012). Recently, a comprehensive grand canonical Monte Carlo (GCMC) simulation study of water adsorption in Na-Mnt predicted stable and metastable hydration states in good agreement with experiments (Brochard, 2021). However, only standard models consisting of planar (flat) clay layers in stable hydration states (1W–3W) were used. To our knowledge, the hydration behavior of interlayer regions associated with bent clay layers has not been reported previously. Therefore, in this work we compare simulated water adsorption isotherms in Mnt models composed of bent clay layers in mixed hydration states with those consisting of planar clay layers and interstratified hydration states (i.e. 0W–1W). The effect of interlayer cations is included by comparing results for Na- and Ca-exchanged Mnt. Results are also shown for a large-scale MD simulation of a more realistic clay stack that includes dehydrated and hydrated regions in the same interlayer.

Materials and methods

Water adsorption isotherms at 300 K were obtained from GCMC simulations using the Towhee code (Martin, 2013) in the grand canonical ensemble. The Clayff parameter set (Cygan et al., 2004) together with the rigid extended simple point charge (SPC/E) water model (Berendsen et al., 1987) was used for all Coulombic and van der Waals interactions. Clayff has been used extensively for molecular simulation of clay minerals (Cygan et al., 2021), including water adsorption in Mnt (Brochard, 2021; Teich-McGoldrick et al., 2015; Honorio et al., 2017). Modified interaction parameters between interlayer Na⁺ ions and surface O atoms (ob atom type) were used as these modifications resulted in a more realistic description of energy barriers between hydration states (Ho et al., 2019). Three-dimensional periodic boundary conditions (3D PBCs) were used with a short-range cut-off distance of 10.0 Å, and long-range electrostatic interactions were evaluated using Ewald summation with a precision of 1×10⁻⁴.

The Mnt supercell (Fig. 1) consisted of two clay layers each containing 12 unit cells (4×3×2 repeats in a, b, and c) generated from the pyrophyllite structure (Lee and Guggenheim, 1981). After orthogonalization, the lateral (xy) dimensions of the supercell were 20.8 Å×26.9 Å. A layer charge of −0.83 e/u.c. (elementary charge per unit cell) was introduced in each octahedral sheet by replacing 10 randomly chosen Al atoms with Mg atoms while avoiding Mg-O-Mg close contacts. This layer charge is slightly larger than typically used for Mnt simulations (−0.75 e/u.c.) due to the slightly larger supercell (3b). The layer spacing was expanded to specified values (discussed below), and counterions (Na⁺ or Ca²⁺) were randomly inserted in each midplane region for charge balance. The Mnt u.c. formula for Na⁺ was Na_0.833[Si₈](Mg_0.833Al_3.167)O₂₀(OH)₄.

Figure 1. Initial 0W–1W models (yz plane) for GCMC simulations of Na-Mnt in the flat (A) and bent (B) configurations. Dimensions are 20.8×26.9×22.0 Å³. Atoms are colored as follows: red (O), white (H), yellow (Si), magenta (Al), green (Mg), and purple (Na). Periodic boundaries are shown in y as black lines (simulations were performed using 3D PBCs). Layer spacings are indicated with arrows.

Adjustments were made to the cell z-parameter and z-coordinate of clay layers so that models corresponded to mixed hydration states (0W–1W, 1W–2W, 2W–3W) with the clay layers in either planar (flat) or bent configurations. Except for the 0W layer spacings, which were determined from separate simulations to be 9.4 Å and 10.0 Å for Ca-Mnt and Na-Mnt, respectively, uniform layer spacing values of 12.0 Å (1W), 15.0 Å (2W), and 18.0 Å (3W) were used to represent hydration states regardless of interlayer cation. The z-dimension of each model system is the sum of the layer spacings of each single hydration state.

Models with flat clay layers were generated by shifting the layers in z with layer spacings corresponding to each stable hydration state (e.g. 1W and 2W). Models with bent clay layers were generated as follows: starting from a layer spacing midway between the stable hydration states (e.g. 13.5 Å for 1W–2W), the z-coordinate of all layer atoms were shifted up or down in z as a function of the atom’s y-coordinate so that the layer spacing varied continuously between the stable hydration states. Using the 0W–1W hydration state as an example (Fig. 1), the layer spacing of the bent configuration varied from 9.4 Å at y=0 to 12.0 Å at the midpoint (y=13.45 Å). The other layer was bent in the opposite direction. In this study bending was arbitrarily chosen to occur only in the y direction, although bending in x is plausible based on similar elastic moduli in lateral dimensions of clay particles (Teich-McGoldrick et al., 2012; Zartman et al., 2010). This type of bent geometry is consistent with a recent theoretical analysis of clay layer bending (Honorio et al., 2018) which showed that localized bending between hydration states can occur over bending lengths as small as 1 nm. Different degrees of bending would be expected in larger models (see below), but this bending scheme representing two hydration states over a nanometer-scale lateral distance is consistent with conceptual models in the literature (Ferrage et al., 2007). Layer spacings and porosities of all models are given in Table 1.

Table 1. Layer spacing and porosity for GCMC models

^a For bent models, values represent the average as layer spacing varies with y (as described in the text). For flat models, layer spacings for smaller and larger interlayers are shown.

^b Ratio of free volume to total volume using the Connolly surface method (Connolly, 1983) in BIOVIA Materials Studio (Dassault Systèmes SE) with a probe diameter (3.2 Å) corresponding to a water molecule.

GCMC moves and probabilities included configurational-bias (Siepmann and Frenkel, 1992) water insertion (20%) and deletion (20%), configurational-bias intrabox water transfer (15%), configurational-bias water molecule regrowth (15%), translation of interlayer water molecules and cations (20%), and rotation of water molecules (10%). No moves involved atoms in the clay lattice as these atoms were held fixed throughout the simulations. Each system was simulated for 2.5×10⁸ steps, with statistics for water loading computed over the final 1.0×10⁸ steps. Values of the water chemical potential (μ) corresponding to RH from 0.5 to 93% were based on previous GCMC simulations of pure SPC/E water (Brochard, 2021).

Large-scale MD simulations at 300 K were performed on Na-Mnt using 3D PBCs in both flat and bent configurations. The simulation was performed with the LAMMPS code (Thompson et al., 2022) using Clayff parameters with the flexible SPC water model (Teleman et al., 1987). Temperature and pressure were controlled with a Nose–Hoover thermostat and barostat with relaxation times of 100 and 1000 fs, respectively. A timestep of 1 fs was used with a 10.0 Å cut-off distance for short-range interactions and a particle-mesh Ewald algorithm (Plimpton et al., 1997) with a precision of 1×10⁻⁴ for long-range Coulomb interactions. To construct the simulation cell, an orthogonalized Mnt layer with xy dimension 4.16×3.58 nm² was created in a 4.16×3.58×1.8 nm³ simulation box (8×4×1 u.c. repeats in a, b, and c) as shown in Fig. 2A. The layer charge of the Mnt layer was −0.75 e/u.c. This structure was then replicated 10×2×20 times to create a large clay structure containing 521,600 atoms (Fig. 2B). Water molecules were then randomly inserted in each interlayer, resulting in a water content of 2.34 H₂O/u.c. (1W). After minimization, the system was equilibrated for 50 ps at constant volume followed by 1 ns at a constant pressure (P_zz) of 1 atm to equilibrate the layer spacing. Here, P_zz refers to pressure (stress) applied normal to the basal surfaces, so that volume changes only involved changes in the supercell c parameter. The final supercell with dimensions 41.6×71.6×23.8 nm³ is shown in Fig. 2C. The average layer spacing was ~11.9 Å, confirming the monolayer (1W) hydration state. Water uniformly distributes in each interlayer and each clay layer is flat.

Figure 2. Model preparation for large-scale MD simulation of Na-Mnt with heterogeneous interlayers. (A) Initial simulation box with one clay layer (gray) and Na⁺ ions (yellow). (B) Large clay structure by replicating clay structure in A. (C) Homogeneous 1W configuration (water O red) with dimensions ~ 41.6×7.1×23.3 nm³.

To create a bent clay configuration, the simulation box in Fig. 2A was filled randomly with water molecules (i.e. 2.36 H₂O/u.c.) and five Lennard-Jones (LJ) particles ( $ \unicode{x025B} $ =0.294076 kcal mol^–1, $ \unicode{x03C3} $ =8.1656 Å). The large $ \unicode{x03C3} $ value is required to locally expand the interlayer to accommodate the 2W state. This system was then replicated 10×2×20 times to create a large structure similar to that shown in Fig. 2C. The large structure was then equilibrated for 10 ps at constant volume, followed by 500 ps at constant pressure (P_zz=100 atm), and finally 2 ns at constant pressure (P_zz=1 atm). As seen in Fig. 3, the LJ particles aggregated in the interlayer, acting as a spacer that caused localized bending of the clay layers. Note that interactions between water molecules and LJ particles were turned off, so water molecules can coexist with the LJ particle aggregates. The LJ particles were then removed, and the system was equilibrated during a 100 ps constant pressure simulation (P_zz=100 atm) followed by 3 ns P_zz=1 atm. The final box dimensions are 41.6×71.6×23.4 nm³.

Figure 3. Equilibrium snapshot from a large-scale MD simulation of Na-Mnt showing interlayers containing spacer particles (pink spheres) to create the bent structure.

XRD patterns were computed from the MD-equilibrated large-scale models in both flat and bent configurations. The positions of atoms obtained from MD simulation were used to compute the XRD intensity I (Coleman et al., 2013):

(1) $$ I= Lp\left(\unicode{x03B8} \right)\frac{F\ast F}{N}, $$

where $ Lp\left(\unicode{x03B8} \right) $ is the Lorentz-polarization factor, $ \unicode{x03B8} $ is the scattering angle of diffraction, N is the number of atoms. The structural factor F was calculated as:

(2) $$ F(k)=\sum \limits_{\mathrm{j}=1}^N{f}_{\mathrm{j}}\left(\theta \right)\exp \left(2\pi ik{r}_{\mathrm{j}}\right) $$

where k is the location of the reciprocal lattice node, f _j is the atomic scattering factor (Peng et al., 1996), and r _j is the position of each atom obtained from the MD simulation. The XRD reflection was calculated with the spacing of the reciprocal lattice node of 0.1 Å in the z direction and the wavelength is 1.541838 Å. The XRD profiles were recorded from 3 to 33° 2 $ \unicode{x03B8} $ with a resolution of 0.1° 2 $ \unicode{x03B8} $ . Note that the XRD calculation was computationally expensive.

Results and Discussion

The effect of clay-layer bending on interlayer hydration was evaluated from GCMC simulations. Water adsorption isotherms for Na-Mnt and Ca-Mnt in mixed hydration states are shown in Fig. 4. As expected, water loading increases with accessible volume (i.e. 0W–1W<1W–2W<2W–3W; see Table 1 for porosity). Curvature of the isotherm profiles as a function of RH is consistent with previous GCMC simulations of water adsorption in stable hydration states (i.e. 1W, 2W, 3W) (Hensen et al., 2001; Smith et al., 2006; Ferrage et al., 2011; Dazas et al., 2015; Teich-McGoldrick et al., 2015). Specifically, the RH required to achieve maximum water loading increases with accessible interlayer volume. Models with the smallest accessible volume (0W–1W) achieve maximum loading at very low RH, while models with the largest accessible volume (2W–3W) do not achieve maximum water loading until ~60% RH. However, isotherms for the Ca-Mnt models for these states show a more gradual increase in water loading at low RH, because there are fewer interlayer cations that serve as initial adsorption sites.

Figure 4. Water adsorption isotherms from GCMC simulations at 300 K for mixed hydration states of (A) Na-Mnt and (B) Ca-Mnt in bent and flat configurations.

Models with bent clay layers (representative of transition hydration states) show similar or identical water adsorption behavior compared with the standard models (flat clay layers with single hydration states). However, bent models with the smallest layer spacing (0W–1W) show slightly lower water loading than the corresponding flat model at all RH values. As seen in Table 1, the 0W–1W bent models have significantly lower porosities than the flat models. For the 1W–2W and 2W–3W systems, the bent and flat models have similar porosities and water loadings at all RH.

GCMC snapshots for the 0W–1W models (Fig. 5) illustrate the effect of reduced porosity on water adsorption. A large proportion of the bent interlayer region has little or no adsorbed water, indicating that there is not enough accessible volume for a water layer to form. Additional snapshots from all GCMC simulations are shown in the Supplementary material (Figs S1–S3). All of the interlayer regions are accessible to water in the expanded states (1W–2W and 2W–3W), so clay-layer bending has little effect on water loading at higher RH when these hydration states are thermodynamically stable (Brochard, 2021).

Figure 5. Equilibrium 0W–1W snapshots (yz plane) of Na-Mnt in the flat (A) and bent (B) configurations. Atoms are colored as follows: red (O), white (H), yellow (Si), magenta (Al), green (Mg), and purple (Na).

As further evidence of the minimal impact of clay-layer bending on montmorillonite structure, a comparison of lattice bond lengths and angles (see Supplementary material, Figs. S4 and S5) shows nearly identical distributions for bent and flat models. The small degree of layer bending in this study results in a small vertical (z) displacement of 1.0–1.5 Å over a 1.3 nm lateral dimension, which has no discernable effect on the bond or angle distributions. Of course, more significant bending of the layers would cause more measurable increases in bonds and angles, resulting in lattice strain. In addition, short 50 ps MD simulations with identical water loadings in each bent/flat pair (see Supplementary material, Fig. S6) show that layer bending has no effect on potential energy in the 1W–2W and 2W–3W models and only increases slightly in the potential energy in the 0W–1W models. Note that these potential energy comparisons are limited to the interlayer species as the clay layers are held fixed.

This depiction of variable layer spacing within a single interlayer is consistent with a recent study of bending moduli for clay layers (Honorio et al., 2018), which showed that isolated clay layers are quite flexible even during stirring in an exfoliated state, with bending moduli in the lateral dimensions ~2×10⁻¹⁷ N m. For stacks of clay particles with large aspect ratios, the authors combined plate theory with molecular simulations to show that hydration energy barriers (<10 kT nm^–2) are overcome by bending of clay layers over a relatively small bending length of ~1 nm (Honorio et al., 2018).

The bent GCMC models reformed their original flat configurations during subsequent constant-P_zz MD simulations (1 atm), resulting in homogeneous interlayers with layer spacings intermediate between stable hydration states. However, the large bent model (41.6×71.6×23.8 nm³) maintained its configuration throughout the final 3 ns simulation (Fig. 6). Removal of the LJ spacer particles in Fig. 3 resulted in some local rearrangement of clay layers and interlayer water, but dry interlayer regions are clearly visible. Interlayer regions corresponding to 0W, 1W, and even 2W can exist within a single interlayer of the structure. Note that the average water content in each interlayer was the same (i.e. 2.36 H₂O/u.c.). Results from the large-scale MD simulation indicate that (1) the persistence of bent clay layers depends on simulation system size, and (2) localized layer bending in larger clay stacks appears to be long lived. Similar layer bending was also seen during large-scale simulations of organically modified smectites with interlayer ethanol (Metz et al., 2015). Additionally, thermal undulations in layered double hydroxides seen in large-scale MD simulations enabled the calculation of bending moduli (Thyveetil et al., 2007).

Figure 6. Equilibrium snapshot from a large-scale MD simulation of Na-Mnt showing heterogeneous 0W–1W–2W interlayers after removal of spacer particles.

In Fig. 7 we compare the XRD profile for the bent (heterogeneous) model (Fig. 6) with that of the flat homogeneous model (Fig. 2C; 2.34 H₂O/u.c.). The results indicate that heterogeneity in the same interlayer at the length scale studied here does not result in any significant difference in the (001) XRD reflection, but such heterogeneity does result in a broadening of the peaks at larger values of 2θ. This broadening has been attributed to crystal strain due to variability in atomic z coordinates in XRD samples (i.e. layer bending), and can be quantified with a variance parameter s_z (Ferrage et al., 2005). Importantly, larger values of s_z are associated with incomplete hydration transitions, which is accommodated by layer flexibility (Ferrage et al., 2005).

Figure 7. Comparison of XRD profiles from large-scale MD simulations of Na-Mnt with heterogeneous (black) and homogeneous (red) interlayer hydration.

It is possible that water distribution in the same interlayer is not uniform (Wang and Xu, 2006). A conceptual model of smectite dehydration based on kinetic analysis of XRD patterns at different temperatures includes localized layer bending as water molecules diffuse out of the interlayers (Ferrage et al., 2007). Similarly, analysis of X-ray scattering and cryogenic microscopy indicates a two-step dehydration process involving diffusional mixing and localized layer bending (Whittaker et al., 2019). As discussed above, XRD patterns of heterogeneous hydration structures are characterized by irrational, asymmetrical, and broad reflections (Ferrage et al., 2005). The issue of heterogeneous hydration states in smectites becomes more relevant when considering continuous hydration processes (i.e. not stepwise swelling mechanism) of natural smectite due to heterogeneous charge distribution on the clay layers (Chipera et al., 1995; Tamura et al., 2000; Morodome and Kawamura, 2011; Schleicher et al., 2013).

Conclusions

We used molecular simulations to compare water adsorption isotherms for Mnt models corresponding to transitions between standard hydration states. Hydration transitions may cause bending and undue forces in the near field of clay phases, so it is important to quantify such effects in relevant scenarios such as underground repository settings.

GCMC simulations on small model systems (nanometer dimensions) were performed to compare the effect of layer bending on water adsorption. For models with expanded interlayers (1W–2W, 2W–3W), porosities and simulated water adsorption isotherms were nearly identical between models with bent clay layers representing mixed hydration states and those with flat clay layers with (interstratified) standard hydration states. Bent models involving dry states (0W–1W) have reduced porosities and water adsorption compared with the corresponding flat models. Similar results were obtained for models with monovalent (Na⁺) and divalent (Ca²⁺) interlayer cations, indicating that ion hydration properties have a minor influence layer bending, within the models studied. However, an analysis of (001) peak irrationality suggests that more fluctuations in layer thickness in the presence of divalent cations may be due to the lower cation density in the interlayer (Ferrage et al., 2005).

Bond lengths and angles in the bent models with a small vertical displacement (up to 1.5 Å) are essentially unchanged compared with the standard flat models. Although the small models used in the GCMC simulations do not retain the bent configuration during unconstrained MD simulation, large-scale simulations (with dimensions similar to actual clay particles) indicate that models with localized bending of clay layers are stable over MD timescales. The tendency for bent layers to flatten clearly depends on system size. Consistent with a previous analysis of mechanical properties of clay-layer bending (Honorio et al., 2018), our results indicate that stacks of Mnt layers are able to exist in mixed hydration states (and localized bending) with minimal effect on macroscopic clay structure or water adsorption properties. Additional research is needed to confirm the coexistence of these configurations under hydrating conditions. Finally, this study has established techniques for creating model systems at various length scales which can be used to explore the effect of layer bending on smectite properties.