System Setup for Molecular Simulations

A periodic 3D supercell containing two fully flexible Mnt layers was set up according to the previous modeling studies of hydrated Mnt (Skipper et al. Reference Skipper, Chang and Sposito1995a,Reference Skipper, Chang and Spositob; Tambach et al. Reference Tambach, Hensen and Smit2004; Tournassat et al. Reference Tournassat, Chapron, Leroy, Bizi and Boulahya2009; Bourg and Sposito Reference Bourg and Sposito2011; Holmboe & Bourg Reference Holmboe and Bourg2014; Marry et al. Reference Marry, Rotenberg and Turq2008). A single TOT (2:1 tetrahedral-octahedral-tetrahedral) layer in the simulation supercell consisted of 6×6×1 unit cells of Mnt. The supercell dimensions in the direction perpendicular to the basal plane were adjusted to accommodate the appropriate number of water molecules in the interlayer to simulate the hydration state of Mnt. The Mg–for-Al isomorphic substitutions in the octahedral sheet were distributed randomly, but excluded Mg-Mg pairs to maintain a mean layer surface charge density of –0.116 C/m², consistent with previous MD simulations (Tournassat et al. Reference Tournassat, Chapron, Leroy, Bizi and Boulahya2009). A representative snapshot (Fig. 1) shows the simulation supercell of Na-Mnt in the 5W hydration state with z dimension 50.8 Å, containing 52 Na ions, 4 Cl ions, and 1800 H₂O water molecules. In the following discussion, the interlayer size is defined as the distance between the basal planes of the outermost oxygen atoms of Mnt.

Fig. 1. Snapshot of a MD simulation cell showing Na-Mnt in the 5W hydration state with z dimension 50.8 Å containing Na ions (dark blue), Cl ions (light blue), and water molecules (O and H atoms are in red and white, respectively). The clay layer consists of an octahedral sheet (pink) sandwiched by tetrahedral sheets (gold yellow)

Classical MD Simulations

MD simulations were performed using the open source package GROMACS (Berendsen et al. Reference Berendsen, van der Spoel and van Drunen1995; van der Spoel et al. Reference van der Spoel, Lindahl, Hess, Groenhof, Mark and Berendsen2005) which uses the CLAYFF force field (Cygan et al. Reference Cygan, Liang and Kalinichev2004) for montmorillonite, the ion-potentials reported by Åqvist (Reference Åqvist1990) for Ca²⁺, Joung-Cheatham parameters (Joung & Cheatham III Reference Joung and Cheatham III2008) for monovalent ions (Na⁺/Cl^–), and the SPC/E water model (Berendsen et al. Reference Berendsen, Grigera and Straatsma1987) for all interatomic interactions. The Lennard-Jones shift potential was truncated beyond a cutoff value of 1.0 nm. The long-range electrostatic interactions were computed using the Particle-Mesh-Ewald (PME) method (Darden et al. Reference Darden, York and Pederson1993). All atomic simulations were conducted at 298 K using a time step of 1 fs. The integration of the equations of motion was performed using the leap-frog algorithm. Energy minimization with the steepest descent algorithm was used to prepare the initial configurations. The equilibration was performed in two steps. First, the system was pre-equilibrated in the NPT ensemble (pressure P, number of particles N, and temperature T were fixed) for 5 ns using the Berendsen barostat (Berendsen et al. Reference Berendsen, Postma, van Gunsteren, DiNola and Haak1984) at 298 K and 1 bar to relax the Mnt layers and interlayer electrolytes in all the dimensions. The 5 ns-long NVT (fixed volume V, N, and T) runs were performed using position restraints of the clay lattices to generate equilibrated systems. The x and y dimensions of the NVT supercell in all simulated hydration states were kept at 31.0 Å×54.2 Å. The dimensions of the supercell in the z direction and the composition of the interlayer are given in Table 1. The production runs were performed for 125 ns in the NVT ensemble at 298 K.

Table 1. Simulation parameters of the supercell for the MD simulations

Fluid DFT Calculations

The f-DFT simulations were conducted using the open source package Tramonto developed at Sandia National Laboratories (https://software.sandia.gov/tramonto/). The f-DFT simulations were performed in grand canonical ensemble (μVT), where the total volume V, the temperature T, and the chemical potentials μ of the system were kept constant. Classical MD used in this study was implemented in the NVT ensemble. Therefore, the assignment of a consistent setup for f-DFT and MD simulations and μVT and NVT ensembles was not trivial. To match the condition of f-DFT and MD simulations two different approaches for the system were used at low (<5 water layers) and high (≥5 water layers) degrees of hydration. At low degrees of hydration, the total number of ionic species and water molecules in the interlayer space were set to be the same in both f-DFT and MD systems. At high degrees of hydration the chemical potential in f-DFT was set to match the density representing the regions of homogeneous fluid in the central part of the interlayer space.

In f-DFT, the density distributions of molecular species of fluid at equilibrium, ρ _i (r), were determined by minimization of the grand potential energy at equilibrium described by the Euler-Lagrange equation. The chemical potential was expressed as a function of the external potential, Coulomb forces, ideal gas components, and excess energy due to hard-sphere repulsions ( $f_i^{\mathrm{hs}}$ ), Lennard-Jones attractions ( $f_i^{\mathrm{LJ}}$ ), and short-range electrostatic interactions ( $f_i^{\mathrm{MSA}}$ ) were evaluated with the Mean Spherical Approximation (MSA) (Waisman & Lebowitz Reference Waisman and Lebowitz1972):

(1) ${\mu }_i\left(r\right)=q_{i}e\psi \left(r\right)+u_i\left(r\right)+kTln{\rho }_i\left(r\right)+f_i^{\mathrm{hs}}\left(r\right)+f_i^{\mathrm{LJ}}\left(r\right)+f_i^{\mathrm{MSA}}\left(r\right)$

in which ψ(r) is the local electrostatic potential resulting from the Coulombic interaction between point charges and that from the electrostatic part of the external potential, which was obtained by solving Poisson’s equation, q _i is the valence of the ionic species i using the expression

(2) ${\nabla }^2\psi \left(r\right)=-\frac{1}{{\epsilon }_0{\epsilon }_r}{\sum }_i{\rho }_i\left(r\right)e{z}_i$

where ε₀ is the vacuum permittivity, ε _r is the relative permittivity or dielectric constant, z _i is the valence of the ionic species i, and e is the elementary charge. The external field, u _i (r), is induced by the clay surface, which is described as a flat wall with surface charge density σ uniformly distributed on the surface, k is Boltzmann’s constant, T is the absolute temperature, and f _i (r) is the Helmholtz free energy per molecule of the particle i.

The one-dimensional external potential due to the wall–fluid interaction was represented by the Lennard-Jones 9-3 potential:

(3) ${\mu }_i\left(z\right)={\epsilon }_{\mathrm{wf}}^{\mathrm{LJ}}\left(\frac{2}{15}{\left(\frac{d}{z}\right)}^9-{\left(\frac{d}{z}\right)}^3\right)$

where ${\epsilon }_{\mathrm{wf}}^{\mathrm{LJ}}$ is the LJ energy interaction parameter for wall–fluid interaction potential, d is the LJ particle diameter, and z is the position of molecular particle i from the wall in the direction perpendicular to the wall. The wall-fluid potential cutoff was set at 3.5 d.

The LJ interaction is described in a spirit of perturbation theory as hard-sphere interaction with an effective hard-sphere parameter d and the correction due to long-range interactions (Tang et al. Reference Tang, Scriven and Davis1991). The hard-sphere component, $f_i^{\mathrm{hs}}\left(r\right)$ , was calculated using the fundamental measure theory (FMT) developed by Rosenfeld (Reference Rosenfeld1989). The short-range interaction cutoff for $f_i^{\mathrm{hs}}\left(r\right)$ was 1d. The excess free energy was based on the following formula,

(4) $f_i^{\mathrm{hs}}\left(r\right)=-n_{0}ln\left(1-n_3\right)+\frac{n_1{n}_2-n_{v1}.n_{v1}}{1-n_3}+\frac{1}{24\pi }\frac{n_2^3}{{\left(1-n_3\right)}^2}-\frac{1}{8\pi }\frac{n_2\left(n_{v2}.n_{v2}\right)}{{\left(1-n_3\right)}^2}$

Note that for a complete description including non-local densities, n and n_v, as well as an analysis of the functionals, the reader is referred to the original references (Rosenfeld Reference Rosenfeld1989).

The excess part of free energy due to the LJ attraction with respect to a hard-sphere reference fluid is defined by a density-weighted integral of a pair potential function over the surrounding fluids, separately summing the contributions from each species,

(5) $f_i^{\mathrm{LJ}}\left(r\right)={\sum }_j{\int }_{1d}^{3.5d}4{\epsilon }_{\mathrm{ij}}^{\mathrm{LJ}}\left({\left(\frac{d}{r}\right)}^{12}-{\left(\frac{d}{r}\right)}^6\right){\rho }_j\left(r\prime \right)dr\prime$

where ${\epsilon }_{\mathrm{ij}}^{\mathrm{LJ}}$ is the LJ energy interaction parameter for fluid–fluid interaction potential, r = |s − r ^′| refers to the distance between particles i at position s and j at r ^′, and r is the cutoff distance. This contribution was considered to be within 1 to 3.5 molecular diameters, while it is set to zero outside this region.

The contribution to the free energy accounting for the short-range electrostatic interactions is expressed as

(6) $f_i^{\mathrm{MSA}}\left(r,{\rho }_i\right)={\sum }_j{\int }_0^{1d}{c}_{\mathrm{ij}}^{\left(2\right),\mathrm{el}}\left(r,r\prime \right){\rho }_j\left(r\prime \right)dr\prime$

in which $c_{\mathrm{ij}}^{\left(2\right),\mathrm{el}}\left(r,r\prime \right)$ is the cross-correlation between the Columbic interaction and the hard-sphere exclusion, i.e. the electric residual part of the pair direct correlation function defined as the second functional derivative of the excess part of the free energy functional with respect to the single-particle density distribution (Hansen & McDonald Reference Hansen and McDonald2006). It can be modeled with MSA (Waisman & Lebowitz Reference Waisman and Lebowitz1972) for homogeneous systems, leading to

(7) $c_{\mathrm{ij},\mathrm{MSA}}^{\left(2\right),\mathrm{el}}\left(r\right)=\frac{q_i{q}_j{e}^2}{2\pi {\epsilon }_0{\epsilon }_r}\left(\frac{B}{d}-{\left(\frac{B}{d}\right)}^2\frac{r}{2}-\frac{1}{2r}\right),r<d$

with the B parameter written as B = Γd/(1 + Γd), where $\Gamma =\frac{\sqrt{1+2\mathrm{\kappa d}}-d}{2d}$ , ${\kappa }^2=\frac{e^2}{{\epsilon }_0{\epsilon }_r{k}_{B}T}{\sum }_i{q}_i^2{\rho }_i^b$ , ${\rho }_i^b$ is the density of species i in bulk reservoir, q _i is the valence of the ionic species i, and e is the elementary charge. The relative permittivity depends, in general, on the system temperature and composition.

Fluid Structure at the Interface

The interaction parameters used in the f-DFT calculations are reported in Table 2. In the DFT/LJ-3CM results, the structural and thermodynamic properties such as density profiles, cation selectivity of ion exchange equilibrium, and free energy were fitted to reproduce atomistic simulations (Marry et al. Reference Marry, Rotenberg and Turq2008; Yang et al. Reference Yang, Neretnieks and Holmboe2017b), and the LJ energy parameters for fluid–fluid and wall–fluid were taken from molecular simulations (Magda et al. Reference Magda, Tirrell and Davis1985; Lee et al. Reference Lee, Nilson, Templeton, Griffiths, Kung and Wong2012). The uniform surface charge density in DFT//LJ-3CM was set to be the same as in MD simulations, –0.116 C/m². To test the performance of the DFT/LJ-3CM in comparison with the DFT/HS-3CM, i.e. the solvent primitive model, the results for both approaches were compared with MD simulation for Mnt in the 10W hydrate state. HS radii used in the DFT/HS-3CM were the same as the LJ parameters in Table 2.

Table 2. The interaction energy parameters for the DFT/LJ-3CM simulation

The DFT/LJ-3CM results taking into account the solvent effects were compared first with the atomic simulations for the case when Na-Mnt was saturated in a 0.62 M NaCl bulk solution. Density distributions were compared among water H₂O (i.e. water oxygen OW in MD simulations), Na⁺ and Cl^– in the interlayer obtained with DFT/LJ-3CM, DFT/HS-3CM models, and atomic simulations (Fig. 2). The structured density profiles of interlayer water and ions from atomic simulations were predicted accurately by the DFT/LJ-3CM results. DFT/HS-3CM calculations overestimated, by nearly four times, the contact cation (Na⁺) and water concentrations. The excessive cation density at the contact resulted in overcharging of the surface and values which were too large for co-adsorption of anions (e.g. Cl^–) at the interface. In the DFT/HS-3CM model, the solvent was modeled by neutral spheres instead of a dielectric continuum as used in DFT/HS-2CM. Therefore, the wall/ion electrostatic interaction would be overestimated in the DFT/HS-3CM model and this could explain the overestimation of the densities at contact.

Fig. 2. Density distributions of Na⁺ and Cl^– ions and H₂O from DFT/LJ-3CM calculations as a function of distance from the surface in the 10W hydration state. The DFT/LJ-3CM results are compared with atomic simulations (gray curves) and with calculations using the DFT/HS-3CM (blue curves) and DFT/HS-2CM (dark blue curves) models

A similar comparison was conducted using the DFT/LJ-3CM model for Mnt saturated in a 0.23 M CaCl₂ bulk solutions for Ca-Mnt in the 10W hydration state. The density distributions of calcium and chloride ions and water molecules as a function of distance from the surface were obtained with the DFT/LJ-3CM model and from atomistic simulations. The DFT/LJ-3CM results agreed well with the atomistic simulations. The DFT/LJ-3CM model successfully captured the damped oscillations of Ca²⁺ and H₂O density distributions in the vicinity of the surface (Fig. 3). These important features were missing from the density profile obtained with the DFT/HS-2CM model. Furthermore, the number of density profiles of ionic and solvent molecules in the vicinity of the charged surface were substantially overestimated by the DFT/HS-3CM model. The mismatch of anion density predicted by the DFT/HS-3CM model in the 2:1 (e.g. CaCl₂) was even larger than in the 1:1 (e.g. NaCl) electrolyte. This can be explained by the stronger electrostatic effect of the Ca²⁺ ion compared to Na⁺. For the same (overestimated) cation atom density, Ca²⁺ resulted in stronger co-adsorption of Cl^– anions at the surface. The comparison demonstrated that the DFT/LJ-3CM model clearly performed better than the DFT/HS-3CM model near the fluid–wall interface, especially for 2:1 electrolytes. Accordingly, further comparison of results from the molecular simulation with the DFT/LJ-3CM model was restricted.

Fig. 3. Density distributions of Ca²⁺, Cl^– ions, and H₂O from the DFT/LJ-3CM calculations as a function of distance from the surface in the 10W hydration state. The DFT/LJ-3CM results are compared with atomic simulations (gray curves) and with calculations using the DFT/HS-3CM (blue curves) and the DFT/HS-2CM (dark blue curves) models

To test further the accuracy of the DFT/LJ-3CM model, the f-DFT results were compared with atomistic simulations of Na-Mnt in NaCl solution with various water contents, i.e. 2W, 3W, 4W, and 5W hydrated states (Holmboe & Bourg Reference Holmboe and Bourg2014; Yang et al. Reference Yang, Neretnieks and Holmboe2017b). The anions were essentially excluded from the nanopore when the distance between Mnt layers was <10 Å. The f-DFT results were insensitive to the bulk concentration of NaCl when the nanopore between Mnt layers was <10 Å. Accordingly, the concentration of bulk NaCl solutions in DFT calculations was set to be 0.6 M for Mnt in 2W, 3W, 4W, and 5W hydrated states.

The density profiles in a ~6.4 Å interlayer pore (Fig. 4 ) were obtained from the DFT/LJ-3CM model and atomistic simulations corresponding to the 2W hydration state. The DFT/LJ-3CM model predicted two symmetrically peaked, cation-density profiles whereas atomistic simulations predicted single-peaked Na⁺ density in the middle of the interlayer (Fig. 4 ). At such a low hydration state, the observed structuring of the fluid was controlled by the strong orientation-dependent ion–solvent and solvent–surface interactions. The observed difference was clear evidence that important features such as H-bonding and rotation of water molecules were not taken into account by the DFT/LJ-3CM model.

Fig. 4. Comparison of density distributions of Na⁺ (dark blue) and water (light blue) predicted by the DFT/LJ-3CM model with atomic simulations (black curves) as a function of distance from the Mnt surface in 2W hydration state. The dashed gray lines represent the position of the outermost oxygen sites on the Mnt surface. The concentration of Cl^– ions was negligibly low and is not shown here for clarity

Low temperature ion-dipole interaction cannot be captured by the orientation-independent LJ potential. However, the main features of water density distributions were captured successfully by the DFT/LJ-3CM results (Fig. 4). Atomistic simulations also indicated a small fraction of Na⁺ forming a small local density maximum at ~1.5 Å from the surface. These were inner-sphere complexes entering the hexagonal cavities of the tetrahedral sheet of Mnt. Such complexes were also present at the greater degrees of hydration discussed later. These features cannot be considered in the current DFT/LJ-3CM model which uses a simple planar interface.

The density distributions of calcium ions and water molecules in a ~9 Å interlayer pore were predicted by the DFT/LJ-3CM model and by atomistic simulations for Mnt in the 3W hydration state (Fig. 5). The DFT/LJ-3CM model captured the features of the density profiles of cations, especially the structure of interlayer water molecules in the 3W hydration state.

Fig. 5. Comparison of density distributions as a function of distance from the surface in the 3W hydration state from the DFT/LJ-3CM calculations (Na ion density is dark blue and H₂O density is light blue) with atomic simulations (black curves)

Again, Na and H₂O density profiles were predicted by the DFT/LJ-3CM model and the atomistic simulations of Mnt in the 4W hydrate state (Fig. 6). Similar Na and H₂O density profiles were predicted by the DFT/LJ-3CM model and the atomistic simulations for the 5W hydration state (Fig. 7). The density profiles of sodium ions predicted by the DFT/LJ-3CM model were more structured than profiles from MD simulations. Both the DFT/LJ-3CM results (blue curves) and the atomistic simulations (black curves) showed that Na⁺ ions were enriched at a distance of approximately 0.5d from the surface. This is comparable to the typical ion–surface distance of Na-outer-sphere complexes. As for the 3W state, DFT/LJ-3CM slightly over-predicted the structuring of the ion and water in the middle of the interlayer. The secondary features in the ions and water density distributions related to the non-planar surface topography were not reproduced by the f-DFT model. These secondary density maxima are related to partially hydrated inner-sphere Na ions entering the hexagonal cavities on the Mnt surface. Such a complex solvation behavior cannot be captured by the DFT/LJ-3CM model due to the over-simplistic planar representation of the surface and the lack of orientation-dependent ion–dipole interaction.

Fig. 6. Comparison of density distributions as a function of interlayer separations in the 4W hydration state obtained using the DFT/LJ-3CM model (Na ions are dark blue and H₂O is light blue) with atomic simulations (black curves)

Fig. 7. Comparison of density distributions as a function of interlayer separations in the 5W hydration state from the DFT/LJ-3CM model (Na ions in dark blue and H₂O in light blue) with atomic simulations (black curves)

Overall analysis of the data demonstrated that the DFT/LJ-3CM model could capture the primary structural properties of ions and water in the vicinity of a fluid surface interface for very different hydration states of Mnt. The model tended to over-predict the structuring at greater distances from the surface in comparison with the MD simulation results. The features related to the non-planar interface and the changes in the hydration state of the ions, because of the changes from outer-sphere to the inner-sphere complexes, cannot be captured within the currently used DFT/LJ-3CM model.

Swelling Behavior

Depending on the chemical potential of water, smectite minerals are known to contain a distinct number of water molecules in the interlayer, commonly referred to as 1W, 2W, and 3W hydration states (Mooney et al. Reference Mooney, Keenan and Wood1952; Karaborni et al. Reference Karaborni, Smit, Heidug, Urai and Oort1996; Kawamura et al. Reference Kawamura, Ichikawa, Nakano, Kitayama and Kawamura1999). The stability of a particular hydration state and the exact interlayer thickness depends on temperature, composition of the electrolyte solution (both concentration and ion speciation in particular), and the distribution of isomorphic substitution in the TOT layer (Boek et al. Reference Boek, Coveney and Skipper1995a,Reference Boek, Coveney and Skipperb; Young & Smith Reference Young and Smith2000; Hensen & Smit Reference Hensen and Smit2002; Smith et al. Reference Smith, Wang, Chaturvedi and Whitley2006; Tambach et al. Reference Tambach, Bolhuis, Hensen and Smit2006). A model-based description of smectite swelling is essential for prediction of a clay’s mechanical properties (Hoang & Galliero Reference Hoang and Galliero2015; Sun et al. Reference Sun, Hirvi, Schatz, Kasa and Pakkanen2015). The stable hydration state corresponds to the minimum value of the free energy of the system. The dependence of surface free energy, Ω(h), on the interlayer separation, h, is calculated by an integration of the solvation force via (Frink & van Swol Reference Frink and van Swol1998):

(8) $\Omega \left(h\right)-2\ast \Omega \left(\infty \right)=-{\int }_{\infty }^{h}d{h}^{\prime }{f}_s\left(h\prime \right)$

The constant Ω(∞) is equal to the surface free energy of an isolated wall in equilibrium with bulk fluid. The factor 2 in the left hand side of Eq. 8 accounts for the presence of two interfaces in the slit pore. The absolute minimum in Ω(h) determines the thermodynamically stable equilibrium separation of two walls, separated by h, for an unconstrained swelling condition, where the solvation force f _s is defined as f _s = p − p _b , and p is the total normal force per unit area on each wall, whereas p _b is the bulk pressure at the chosen chemical potential and temperature.

The dimensionless free energy of the system for Na- and Ca-Mnt was plotted as a function of slit pore size for various ionic strengths (Fig. 8). The symbol A in the vertical axis represents the surface area. Depending on the ionic strength of the external electrolyte, the stable states of the Ca and Na-Mnt predicted by DFT/LJ-3CM were either a dispersed one at lower ionic strength or a compacted “3W-Layer” state at higher ionic strength. The model for Na-Mnt predicted the dispersion to be stable if the bulk concentration of monovalent ions was <2 M (Fig. 8a). The model for Ca-Mnt predicted the dispersion state to be stable at bulk concentrations of 2:1 electrolyte below 0.01 M (Fig. 8b). The metastable states with 2W and 4W layers were also apparent in the free energy profiles as local minima. The prediction of stable interlayer distances was based on changes in the chemical potential of each ion in the bulk electrolyte. The 2W state can become stable over the 3W state if the chemical potential of solute is further reduced, e.g. partial saturation (Churakov Reference Churakov2013). The behavior of the DFT/LJ-3CM for Mnt agreed very well qualitatively with the behavior of natural systems. The quantitative differences were related to the over-simplistic description of solvent and solvent–solute interactions, which were not captured in the model. Contrary to the experimental observations, Poisson-Boltzmann approximations predict repulsive interaction of Mnt platelets in NaCl and CaCl₂ solutions (Segad et al. Reference Segad, Jonsson, Akesson and Cabane2010). The 2-component RPM model improved the predictions, suggesting an attractive force between TOT layers in Ca-Mnt at ~10 Å separation. In contrast to the DFT/LJ-3CM, the two-component RPM did not predict the existence of other, eventually metastable, hydration states. In this context the DFT/LJ-3CM model used in the present study clearly performed better than the widely used DFT-based description of swelling in clay minerals.

Fig. 8. Free energy per unit surface area as a function of the equilibrium interlayer separation between two smectite TOT layers in a NaCl, and b CaCl₂ bulk solutions. The ionic concentrations are indicated in the plots

Thermodynamics of Ion Sorption

In this section, the DFT/LJ-3CM model was applied to evaluate ion exchange equilibrium and to estimate the effect of temperature on cation selectivity. The cation exchange between Ca²⁺ and Na⁺ ions in the interlayer of Mnt can be expressed by the following reaction (Appelo & Postma Reference Appelo and Postma2005):

(9) $\frac{1}{2}{\mathrm{Ca}}^{2^{+}}+\mathrm{Na}-X\leftrightarrow \frac{1}{2}\mathrm{Ca}-X_2+{\mathrm{Na}}^{+}$

where X indicates the clay mineral exchanger. The selectivity coefficient of ion exchange K _GT is given by the Gaines-Thomas convention (Gaines & Thomas Reference Gaines and Thomas1953; Appelo & Postma Reference Appelo and Postma2005) as

(10) $K_{\mathrm{GT}}=\frac{\left(N{a}^{+}\right){\left(\mathrm{Ca}-X_2\right)}^{0.5}}{{\left(C{a}^{2+}\right)}^{0.5}\left(\mathrm{Na}-X\right)}$

where square brackets denote activities. [Na ⁺] and [Ca ²⁺] represent the activities of ions Na⁺ and Ca²⁺ in water, which can be related to the ion concentrations in bulk aqueous solutions. The surface activities [Na − X] and [Ca − X ₂] are approximated by the equivalent charge fractions of Na⁺ and Ca²⁺ ions in the adsorbed phase in the DFT calculations. More detailed discussions of surface activities and ion exchange properties of clay minerals is available elsewhere (Yang et al. Reference Yang, Neretnieks, Moreno and Wold2017a).

The selectivity coefficient K _GT as defined with the Gaines-Thomas convention was obtained by the DFT/LJ-3CM and the experimental results (Sposito et al. Reference Sposito, Holtzclaw, Charlet, Jouany and Page1983) of the surface activities [Na − X] i.e. β_Na = τ _Na /(τ _Na + 2τ _Ca ), where ${\tau }_i={\int }_0^h\left({\rho }_i\left(x\right)-{\rho }_i^b\right)\mathrm{dx}$ (Table 3). The selectivity coefficients K _GT obtained previously using the DFT/HS-3CM model (Yang et al. Reference Yang, Neretnieks, Moreno and Wold2017a) are shown for comparison.

Table 3. Thermodynamics of Ca-Na exchange in Wyoming Mnt predicted using the DFT/LJ-3CM model (present study), the DFT/HS-3CM model (Yang et al. Reference Yang, Neretnieks, Moreno and Wold2017a), and the experimental data of Sposito et al. (Reference Sposito, Holtzclaw, Charlet, Jouany and Page1983)

To demonstrate the ability of the DFT/LJ-3CM to describe the cation exchange in Ca-Na-Mnt, the corresponding data in Table 3 were plotted (Figs. 9 and 10), and revealed (Fig. 9) the equivalent charge fractions of Na ions in the adsorbed phase, i.e. surface activities of Na-X (β_Na) from the experiments and the DFT/LJ-3CM calculations as a function of the mole fraction of Na ions in the mixture of NaCl and CaCl₂ bulk solutions. The calculations were performed at room and higher temperatures taking into account the change in the relative dielectric permittivity of water (Malmberg & Maryott Reference Malmberg and Maryott1956). The dielectric constant in the interlayer and in the bulk solution was assumed to have the same temperature dependence. The model predicted that temperature has a very limited effect on the cation exchange in Na-Ca-Mnt. Na sorption decreased only slightly with increasing temperature. The comparison of f-DFT calculations and experiments showed that the DFT/LJ-3CM model can quantitatively and qualitatively reproduce the measurements of cation exchange for Mnt in the Na and Ca electrolyte solutions.

Fig. 9. Surface activities as a function of the mole fraction of Na⁺ ions in the mixture of NaCl and CaCl₂ bulk electrolytes

Fig. 10. The selectivity coefficient K _GT as a function of the mole fraction of Na ions in the mixture of NaCl and CaCl₂ bulk solution

The selectivity coefficient K _GT from the DFT calculations and measurements, represented as a function of the mole fraction of NaCl in the mixture of NaCl and CaCl₂ bulk solutions using the DFT/LJ-3CM (Fig. 10), revealed a distinct advantage over the RPM calculations in comparison with the experimental results (Sposito et al. Reference Sposito, Holtzclaw, Charlet, Jouany and Page1983). The DFT/LJ-3CM results agreed quite well with the experiments when the mole fraction of Na ions in the bulk was >0.4. The interlayer separations in the f-DFT calculations of the selectivity coefficient were set at 30 Å when Na was dominant (i.e. Na mole fraction >0.7), which is commonly used in simulations for pure Na electrolyte solutions because the Mnt layers tend to disperse. However, the Mnt layers collapsed and formed a gel with ~9 Å interlayer separation (the compacted ‘3W-Layer’ state) for Ca salts with bulk ion concentration ≥10 mM (Fig. 8). Accordingly, the interlayer separations in the DFT calculations were set to be 9 Å when Ca was dominant in the bulk solutions (i.e. Na mole fraction <0.7). The selectivity coefficient log₁₀ K _GT (Fig. 11) was plotted as a function of the Na charge fraction at the surface, i.e. the surface activities in the adsorbed phase. A linear relationship consistent with the experimental observations was recovered between the selectivity coefficient and the surface activities.

Fig. 11. The selectivity coefficient, log₁₀ K _GT , as a function of the charge fraction of the Na⁺ ions in the adsorbed phase. The curve shows the DFT results. The solid dots represent the experimental data from Sposito et al. (Reference Sposito, Holtzclaw, Charlet, Jouany and Page1983)