Optimized Geometry and Total Potential Energy of Different Configurations

Atomic and bond details of the resultant 15 montmorillonite structural models from this DFT study are shown in the Supplementary Materials section. Configuration A with alternate locations for Fe atoms is shown as an example of an optimized structure in Fig. 1. Slight variations were noted in the clay structures with different octahedral cation distributions (Table 2). The a dimension of cv configurations (10.47 ± 0.02 Å) were smaller (p < 0.001) than the tv ones (10.54 ± 0.01 Å). The b dimension of cv configurations (9.10 ± 0.01 Å) were larger (p < 0.001) than the tv ones (9.06 ± 0.02 Å), but the c dimension (9.57 ± 0.06 Å) was not significantly different for all configurations. The α angles of cis-Mg configurations (90.08 ± 0.14^ο) were larger (p < 0.001) than the trans-Mg ones (89.49 ± 0.05^ο). The β (91.01 ± 0.26^ο) and γ (89.98 ± 0.10^ο) angles showed no significant differences among the configurations. The slight variations observed were due to the different ionic radii of the octahedral cations in each clay structure: Al³⁺ (0.61 Å), Mg²⁺ (0.80 Å,) and Fe³⁺ (0.73 Å). Similar results were also obtained in DFT simulations of smectites where the a, c, and β parameters of cv structures were smaller than tv ones (Sainz-Diaz et al. Reference Sainz-Diaz, Hernández-Laguna and Dove2001).

Fig. 1. A [001] (A) and [100] (B) projection of cis-vacant cis-Mg montmorillonite supercell model with the top plane of silica removed for clarity for the [001] projection. The positions of Mg and H⁺ are kept fixed in all models; Fe replaces Al in the octahedral sheet and is indexed by the number of octahedral jumps (j) separating Fe from Mg. The digits (1 to 3, 1* and 2*) indicate the several octahedral jumps of Fe positioned with respect to the Mg within the clay structure.

Table 2. The optimized main lattice parameters (length in Ångstroms and angles in degrees)

These structural and elemental variations resulted in internal energy stresses expressed by the total potential energy of the clay (Table 3). The average of all the total potential energies of the structures with and without charge-compensating H⁺ were –587.90 eV and –583.62 eV, respectively. The changes of total potential energy were within 0.42 eV and 0.52 eV for configurations with (E _w/H) and without (E _w/o H) charge-compensating H⁺. These values were within 0.07% and 0.09% of the average total potential energy of structures with and without charge-compensating H⁺, respectively. These energy differences of different configurations were very small considering the size of the supercell structure studied (80 or 81 atoms). These small energy differences came from the interactions among Al³⁺, Mg²⁺, and Fe³⁺ cations as first, second, and third neighbors in the octahedral sheet. Small energy differences of simulated clay structures with different dioctahedral cation distributions were also observed in other DFT studies (Sainz-Diaz et al. Reference Sainz-Diaz, Hernández-Laguna and Dove2001; Hernández-Laguna et al. Reference Hernández-Laguna, Escamilla-Roa, Timón, Dove and Sainz-Díaz2006; Tunega et al. Reference Tunega, Goodman, Haberhauer, Reichenauer, Gerzabek and Lischka2007; Ortega-Castro et al. Reference Ortega-Castro, Hernández-Haro, Dove, HernándezLaguna and Sainz-Díaz2010; Hernández-Haro et al. Reference Hernández-Haro, Ortega-Castro, Pruneda, Sainz-Díaz and Hernández-Laguna2014; Ohkubo et al. Reference Ohkubo, Okamoto, Kawamura, Guégan, Deguchi, Ohki, Shimizu, Tachi and Iwadate2018).

Table 3 The optimized total potential energy of configurations with (E _w/H) and without (E _w/o H) the charge-compensating H⁺. The relative total potential energy (∆E _w/H) of configurations with H⁺ cations was calculated by using the minimum energy of configuration R as 0 eV. The HAE indicate the H⁺ adsorption energy in H-rich or H-poor environments for each configuration. Configurations correspond to those shown in Table 1.

The relative energy was obtained by setting the minimum energy of structure R as 0 eV. The average of the relative energy change (ΔE _w/H) of cis-vacant cis-Mg structures (0.34 ± 0.05 eV) revealed no significant difference from the cis-vacant trans-Mg structures (0.30 ± 0.07 eV), but both were significantly larger (p < 0.001) than the trans-vacant cis-Mg structures (0.13 ± 0.08 eV) (Fig. 2). Therefore, the tv configurations were more stable than the cv ones, and the relative position of octahedral hydroxyls (cis-Mg vs trans-Mg) did not influence the stability of the structures. Using DFT simulations, Tunega et al. (Reference Tunega, Goodman, Haberhauer, Reichenauer, Gerzabek and Lischka2007) also showed that the tv structures were more stable than the cv structures when isomorphic substitution occurred in the octahedral sheet.

Fig. 2. The optimized total relative potential energy ΔE _w/H with the charge-compensating H⁺ for each configuration. Configurations correspond to those shown in Table 1. The ▲, ■, and ♦ symbols represent configurations with octahedral jumps from Fe to Mg as 1, 2, and 3, respectively.

Additionally, these small energy differences suggest that cv and tv structures can be transformed from one to the other, which is substantiated in the literature. Dehydroxylation-rehydroxylation reactions of two dioctahedral micas (celadonite and glauconite) resulted in the transformation of the tv structure into the cv structure (Muller et al. Reference Muller, Drits, Plançon and Robert2000). Similarly, a cv montmorillonite was transformed into a tv structure by dehydroxylation at high temperature (510–630°C) and subsequent rehydroxylation under steam between 1.5 and 8.0 MPa (Emmerich and Kahr Reference Emmerich and Kahr2001). Using hydrothermal treatments, the proportion of cv and tv layers in synthesized low-charge beidellites can be modified by controlling the synthesis pressure, temperature, and starting materials (Lantenois et al. Reference Lantenois, Muller, Bény, Mahiaoui and Champallier2008). This also suggests that the probability of occurrence of natural smectites with similar composition and structure should be nearly equal.

The octahedral Fe-Mg distance did not significantly influence the stability of the various configurations, where the relative total potential energies were 0.21 ± 0.12, 0.26 ± 0.15, and 0.29 ± 0.08 eV for j = 3, 2, and 1, respectively (Table 3, Fig. 2). Therefore, no significant preference should occur for the formation of any specific structure with a specific octahedral Fe-Mg distance. That is, the octahedral Al, Fe, and Mg atoms in this mineral should be randomly distributed. Accordingly, the relative possibilities of the formation of configurations with j = 3, 2, and 1 were 1/7, 3/7, and 3/7, respectively, regardless of the presence of cv or tv, and cis- or trans-Mg structures.

Influence of Fe Octahedra on the H⁺ Adsorption Strength

The H⁺ adsorption energies (HAE) were calculated for two H⁺ environments: H-rich and H-poor conditions (Table 3). The HAE of configurations in the H-poor environment were 1.98 eV higher than the HAE values in the H-rich environment because of the difference in the μ_H value for H-rich versus H-poor environments. The negative values of HAE in the H-rich environment indicated that the H⁺ adsorption was an exothermic reaction, and the clay preferred to adsorb the added H⁺. The positive HAE values in the H-poor environment indicated that the H⁺ adsorption was an endothermic reaction and the H⁺ adsorption strength was very weak. The more negative HAE values corresponded to higher H⁺ adsorption stability.

A strong linear correlation of HAE with the direct Fe-Mg distance was found for both H-rich and H-poor environments (R ² = 0.83, Fig. 3). The results were separated into three distinct groups (labeled 1, 2, and 3) based on the number of octahedral jumps (j) from Fe to Mg and HAE values. The average HAE values were –0.78 ± 0.05, –0.99 ± 0.05, and –1.12 ± 0.03 eV for groups 1, 2, and 3 in the H-rich environment, respectively. These average HAE values were significantly different from each other (p < 0.001). Accordingly, the H⁺ adsorption stability increased significantly as a function of Mg-Fe distance: group 1 < group 2 < group 3.

Fig. 3. The H⁺ adsorption energy as a function of the Fe-Mg distance under H-rich and H-poor conditions. The letters beside each symbol correspond to the configurations in Table 1. Three groups (labeled 1, 2, and 3) correspond to the number of octahedral jumps from Fe to Mg (noted by j values in Table 1). The dashed line is the linear regression line with R ² = 0.83.

The influence of the relative distance of octahedral Fe to Mg on the H⁺ adsorption strength (Fig. 3) can be ascribed to the electronegativity of these cations. The electronegativity of Fe³⁺ (1.83) is larger than Al³⁺ (1.61) and Mg²⁺ (1.31) (Pauling Reference Pauling1960). Thus, the Fe³⁺ cation tends to pull electrons toward itself and influences the polarity of bonds within these ions. The closer the Fe³⁺ cation was to the Mg²⁺ cation, the more electrons were pulled away from Mg²⁺, causing diffusion of the electron charge density of the Mg²⁺ cation and a weaker adsorption strength of the H⁺ cation. Similarly, Lavikainen et al. (Reference Lavikainen, Tanskanen, Schatz, Kasa and Pakkanen2015) observed that the distance of Mg-substitutions on montmorillonite affected the electron density around the surface oxygens. An analogous example to this phenomenon was the inductive effect observed in organic compounds, where the distance of the functional group substitutions influences the strength of the bonds (Carey and Sundberg Reference Carey and Sundberg2007). For example, if X withdraws electrons in X←COO^-, the anion XCOO^- is stabilized and the acid is strengthened (that is, the ratio of XCOO^-/XCOOH increases). Conversely, if X releases electrons in X→COO^-, the anion XCOO^- is destabilized and the acid is weakened (that is, the ratio of XCOO^-/XCOOH decreases).

The three groups of configurations with significantly different HAE identified from this DFT simulation were consistent with the observations from published cation adsorption data on Wyoming SWy-1 montmorillonite (Schulthess and Huang Reference Schulthess and Huang1990; Nolin Reference Nolin1997), which had Mg and Fe concentrations similar to those used in this DFT study. The adsorption capacities of these three types of adsorption sites on Wyoming SWy-1 montmorillonite were 43, 41, and 16% of the total CEC, respectively (Nolin Reference Nolin1997). These values agreed with the relative possibilities of the formation of configurations with j = 1, 2, and 3 as noted above, which were 43, 43, and 14%. These ratios will not be influenced by the formation of a certain group of structures. For example, if both cis- and trans- structures were formed, nine structures out of a total of 21 will be formed with j = 1, and the possibility is 43% (0.43 = 9/21). If only the cv was present, this resulted in 6 out of 14 structures or 43% (= 6/14). If only the tv was present, this resulted in 3 out of 7 structures or 43% (= 3/7). Other combinations will also lead to the same answer, such as cis-Mg or trans-Mg structures. Similarly, j = 2 has 43% probability (= 3/7), and j = 3 has 14% probability (= 1/7). Given the similarity of the chemical composition of Wyoming SWy-1 (SWy-2 or SWy-3) montmorillonite with the composition used in these DFT simulations, these percentages should also apply to the natural SWy montmorillonite. This DFT study offered an entirely new explanation of the location and proportion of these three adsorption sites, namely that all three sites could reside in the interlayer, and that the strength of these adsorption sites is a function of the direct octahedral Fe to Mg distance (or number of octahedral jumps between them).

Naturally, external edge and planar sites also exist, but their contributions to the ion-exchange reactions are small relative to the interlayer sites (Schulthess & Huang Reference Schulthess and Huang1990; McKinley et al. Reference McKinley, Zachara, Smith and Turner1995; Bradbury and Baeyens Reference Bradbury and Baeyens1997; Missana et al. Reference Missana, Benedicto, García-Gutiérrez and Alonso2014). Furthermore, if the NISE effect (Ferreira et al. Reference Ferreira, Schulthess and Giotto2011; Schulthess et al. Reference Schulthess, Taylor and Ferreira2011) can be applied to clay minerals, then the strength of adsorption sites in constrained environments (i.e. interlayer sites) can be much stronger than those of non-constrained environments (i.e. edge and planar sites). This is particularly relevant with weakly adsorbing alkali cations, such as Na⁺, that are observed to adsorb surprisingly strongly on montmorillonite at low pH values (pH < 2.5; Nolin Reference Nolin1997). Similarly, Cs⁺ cations adsorbed stronger in the constrained interlayer sites of vermiculite than the external (edge and basal) sites (Dzene et al. Reference Dzene, Tertre, Hubert and Ferrage2015).

This study can be applied to other 2:1 expansive clay minerals with different isomorphic substitutions. For example, the Fe³⁺:Mg²⁺ ratio of smectites can vary significantly. Since the electronegativity of Fe³⁺ is greater than Al³⁺, increasing the number of octahedral Fe³⁺ atoms around the Mg²⁺ atoms should also increase the electronegativity of the environment surrounding the Mg²⁺ atom. Based on the DFT results (Fig. 3), the smectite clays with higher Fe³⁺:Mg²⁺ ratios are presumed to have a weaker affinity toward the H⁺ counter cations. The presence of trace levels of other metal substitutions in montmorillonite’s octahedral sheet (e.g. Ti⁴⁺ or Mn²⁺) could also account for additional adsorption sites, albeit very small. For example, trace levels of a strong fourth site have been proposed in order to improve adsorption model fit, but the location and chemical structure of this fourth site was not specified (Tertre et al. Reference Tertre, Beaucaire, Coreau and Juery2009; Martin et al. Reference Martin, Wissocq, Benedetti and Latrille2018). Macroscopic data would have trouble distinguishing trace levels of adsorption on trace levels of octahedral substitutions, but statistical correlations of the amount of isomorphic substitution could help resolve this problem. Furthermore, rather than optimizing adsorption site densities and adsorption strengths based only on modeling adsorption data, the adsorption site densities calculated from the chemical structures of the clay minerals would help to simplify and verify the adsorption modeling of clay minerals.