## 2 Model and Methods
### Simulation Protocol  
MD simulations were conducted in Forcite module in the Material studio 8.0 simulation package and the CLAYFF force field, developed by [36]. The latter is a flexible model which describes hydrated mineral systems through nonbonded electrostatic and Lennard-Jones terms. It can accurately reproduce structural and spectroscopic properties of clays, as well as dynamical and energetic properties of clay interlayer and aqueous interfaces. The interlayer water molecules were described with the flexible version of the simple point charge (SPC) potential. The potential energy of the system was evaluated with an 8.5 A cutoff, a 0.5 A spline, and buffer width for short-range van der Waals interaction. The Ewald summation for the electrostatic interaction was calculated with a precision of  
Figure 1: Super-cell model of Na-montmorillonite (Mg octahedral are green, Al octahedral and tetrahedral purple, Si tetrahedral brown, O atoms red, H atoms white, and Na ions cyan).  
\(10^{-3}\) kcal/mol. Periodic boundary conditions were imposed on three dimensions to avoid interface effects.  
The initial conformation of hydrated MMT needs optimization of the geometry first by searching for the minimum potential energy of the system. A smart algorithm was chosen for the optimization with a maximum number of 5000 cycles, and the convergence thresholds for the specified maximum energy and displacement changes were set as \(1.0\times 10^{-4}\) kcal/mol and \(5.0\times 10^{-5}\) A, respectively. Following this, the final configuration from optimization was used as the initial configuration for the equilibration stage simulation in the NPT (isothermal-isobaric) ensemble at \(P=1\) atm and \(T=298\) K. Each system was equilibrated for 100 ps (105 steps) with a time step of 1 fs. The temperature and the pressure were controlled by a Nose-Hoover-Langevin thermostat of 1 ps and a Berendsen barostat of 0.1 ps, respectively. Finally, following the thermal equilibration configurations, a microcanonical NVE ensemble was then performed for a 1 ns production simulation, and the time step was still set as 1 fs (NVE: Each system in the ensemble has the same energy, and the number of particles and volume of each system are usually the same). The data from all 1 ns trajectories of the NVE simulation runs were recorded every 500 fs and were used to analyze the microstructure and thermodynamic properties of hydrated MMT. By comparing the relative thermodynamic energy of different systems, the steady-state water content of hydrated MMT can be determined. Furthermore, the mechanical properties of the steady-state system will be further simulated by the constant strain method. The constant strain method obtains the elastic constants via minimizing the energy of the system and deforming the supercell parameters in 12 directions in order. To improve the statistical accuracy, the simulations of mechanical properties are carried out based on the production stage trajectory files obtained above, that is, 20 frames of instant conformation are extracted from the corresponding trajectory files. The strain step is set as 10, and the maximum strain in each strain step is set to 1% according to the linear elastic assumption. The final mechanical parameters are obtained by the geometric average of the 20 frames.