*3.1. Swelling and Thermodynamics Analysis*

To obtain the water content corresponding to the thermodynamic steady-state in the hydration process of MMT, the strategy of simulating 23 different water content systems was adopted. The interlayer space of the MMT lamellae increases with the increase of water content, and this swelling process can be described by using the parameter of basal spacing. The basal spacing can be obtained from the view of water content by averaging the system volume during the production stage of simulation and defined as in Equation (13):

$$d\_{001} = \frac{\langle V \rangle}{\langle a \rangle \langle b \rangle \sin \langle a \rangle} \tag{13}$$

where <*V*> denotes the statistically averaged volume, and <*a*> <*b*> <*α*> are the statistically averaged parameters of the simulation supercell.

Figure 2 illustrates the variation of basal spacing of different hydrated montmorillonite systems as a function of water content from this study and available results from the literature. In the early stage of crystalline swelling, the basal spacing increases from step to step, while the step characteristic (plateau) disappears, and shows approximately linear variation in the following stage of osmotic swelling. The swelling curves of MMT with different compensation cations are different. The basal spacing of Cs-MMT is greater than Na- and Ca-MMT under the same water content. However, there is almost no difference between the latter two. This is mainly due to the fact that Cs+ ions are a heavy metal ion and their ion radius is about 1.7 times that of Na+ and Ca2+ ions, while the latter two's ion radii are almost the same. The simulated results can well penetrate the distribution range of other available measurements.

**Figure 2.** The simulated and experimental basal spacing of montmorillonite as a function of water content (*g* denotes molecular mass).

It has been shown that the energy contribution to free energy dominates clay swelling and the entropy term only played a minor role. The local minima of the swelling energy curve correspond to the stable hydration state. Hence, the relative stabilities of different states can be determined by immersion energy and hydration energy. The immersion energy is defined as in Equation (14):

$$Q = \left< E(N) \right> - \left< E(N^0) \right> - \left( N - N^0 \right) E\_{\text{bulk}} \tag{14}$$

Here <*E*(*N*)> is the average potential energy of the hydrate with water content *N*, and <*E*(*N0*)> is the average energy of the referential hydration state (the state with the highest water content, 49%, was selected). *Ebulk* is the mean interaction potential of the bulk flexible SPC water; *Ebulk* = −41.5392 kJ/mol. The immersion energy *Q* is the energy released when MMT with water content *N* transforms into the referenced state with water content *N*<sup>0</sup> by adsorbing water from bulk water.

The hydration energy is defined as in Equation (15)

$$
\Delta E = \frac{\langle E(N) \rangle - \langle E(0) \rangle}{N} \tag{15}
$$

Here <*E*(0)> is the average potential energy of the dry MMT. The hydration energy Δ*E* evaluates the energy change associated with water uptake by the dry clay.

The simulated hydration energy and immersion energy curves of MMT with different compensation cations are presented in Figure 3. For Na- and Cs-MMT the local minima occur when the interlayer water molecular number is 40 and 80 in the hydration energy curves, while Ca-MMT has only a global minimum at low water content. Correspondingly, in the immersion energy curves of Na- and Cs-MMT, local minima occur when the interlayer water molecular numbers are also 40 and 80. The global minima of Cs-MMT is found at *n* = 80, while Na-MMT global minima is located at *n* = 160, as shown by the enlarged illustration in Figure 3. Meanwhile, the variation of the Ca-MMT immersion energy curve decreases monotonically with the increase of water content, and the global minima occurs also at *n* = 160. Then, since the local minima of the swelling energy corresponds to the thermodynamic stable states, it can be deduced that the Na- and Cs-MMT will form the single- and double-layer stable hydrated states when the number of interlayer water molecules reach 40 and 80 respectively. The double-layer hydrate for Cs-MMT is in the most stable state. Thus, the Cs+ ions can inhibit swelling. The global minima of Na-MMT is located at *n* = 160. After forming the first and second hydration layer, Na-MMT will continue expansion to the osmotic swelling stage, namely that the Na+ ions promote swelling in contrast with Cs+. For Ca-MMT, its role in promoting expansion is even more

pronounced. The intermediate state of the first or second hydration layer will not form during the swelling process, and the osmotic swelling stage will rapidly and directly develop. The nanoscale elastic properties of these stable state water content systems will be further simulated in the following sections.

**Figure 3.** Hydration energy and immersion energy of montmorillonite as a function of increasing interlayer water molecule number.

## *3.2. Nanoscale Elastic Properties*

To fully demonstrate the yield behavior of the MMT mineral, the maximum strain of each direction is expanded to 10% and the corresponding calculation step is increased to 100 steps. The stress–strain curves of the typical MMT system obtained from simulation are shown in Figure 4, where the interlayer compensation cation is Ca2+ ion, and the interlayer water number *n* is 80, corresponding to bilayer hydration.

**Figure 4.** Typical stress–strain behavior for hydrated Ca-montmorillonite.

The micromechanical properties of hydrated MMT show obvious anisotropy. The in-plane uniaxial tensile and compressive strength of MMT is approximately equal (for *σ*<sup>1</sup> and *σ*2), while the tensile strength is slightly less than the compressive strength. When the tensile strain reaches about 7%, the material yields. Meanwhile, the compressive stress– strain relationship shows a good linear relationship within the maximum compressive strain (10%). The tensile and compressive strength (*σ*3) of MMT orthogonal to the mineral plane is much lower than that in-plane strength (*σ*<sup>1</sup> and *σ*2). This is mainly due to the appearance of the interlayer water in MMT lamellae, and the bonding strength between the MMT crystal and water molecule through non-bonding is much lower than that of covalent or ionic bonding in the mineral plane. In addition, the compression of MMT in the Z direction shows obvious non-linear characteristics and hardening under high pressure, which is different from the linear characteristics in the crystal plane.

The hydration molecular layer that mixed within MMT lamellae has a certain shear strength (τ4 and τ5) and also shows anisotropy along with different directions of the MMT crystal plane. The shear stiffness of the interlayer water along the Y direction (τ5) is slightly larger than that in the X direction (τ4), while its shear strength is slightly lower. The shear capability of MMT lamellae is very high (τ6) and shows good linear characteristics in the simulated strain ranges.

By comparing the quantitative values of in-plane strength, it was found that the inplane compressive strength of hydrated MMT is the largest, followed by tensile strength and shear strength. Similarly, this variation also applies to the strengths that are perpendicular to the mineral plane.

All of the stress responses show good linearity in the range of 1% small strain, and this is illustrated in the lower right of each block in Figure 4. This also shows that the linear elastic constants calculated by the small strain of 1% are reliable. The stiffness tensor and linear elastic constants of typical MMT in the 1% strain range corresponding to Figure 4 are listed in Table 1.


**Table 1.** Typical stiffness coefficients and elastic constants of hydrated Ca-MM (GPa).

The simulated stiffness tensor components in Table 1 are in good agreement with simulations and measurements in the literature. The in-plane components *C*<sup>11</sup> and *C*<sup>22</sup> of the MMT with water bilayer at 300 K are 189 GPa and 191 GPa, respectively, from [25], which are almost equal to the results of this study that *C*<sup>11</sup> = 189 GPa and *C*<sup>22</sup> = 212 GPa. Furthermore, there are differences in the in-plane Young's moduli *E*<sup>1</sup> and *E*<sup>2</sup> which are derived from the stiffness tensor, and these differences are induced by the different topologies of trivalent cations in dioctahedral structure of MMT, as shown in Figure 5. However, the in-plane Young's moduli (*E*<sup>1</sup> and *E*2) are all much higher than the Young's modulus perpendicular to the plane (*E*3). On the contrary, the in-plane shear moduli (*G*<sup>23</sup> and *G*31) are lower than the shear modulus perpendicular to the plane (*G*12). The larger difference of Poisson's ratio in different directions once again shows the anisotropy of the MMT structure. The boning

strength in the mineral plane is much higher than that of non-bonding strength between mineral lamellae.

**Figure 5.** Different topologies of Al atom in the dioctahedral structure of montmorillonite.

The evolution of elastic properties of hydrated MMT under different stable state water content will be further illustrated by analyzing the changes of some components in the stiffness tensor. To validate the study, the simulation results from [29] and UPV test results from [17] are adopted for comparison due to the high correlation with the work in this study. The authors in [29] gave the variation of stiffness components during the hydration process of Na-Wyoming MMT with a water content of 0–40%. The authors in [17] obtained the range of elastic stiffness constants of smectite minerals on the microscale based on UPV test results of clay composition in shale.

For compression of clay under high pressure, the most concerning component of the stiffness tensor is *C*33, that is, the ability to resist deformation in the direction of perpendicular to the mineral plane. Figure 6 shows the variation of *C*<sup>33</sup> obtained upon hydration by [29], Na-, Cs-, Ca- MMT in this study under different stable state conditions, and the test range results of [17]. Since the stiffness tensor elements in this study are obtained by a statistical average of 20 frame trajectory files, the statistical standard deviation of data is also given in Figure 6.

**Figure 6.** The *C*<sup>33</sup> elastic constant as a function of basal layer spacing.

Ebrahimi et al. (2012) revealed that *C*<sup>33</sup> decreased sharply from dry MMT to the formation of the first hydrated layer between lamellae, then increased to a local maximum with the formation of the stable hydrated layer, and finally maintained a constant until the basic spacing d > 15.6 Å [29]. The values of *C*<sup>33</sup> obtained in this study are in good agreement with the results of [29], apart from the first hydrated layer. The *C*<sup>33</sup> of Ca-MMT is slightly higher than that of Na-MMT, while the *C*<sup>33</sup> of Cs-MMT is the lowest among all the same hydrated states. Apart from the dry MMT, all of the simulated *C*<sup>33</sup> fall well into the given range of UPV measured by [17] (*C*<sup>33</sup> = 11.5~30.4 GPa).

The in-plane stiffness tensor components *C*11, *C*22, and *C*12, which reflect the Poisson effect in the plane, are shown in Figure 7 with basal spacing. *C*11, *C*22, and *C*<sup>12</sup> all decrease nonlinearly with the increase of water content. These components are mainly controlled by the ionic or covalent bonding strength of mineral atoms and the geometric conformation of the hydrated MMT system. In addition, the simulation results of *C*<sup>11</sup> and *C*<sup>12</sup> are in good agreement with those of [29], but *C*<sup>12</sup> is slightly lower. The reason may lie in the difference between our simulation structure and the octahedral isomorphic substitutions of [29]. However, the measured range of UPV given by [17] (*C*<sup>11</sup> = 13.8~46.1 GPa and *C*<sup>12</sup> = 6.9~17.8 GPa) is much lower than the simulated results. This is because the MMT mineral lamellae are very thin, and the relative slip between lamellae easily occurs in the UPV test. However, the simulated results in [25] (*C*<sup>11</sup> = 189~231 GPa) and test results by Brillouin scattering on mica minerals in [26] (*C*<sup>11</sup> = 181 GPa and *C*<sup>22</sup> = 178 GPa) both show that the results of this study are credible.

**Figure 7.** The *C*11, *C*22, and *C*<sup>12</sup> elastic constants as a function of basal layer spacing.

For the shear of clay under high pressure, the most concerned is the shear resistance of limited interlayer water molecules in the hydrated MMT system, i.e., the in-plane shear stiffness tensor components *C*<sup>44</sup> and *C*55. Their variations with basal spacing are shown in Figure 8; *C*<sup>44</sup> and *C*<sup>55</sup> fluctuate with the increase of water content which is like the variation of *C*33. Therefore, the mechanism of shear strength of interlayer water is closely related to the formation of H-bonds and atom-free volume. The shear strength between dry MMT lamellae is much larger than that of hydrated MMT lamellae, that is to say, the addition of water plays a lubricating role in interlayer friction, which makes *C*<sup>44</sup> and *C*<sup>55</sup> decrease sharply. However, these limited interlayer water molecules do not completely show the zero shear strength characteristics of free water. As water content increases, the in-plane shear strength increases to a local maximum when the system approaches the thermodynamic stable state and minimum atom free volume occurs. In the range of simulated water content, *C*<sup>44</sup> and *C*<sup>55</sup> eventually tend to be constant. Furthermore, the results in Figure 8 shows that the components *C*<sup>44</sup> and *C*<sup>55</sup> are in good agreement with those obtained by [29], and the in-plane shear strength between different MMT systems is close except for the dry MMT system. Meanwhile, the simulated results also fall well into the UPV measurement range of [17] (*C*<sup>44</sup> = 2.9~8.9 GPa).

**Figure 8.** The *C*<sup>44</sup> and *C*<sup>55</sup> elastic constants as a function of basal layer spacing.

The shear of hydrated MMT along the direction perpendicular to the mineral plane is mainly controlled by the relative torsion of bonded atoms. The variations of *C*<sup>66</sup> with basal spacing are shown in Figure 9. With the increase of interlayer water content, *C*<sup>66</sup> decreases continuously, which is mainly due to the increase of interlayer water contribution in the weak shear zone. Considering that the ratio of width to height of a true MMT lamella is very large (about 100:1 to 1000:1), shear along the direction of orthogonal to the mineral plane is almost impossible, so the variations of *C*<sup>66</sup> will not be discussed here.

**Figure 9.** The *C*<sup>66</sup> elastic constant as a function of basal layer spacing.

Young's moduli of the different MMT systems under various hydration can be obtained by further conversion of the stiffness tensor, and the variation of Young's moduli with basal spacing is shown in Figure 10. The in-plane Young's moduli *E*<sup>1</sup> and *E*<sup>2</sup> are consistent with *C*<sup>11</sup> and *C*<sup>12</sup> that decrease nonlinearly with the increase of water content. For the Young's modulus of *E*<sup>3</sup> which is perpendicular to the mineral plane, its value for the dry MMT system is 2–3 times that of the hydrated system for Na-MMT and Ca-MMT. *E*<sup>3</sup> of the hydrated MMT system also decreases nonlinearly with the increase of water content. Interestingly, for Cs-MMT, the Young's modulus, E3, of different hydrated MMT systems, even dry MMT systems, is almost the same. Preliminary analysis shows that this may be related to the larger radius of Cs+ ions. For the dry MMT system, the interlayer spacing is large, which induces a low compression strength. After hydration, the Cs+ ions always form an inner-sphere complex, which has little effect on the strength, thus making all Cs-MMT systems of E3 almost the same within the range of simulated water content. However, this

interesting phenomenon has not yet been reported in the literature, so the above corollary needs further study in the future.

**Figure 10.** Young modulus as a function of basal layer spacing.

From the above comparisons, molecular microscopic testing techniques, such as UPV, can only qualitatively analyze the strength of clay minerals and the test results are not accurate enough. Molecular dynamics simulation provides a good way to understand the properties of hydrated clay minerals on an atomic scale. In particular, unconstrained system techniques and thermodynamic stable-state water content analysis techniques are used to make the calculation results of molecular dynamics more reliable.
