The swelling and cation exchange properties of montmorillonite are fundamental in a wide range of applications ranging from nanocomposites to catalytic cracking of hydrocarbons. The swelling results from several factors and, though widely studied, information on the effects of a single factor at a time is lacking. In this study, density functional theory (DFT) calculations were used to obtain atomic-level information on the swelling of montmorillonite. Molecular dynamics (MD) was used to investigate the swelling properties of montmorillonites with different layer charges and interlayer cationic compositions. Molecular dynamics calculations, with CLAYFF force field, consider three layer charges (–1.0, –0.66 and –0.5 e per unit cell) arising from octahedral substitutions and interlayer counterions of Na, K and Ca. The swelling curves obtained showed that smaller layer charge results in greater swelling but the type of the interlayer cation also has an effect. The DFT calculations were also seen to predict larger d values than MD. The formation of 1, 2 and 3 water molecular layers in the interlayer spaces was observed. Finally, the data from MD calculations were used to predict the self-diffusion coefficients of interlayer water and cations in different montmorillonites and in general the coefficient increased with increasing water content and with decreasing layer charge.
Smectites are 2:1 layered swelling clay minerals consisting of layers formed of two tetrahedral silicate sheets one on either side of an octahedral aluminium sheet. The charged layers have a negative charge and are stacked together with an interlayer space containing hydrated charge-compensating cations. Smectites have the ability to adsorb water in the interlayer space and to modify the cationic composition (Brigatti et al., 2006; Velde & Meunier, 2008).
The properties of smectites are under extensive study due to their important role in a wide variety of technological and biochemical applications such as fluid filtration, catalytic processes and the disposal of high-level radioactive waste (Murray, 2000; Dohrmann et al., 2013). In the Finnish concept of geological disposal of nuclear waste, a multi-barrier system is planned to prevent/retard the release of radionuclides to the biosphere. One of these barrier materials is bentonite which consists largely of the smectite-group mineral, montmorillonite. The retardation of radionuclides is based mainly on diffusional transport and sorption on solid mineral surfaces. Due to the presence of smectites, bentonite swells in contact with water thus blocking transport pathways.
Montmorillonite is a fine-grained and mechanically stable mineral characterized by layer charge (e/unit cell) with the distribution of charge sites between the octahedral and tetrahedral sheets. The swelling and cation exchange properties of the mineral also depend on solution properties such as pH, ionic strength and ligands.
Isotopic difference neutron diffraction and molecular modelling techniques have made it possible to study the interlayer structure of clay and the effect of hydration on the clay structure (Sposito et al., 1999). In sorption studies, molecular modelling has typically been used with surface characterization techniques such as X-ray photoelectron spectroscopy (Ebina et al., 1999) and extended X-ray absorption fine structure spectroscopy (Hattori et al., 2009). The relationship between atomistic phenomena and surface complexation modelling has also been formulated by Wesolowski et al. (2009). Quantum mechanics (QM) techniques have been used to study clay–cation and clay–cation–water interactions (Chatterjee et al., 1999; Sposito et al., 1999; Ren et al., 2012; Tribe et al., 2012). A recent multiscale approach including ab initio molecular dynamics (MD), classical MD and mesoscopic simulations emphasized the challenge in describing the phenomena with different length and time scales in clay studies (Rotenberg et al., 2014). The benefits of using molecular modelling are evident because the understanding of the chemical reactions requires that the description of the electron structure of the system and the dependence of the thermodynamics as a function of time are known.
MD is used to study larger and more complex systems as compared to QM, e.g. to evaluate equilibrium and transport properties that cannot be calculated analytically (Tambach et al., 2004; Tao et al., 2010). While these properties are determined in the microscopic i.e. atomistic level, they are linked to the macroscopic properties of the bulk system through thermodynamics and statistical mechanics. In this work, QM and MD simulations are used to study the effect of layer charge on the swelling properties of montmorillonite. Although the swelling of montmorillonite has been studied extensively using molecular modelling, the structures used in those studies vary in terms of the location of substitutions (tetrahedral and/or octahedral), the layer charges and the types of interlayer cations (all of which may affect the swelling properties). While considering the effect of layer charge on the crystalline swelling of montmorillonite, comparison of results from different sources increases the number of variables leading to uncertainty in the comparison. In the present study, the layer charge stems only from octahedral substitutions and the same structures are applied for different interlayer cations, rendering comparison of the results more reliable. In addition, two different molecular modelling methods (QM and MD) are applied to the same structures, making comparison of the methods straightforward. The QM simulations are performed for a tridimensional periodic system, and MD simulations for a rectangular simulation box containing one clay lamella and the interlayer filled with water and cations. The optimized mineral structures defined by the QM methods are utilized in MD simulations. The results from the small-scale system QM calculations are compared with those from the slightly larger-scale system MD calculations. For the MD simulations, the simulation box was large enough to describe specific chemical interactions in the system in order to compare between different types of montmorillonites. Using larger systems, bending of clay lamellae, for example, may prohibit the detection of nanoscale phenomena.
In order to perform reliable MD calculations, accurate crystal structures of pure montmorillonites are required as a starting point for the calculations. In the present study, the optimized unit-cell structures were determined from QM for Na-, K- and Ca-montmorillonites. As a starting point for QM calculations, experimental atomic coordinates defined by Tsipursky and Drits (1984) for K-smectite were used. The calculations were performed with the density functional CASTEP (CAmbridge Serial Total Energy Package, Clark et al., 2005) code implemented in the Materials Studio versions 6.0 (Accelrys, 2011) and 7.0 (Accelrys, 2013).
Optimization of the unit cell of K-montmorillonite was performed using different exchange-correlation energy potentials, pseudopotentials and kinetic cut-off energy values to validate the results compared to experimental values. The exchange-correlation was described with generalized gradient approximation GGA-PBE or GGA-PW91. In these calculations, the total electronic energy and overall electronic density distribution were solved in order to define the energetically stable montmorillonite structures (Leach, 2001). As a compromise between the accuracy and computational time of calculations, the ultrasoft pseudopotentials (Table 1) were used for each element, and the kinetic cut-off energy for a plane-wave expansion of the wave function was 310 eV.
Optimizations of the unit cells of Na- and Ca-montmorillonites were done using GGA-PBE approximation only. The pseudopotentials are presented in Table 1, and the kinetic cut-off energy for a plane-wave expansion of the wave function was 310 eV as in the case of K-montmorillonite.
The optimized unit cells were used to construct model structures (1 × 3 × 1 unit cells) for MD simulations. The model structures were reoptimized using GGA-PBE approximation with the pseudopotentials presented in Table 1. For oxygen and silicon, O_00PBE.usp and Si_00PBE.usp potentials were used. In these calculations, the kinetic cut-off energy for a plane-wave expansion of the wave function was 370 eV.
The 1 × 3 × 1 model structures were utilized to study enlargement of the montmorillonite structure as a function of water insertion into the interlayer space of montmorillonites. Optimization of the model structures with different water contents was performed using the GGA-PW91 approximation. GGA-PW91 gives better description of the system, and calculations achieve convergence faster than using GGA-PBE. Further, the GGA-PBE overestimates the thickness of the interlayer space for larger systems, because it cannot describe the energy states of the system correctly. Calculations were performed using the pseudopotentials presented in Table 1. The kinetic cut-off energy for a plane-wave expansion of the wave function was 370 eV.
Hydration of cations was investigated using the density functional DMol3 code implemented in Materials Studio version 7.0 (Accelrys, 2013). The coordination sphere of cations (Na+, K+ and Ca2+) was optimized using the GGA-PW91 approximation with DNP basis set and DFT-based semi-core pseudopotentials. Hydration energies were calculated as a difference between the total energy of a hydrated cation and the total energies of a non-hydrated cation and water molecules.
All of the MD simulations in the present work were conducted using the LAMMPS (Plimpton, 1995) simulation software. The interatomic interactions for these hydrated montmorillonite systems were simulated using the CLAYFF (Cygan et al., 2004) force field which was chosen because of its previous successful use in montmorillonite mineral simulations (Greathouse & Cygan, 2005; Tao et al., 2010). The CLAYFF model differentiates an element into different types based on its structural position (tetrahedral or octahedral sheet and the neighbouring elements, e.g. substitution sites) and assigns partial charges and van der Waals parameters to each atom type. It defines a bond-stretch parameter for the structural hydroxyl groups and an angle bend parameter for the octahedral metal and OH. The interlayer water parameters in CLAYFF are from the flexible simple point charge (SPC) water model (Berendsen et al., 1981). Simulations were performed in the isothermal-isobaric (NPT) ensemble with three-dimensional periodic boundary conditions. After pre-equilibration the systems were equilibrated for 0.5 ns, followed by the production run for another 0.5 ns at a temperature of 298 K and at pressure of 1 bar . A time step of 1 fs was used in the production run. To avoid translational drift of the clay layer, the centre of mass of the layer was fixed during the simulation. The data were recorded every 200 fs. The long-range electrostatic interactions were calculated using the Ewald summation method (Karasawa & Goddard, 1989) with a precision of 1e–4. The short-range electrostatic and van der Waals interactions were calculated using a 12.5 Å cut-off.
Results and Discussion
The aim was to have comparable calculated and experimental lattice parameters of the unit-cell structures so that the modelling results could be utilized to interpret experimental findings. The lattice parameters from QM calculations for K-montmorillonite are presented in Table 2. The correspondence between the calculated and experimental parameters was satisfactory (Tsipursky & Drits, 1984; Bérend et al., 1995). The calculations were performed with different approximations. The main difference between the approximations is how energy states of the system are described, especially partial densities originating from p-orbitals.
The Na- and Ca-montmorillonite structures were generated by replacing the K+ ions with Na+ ions to obtain Na-montmorillonite, and by replacing half of the K+ ions with Ca2+ ions and removing the remaining half of K+ ions to obtain Ca-montmorillonite. In this manner, the electrical neutrality of the structures was retained. The optimized structures and the optimized parameters are presented in Table 3.
In montmorillonites, the equivalent amount of exchangeable cations depends on the amount of Mg(II) atoms replacing Al(III) in the negatively charged octahedral sheets and the typical Mg:Al ratio is 1:5. In order to take into account the Mg:Al ratio in our studies and to retain the electrical neutrality of the models, the unit-cell structures had to be multiplied three times (1 × 3 × 1 unit cells in Fig. 2a). The positions of Mg atoms were selected so the distance between them is as long as possible based on the symmetry. This selection agrees well with previous attempts (Lavikainen et al., 2015). These new periodic models (Fig. 2b) were re-optimized and the optimized parameters for K-, Na- and Ca-montmorillonites are presented in Table 3.
Effect of water on the structures
Typically the so-called dry montmorillonite contains ∼8–10 wt.% interlayer water. Based on the molecular mass of water and montmorillonite without interlayer water, the unit cell of the montmorillonite contains 2–4 water molecules. This means that, according to QM and MD simulations, only one water molecular layer exists in the interlayer space of the montmorillonite.
When the water content exceeded 10 wt.%, the interlayer space between the TOT layers began to increase, causing the swelling of the mineral. In this study, swelling was investigated with both QM and MD methods and QM was also used to detect possible atomic-level changes in the structure of the TOT layers.
In the QM calculations, the water content varied from 0 to 10 water molecules in the unit cell corresponding to 0–246 mg of H2O/g of montmorillonite. Water molecules were added randomly into the interlayer space of the montmorillonite models one by one and the structure was re-optimized between each step. A corresponding technique was used in the earlier work based on Monte Carlo simulations (Zheng et al., 2011). The Na-montmorillonite structure with 98.2 mg of H2O/g of montmorillonite (∼10 wt.%) is shown in Fig. 2c, and its lattice parameters are a = 525.7 pm, b = 2 714.1 pm, c = 1 324.1 pm, α = 93.2°, β = 96.4°, γ = 89.9°, and the d value is 1 313.9 pm. In comparison to the structure without water molecules the change of the lattice parameters is most remarkable in the c-axis direction increasing the distance (d value) between the TOT layers. The water molecules appear as molecules in the interlayer space, with no reactions with the TOT layers. When ten water molecules were present in the unit cell of the montmorillonite, the water content was ∼24 wt.%, and the swelling was ∼33% in the c-axis direction (d value = 1747.6 pm) compared to the structure with 10 wt.% water. At least two molecular water layers exist between the TOT layers (Fig. 2d). More detailed analysis of the swelling results together with those obtained from MD calculations are presented later.
Hydration of cations
To understand the role of water in montmorillonite structures, the hydration of cations and clay layers has to be considered. Analysis of the coordination geometry of water molecules with cations reveals that free Na+ has tetrahedral (Fig. 3a) or octahedral (Fig. 3b) coordination geometry with water molecules. K+ also forms four- or six-coordinated complexes with water (Fig. 3c,d). The six-coordinated complex is trigonal prismatic. Free Ca2+ does not form four-coordinated complexes but occurs in octahedral form (Fig. 3e) or eight-coordinated square antiprismatic form (Fig. 3f). The hydration energies of Na+, K+ and Ca2+ cations are presented in Table 4. The coordination of the interlayer cations was analysed by Sposito et al. (1999) using Monte Carlo simulations.
In the interlayer space of montmorillonites, cations also tend to acquire corresponding coordination geometry like their free form in water solutions. The QM and MD simulations in the present study indicated that the interlayer water molecules form hydrogen bonds with the nearest cations and the oxygen atoms of the tetrahedral sheets. Na+ typically has four water molecules in its coordination sphere, and further Na+ can fulfil part of its coordination sites by forming van der Waals bonding with the oxygen atoms of the tetrahedral sheet. K+ favours a six-coordinated complex with five water molecules and one van der Waals bond to the tetrahedral sheet. Ca2+ also favours a six-coordinated complex, so there are three coordinated water molecules and three van der Waals bonds to the tetrahedral sheet. Each of these structures varies case-by-case, however, depending on the orientation of the randomly set water molecules in the interlayer space.
Swelling of montmorillonites relates to the strength of interaction forces of water molecules with other water molecules, cations and tetrahedral sheets within the interlayer spaces. This expansion process is controlled by the balance of attractive vs. repulsive electrostatic forces resulting from the presence of negatively charged montmorillonite layers and interlayer cations. A variety of different intermolecular and electrostatic forces are valid, including hydrogen bonding, van der Waals forces, double layer repulsions and ion-ion correlation effects (Shirazi et al., 2011).
Due to coulombic attractions between the negatively charged TOT layers and positively charged interlayer cations, the layers are initially held in close proximity to their neighbours. Van der Waals forces may also contribute to the total potential energy of attraction. The hydration energy of interlayer cations determines the total potential energy of repulsion, and thus the ease of penetration of water molecules into the interlayer space (Bérend et al., 1995). These sites remain closed until the expansion pressure of the surrounding water phase overcomes the resistance to interlayer space opening. This allows for a single layer of water molecules to enter and reside within the interlayer space, which in turn leads to a reduction in the potential energy of attraction as neighbouring layers are separated by a greater distance.
Following this initial phase, the hydration of interlayer cations may begin. The extent of cation hydration and swelling of montmorillonites is largely determined by the type of interlayer cations and their hydration energies (Bérend et al., 1995). In the case of the monovalent cations (Na+, K+) the degree of swelling is related to the hydration energy of the cations. However, hydration energies can vary significantly according to the degree of hydration of the cations (Table 4). Based on the greater hydration energies of the polyvalent cations (Ca2+), the initial hydration occurs faster than for the monovalent cations, because Ca2+ as a polyvalent cation may take up more water molecules into its coordination sphere (up to eight water molecules) than Na+ and K+ as monovalent cations (Fig. 3). Hence, the hydrogen bonding of the Ca2+ cations with the tetrahedral sheet is stronger than that of the Na+ and K+ cations. In order to add extra water into the interlayer space, the free energy loss of the incoming water must be greater than the work done in increasing the interlayer space. Therefore, Na- and K-montmorillonites can take up more water than Ca-montmorillonite.
The d spacing values of Na-, K- and Ca-montmorillonites with three different layer charges and varying interlayer water content were obtained from MD calculations where water was forced into the interlayer. The d values are given together with the QM results (using a layer charge of –0.66 e per unit cell), XRD results for Na-montmorillonite (Fu et al., 1990), for Ca-montmorillonite (Cases et al., 1997) and for K-montmorillonite (Bérend et al., 1995), and MC results from Zheng et al. (2011) for Na-montmorillonite with respect to the interlayer water content (Fig. 4). In each case, the MD results indicate that montmorillonites with smaller layer charges experience greater swelling. In K-montmorillonite, however, the difference in swelling between the three structures is the smallest.
Comparison of the MD results for Na-, Ca- and K-montmorillonites with the same layer charges and water contents shows that, in general, K-montmorillonite exhibits the largest d value with the exception of Ca-montmorillonite for water contents of ∼0.1–0.15 g/gclay. This Ca result is consistent with the findings based on the hydration energies of different cations (Table 4) and is also supported by the XRD results. Denis et al. (1991) reported the influence of K in montmorillonite as a swelling inhibitor which is not observed in the present study because water is forced into the interlayer (note the absence of external water). The present QM results indicated that Na-montmorillonite exhibits the largest d value of >∼0.2 g/gclay of water content. Also, QM predicts larger d values than MD, but, in general, similarities can be detected (mainly for Ca-montmorillonite). The trends in the swelling graphs are generally in good agreement with the results from XRD and MC (for Na-montmorillonite). Differences in the experimental and simulation results may be due to the differences in the structures, like the presence of tetrahedral substitutions in the experimental samples.
Examples of density distributions of the interlayer water oxygens and cations for Na-, Ca- and K-montmorillonites (layer charge 0.5 e per unit cell) are presented in Figures 5, 6 and 7, respectively. In these plots, the x axis is the basal distance from the bottom of the clay layer to the top of the interlayer space and the clay-layer thickness is ∼0.64 nm (calculated from atom centre to atom centre). All density curves of water oxygens show peaks indicating the formation of water layers. For some water contents these peaks are not as distinct (e.g.Fig. 5 – bottom left hand corner). They are considered here as ‘mixed states’ when the previously existing water layers are ‘full’ and before the formation of yet another clearly separate layer is favoured. For Na-montmorillonite, the second water layer began to form when the water content increased to six molecules per unit cell and the third layer began to form at water contents of 10–11 molecules per unit cell. These regions are represented by a step in the swelling graphs in Fig 4. For Ca-montmorillonite, two water layers began to form when the water content increased to 4–5 molecules per unit cell and three water layers with 13–14 water molecules per unit cell. For K-montmorillonite, two water layers started to form when the water content increased to six molecules per unit cell and three water layers when the water was 12 water molecules per unit cell.
The density distributions of Na in the plots indicate that the cations are surrounded by water molecules for most water contents, except at lower water contents (Fig. 5). The same behaviour was demonstrated by Ca-montmorillonite (Fig. 6). However, K cations are distributed more uniformly in the interlayer space indicating that they are not ‘bound’ between water layers (Fig. 7).
The radial distribution functions (RDF) of cation-water oxygens for one, two and three layer hydrates are presented in Fig. 8 for Na- and Ca-montmorillonites with layer charge of –0.5 e/unit cell and K-montmorillonite with layer charge of –1.0 e/unit cell. For Na-montmorillonite the first-neighbour peak was at ∼0.23–0.24 nm in each hydrate case. For Ca-montmorillonite, this peak was at ∼0.24–0.25 nm and for K-montmorillonite at ∼0.28–0.29 nm in each case, in accordance with Tao et al. (2010). The RDFs of water oxygen-water oxygen for the corresponding cases are shown in Fig. 9. In each case the first-neighbour peak was at ∼0.27–0.28 nm. However, the three-layer hydrate peak for Ca is shifted slightly when compared to the other two states.
Diffusion of water and cations
The MSDs of interlayer water (oxygens) and cations for Na-, Ca- and K-montmorillonites (Al10Mg2) with 10H2O/unit cell as a function of time are shown in Fig. 10. According to Einstein's relation the slope of the MSD curve gives a quantitative measure of the diffusion coefficient. The statistics of the time correlation data deteriorate after half of the 500 ps run and especially for cations the statistics suffer from using a small number of particles (Fig. 10). To obtain the slope, a least-squares fit was done to 40–80% of the data up to 250 ps.
The values of the diffusion coefficients obtained for water and cations were presented together with the MD results from Holmboe & Bourg (2014), Greathouse et al. (2015) and Boek (2014) (results reported for 3D diffusion and, thus, are lower than 2D values) in Figures 11 and 12, respectively. Note however, that calculation of the diffusion coefficients for cations suffers from use of small numbers of particles and, thus, only the values for montmorillonites with the larger layer charges (–0.66 and –1.0 e/unit cell) are reported. In general, the values increase with increasing water content and with decreasing layer charge. However, the increase was not linear. A rapid increase was observed in the values from monohydrate to bi-hydrate montmorillonite and, as the water content increased from approximately a half-saturated to a fully saturated hydrate structure, the diffusion coefficients decreased. In general, diffusion in Ca-montmorillonite is slower than in either the Na- or K-montmorillonites which is probably due to the different coordination and stronger interaction of the divalent Ca2+vs. the monovalent Na+ or K+. The results are in satisfactory agreement with those from previous MD studies considering the different structures used.
The present study indicated that both QM and MD can be used to study the crystalline swelling of montmorillonites. Based on the results for montmorillonites with only octahedral substitutions, smaller layer charge results in greater swelling, and the type of cation present in the interlayer also has a remarkable effect on the swelling properties. The MD simulations indicated that, in general, K-montmorillonite exhibits the largest d value, although the difference in the d value of different K-montmorillonites was minimal. The modelled swelling behaviour of K-montmorillonite contradicts the experimental observations where K acts as a swelling inhibitor. This is a result of the system setup. The simulations considered only the equilibrium state of a montmorillonite system where water is forced into the interlayer and the rate with which water is adsorbed was not taken into account. Similarities between the MD and QM results were detected, though in general QM predicts d values which are larger than those from MD. The likely reason reason is that QM calculations are static, and movement of atoms has not been taken into account as a function of time. Based on the MD calculations, the formation of water layers in the interlayer space of the montmorillonites, which is well known for natural samples, could be reproduced.
The data from MD calculations were also used to predict the self-diffusion coefficients of interlayer water and cations. The coefficient was not constant but increased non-linearly with increasing water content and with decreasing layer charge.
In the future, cation exchange reactions of montmorillonites in different salt concentrations should be studied in order to better describe the transport behaviour of water and cations in (and to) the interlayer space of montmorillonites. The hydration energies of cations in the interlayer space of montmorillonites and swelling pressures cannot be predicted based on these results, however.
The QM/MD techniques together are appropriate for studying the nanoscale interactions which determine macroscopic phenomena related to montmorillonites. In order to develop reliable continuum mechanics models, relevant nanoscale data are needed for input into these models.
The authors acknowledge the funding received from Posiva Oy.
|CASTEP Cambridge serial total energy package|
|DNP Double Numerical plus Polarization|
|GGA-PBE Perdew, Burke and Ernzerhof version of generalized gradient approximation functional|
|GGA-PW91 Generalized gradient approximation (Perdew-Wang)|
|MD Molecular dynamics|
|QM Quantum mechanics|
|CASTEP Cambridge serial total energy package|
|DNP Double Numerical plus Polarization|
|GGA-PBE Perdew, Burke and Ernzerhof version of generalized gradient approximation functional|
|GGA-PW91 Generalized gradient approximation (Perdew-Wang)|
|MD Molecular dynamics|
|QM Quantum mechanics|