Molecular dynamics modelling of Na-montmorillonite subjected to uniaxial compression and unidirectional shearing

Abstract

This paper presents systematic molecular dynamics modelling of Na-montmorillonite subjected to uniaxial compression and unidirectional shearing. An initial 3D model of a single-cell Na-montmorillonite structure is established using the Build Crystal module. The space group is C2/m, and COMPASS force fields are applied. Hydration analysis of Na-montmorillonite has been performed to validate the simulation procedures, where the number of absorbed water molecules varied with respect to the various lattice parameters. A series of uniaxial compression stress σ_zz and unidirectional shear stress τ_xy values are applied to the Na-montmorillonite structure. It is shown that the lattice parameter and hydration degree exhibit significant influence on the stress–strain relationship of Na-montmorillonite. The ultimate strain increases with increases in the lattice parameter but decreases in the number of water molecules. For saturated Na-montmorillonite, more water molecules result in a stiffer clay mineral under uniaxial compression and unidirectional shearing.

Expansive soils are distributed widely in China. These expansive soils can expand or contract to up to 15 times their original volume and are responsible for distress to some structures (Buzzi et al., 2010; Katti et al., 2015). A good understanding of the corresponding mechanical properties of expansive soils plays an important role in engineering. Continuum mechanics methods and granular mechanics methods have been applied widely to the study of the macroscopic properties of granular materials (Gökalp et al., 2009; Valter & Plötze, 2013; Wolters et al., 2015). However, most of these approaches are not able to predict accurately the nature of soils due to the variability in the particle size of soils, their irregular geometric configuration, their chemical composition, their environmental history and other factors (Anderson & Lu, 2001). At the same time, it is recognized that molecular-scale interactions play a key role in the evolution of the microstructure and macroscale properties in soils (Amarasinghe et al., 2012; Yang et al., 2019). Therefore, we must study these relevant properties of expanded clays at the molecular (atomic) scale.

The major mineral components of expansive clay minerals are montmorillonite, kaolinite and illite, of which montmorillonite is of interest due to its tendency to expand when absorbing water and to shrinkage in the absence of water (Uddin, 2008). Montmorillonite is an aluminosilicate mineral with a layered structure at the nanoscale. Each layer is composed mainly of silica–oxygen tetrahedral sheets and aluminium–oxygen octahedral sheets in a 2:1 manner. At the macroscale, saturated soil mechanics and unsaturated soil mechanics represent two parallel mechanisms for understanding soil behaviour. However, the mechanical behaviour of both saturated and unsaturated soils can be unified at the molecular scale by considering the water molecular content. After interacting with water, expansive soils demonstrate increases in their interlayer spaces (Wilson & Wilson, 2014). This swelling process impacts significantly the structural and mechanical properties of the soils (Ebrahimi et al., 2012). The water in expansive soil minerals can be divided into structural water, adsorbed water or interlayer water according to its state. The water in montmorillonite minerals is mainly in the form of interlayer water. The expansion of the soil is classified into either crystalline swelling or osmotic swelling. Crystalline swelling increases the volume of the mineral (Wilson et al., 2014), whilst osmotic swelling can increase the interlayer spacing from 20 to 130 Å and can cause macroscopic expansion of the soil (Schmidt et al., 2005; Nehdi, 2014). Osmotic swelling results in more massive overall volume increases than crystalline swelling (Swai, 2020).

Molecular simulation is a numerical simulation technique focused on the molecular (atomic) scale that can be applied to the study of the microscopic properties of matter. This method allows us not only to observe the microstructural characteristics of matter, but also to simulate the dynamics of molecules (atoms) of matter at various moments to predict the properties of materials. Thus, molecular simulation methods facilitate the study of the mineral structure and complex molecular dynamics (MD) properties of expansive soils.

MD modelling is a useful tool to understand the microevolution of materials at the atomic scale; hence, it has been used widely in chemistry, biomedicine, materials science and physics studies. MD modelling can model organic matter (e.g. proteins) successfully because the force field used is easy to determine and exhibits clear physical meaning regarding organic matter. The application of MD in geotechnical fields has become accepted by researchers only in recent years since advanced force fields have been introduced, because in the past the establishment of the clay model and selection of force fields were not satisfactory. As pointed out by Han et al. (2021), ‘understanding the physical and mechanical characteristics of clay minerals is essential in analysing geophysical seismic data, interpreting acoustic scatting measurements, and evaluating changes caused by chemical interactions between rocks and carbon dioxide, chemical contaminants, and radioactive wastes stored in underground repositories’. However, measuring the mechanical characteristics of montmorillonite is difficult due to its nanometre-scale properties and high affinity to water (Wang et al., 2001; Zhang et al., 2015).

Delville (1991) and Skipper et al. (1991) first studied clay minerals using molecular simulations in 1991. At that early stage, most researches tried to apply MD to evaluate the physical properties of clay minerals. A particular example is the study of changes in seismic-wave propagation, as hydrous phases are regarded as efficient factors in the dynamic processes of water regulation in the Earth's mantle and may play a role in earthquakes (Sato et al., 2005). In addition, the clay minerals, including montmorillonite, in the geological formation impact seismic-wave propagation (Benazzouz & Zaoui, 2012). Other investigations were conducted on the microscale properties of clay minerals, such as their expansive behaviours (Boek et al., 1995a, 1995b; Smith, 1998), the hydration and thermal properties of montmorillonite (Comin et al., 2022), homogeneous substitution (Tunega et al., 2007; Shi et al., 2013) and solute adsorption, amongst others. It should be noted that the strength and displacement of clays subjected to loading have not been studied fully at the macroscale. It is important to analyse such mechanical characteristics in terms of strength and displacement for geotechnical engineering.

Recent research has demonstrated that MD modelling provides a possible way by which to investigate the mechanical properties of clay minerals at the microscale. With respect to material mechanics, uniaxial compression and unidirectional shearing are very simple and representative loading paths to study. Understanding the stress–strain correspondence induced by loading is significant for obtaining complete information on the mechanical responses of a target material (e.g. elastic limit, yield point, ultimate stress point, fracture or breaking point, etc.). The stress–strain curve can guide the development of constitutive models, which are applicable for investigating the mechanical properties of materials. Mazo et al. (2008) and Zheng & Zaoui (2018) applied MD to analyse the elasticity of dry Na-montmorillonite in three directions. Suter et al. (2017) simulated the deformation of montmorillonite with a monolayer of water, calculated the stress–strain relationships for the four smallest models along the x and r directions and estimated the in-plane Young's modulus and bending moduli. Since then, a number of studies have investigated the impacts of hydration levels, cation-exchange capacity, anisotropy and temperature on the mechanical characteristics (e.g. elastic behaviours, in-plane Young's modulus, bulk modulus) of montmorillonite, kaolinite and pyrophyllite. However, the mechanical characteristics of the elasto-plastic stress–strain for montmorillonite subjected to loading have not been determined in the MD simulation results available at present. In addition, although most of the previous research focused on the mechanical responses of clay minerals under compression, there is a lack of research evaluating the mechanical responses of montmorillonite under shearing. However, it is understood that soils are prone to shear failure due to their structure.

Hence, in this paper, a Na-montmorillonite structure is established using the MD platform Material Studio. Uniaxial compression pressure and unidirectional shearing are applied to the micro-model of the soil clay to investigate the micro-mechanical responses of partially saturated to fully saturated Na-montmorillonite. The effects of the hydration degree and lattice parameters on the stress–strain responses of Na-montmorillonite are studied systematically. The results provide microscale information on Na-montmorillonite and could guide the constitutive modelling of expansive soils from the microscale to the macroscale.

Concepts of MD simulation

MD simulation uses Newton's Second Law of Motion to calculate the microscale information regarding atoms and molecules in various states (e.g. their position, velocity, acceleration, etc.). It can describe macroscopic information through statistical physics. The set of equations used in MD are as follows:

(1)$$m_i\displaystyle{{d^2r_i} \over {dt}} = F_i$$

where r_i and m_i denote the position and mass of the ith atom, respectively (where ith refers to the ordinal number of the number i atom), t denotes an arbitrary time instant, F_i is the force acting on atom i, and this force is calculated from the interatomic potential function U(r ₁, r ₂, …, r_N), which is also referred to as the force field representing the interactions of particles:

(2)$$F_i = {-}\displaystyle{{ {\rm \partial} U ( {r_1, \;r_2, \;\ldots r_N} ) } \over { {\rm \partial} r_i}}, \;i = 1, \;2, \;\ldots , \;N$$

(3)$$U = E_{\rm bond} + E_{\rm angle} + E_{\rm dihedral} + E_{\rm nonbond}$$

where bond energy (E _bond) is the energy of the bond stretch between two atoms (this energy will fluctuate when the actual bond length deviates from the equilibrium length or reference length), E _angle is the angle energy stored in a system of three atoms and E _dihedral is the dihedral energy stored in a system of four atoms. Both E _angle and E _dihedral are induced by the changing of the angles between three atoms and dihedrals between four atoms when the atoms move away from the equilibrium state (Katti et al., 2007). Nonbonded energy (E _nonbond) represents van der Waals and electrostatic interactions due to the attraction and repulsion between any two atoms (Katti et al., 2007).^.

Model information

Development of the Na-montmorillonite model

Building on previous research (Boek et al., 1995a, 1995b; Chang et al., 1995; Jin, 2005), an initial 3D model of the single Na-montmorillonite cell structure was developed using the Build Crystal module of Material Studio, as is shown in Fig. 1. This model belongs to the space group C2/m and symmetrical type L²PC. The lattice parameters are listed as α = 90°, β = 99°, γ = 90°, a ≈ 5.23 Å and b ≈ 9.06 Å (Fig. 2). The lattice parameter c is a variable that is influenced by the number of water molecules and the type of exchangeable cations between the layers. It should be noted that the lattice parameter c denotes the initial interlayer spacing of the Na-montmorillonite cell structure before stresses are applied. When there is no water between the layers, c ≈ 9.60 Å. To investigate the effects of various c values on the mechanical properties of montmorillonite, initial Na-montmorillonite models with c = 12.5, 15.5 or 18.5 Å were developed with reference to c of the hydration-swelling stabilization stage (Kararborni et al., 1996). The parameters and spatial coordinates of the unit cell are shown in Tables 1 & 2, respectively.

Fig. 1. Singe-crystal cell of Na-montmorillonite: (a) ball and stick; (b) polyhedron.

Fig. 2. Geometric definition of the lattice parameters.

Table 1. Parameters of the unit cell.

Table 2. Spatial coordinates of atoms and interlayered ions of the unit cell of Na-montmorillonite.

The initial lengths of single-crystal cells were 5.23 Å × 9.06 Å along the a and b directions, respectively. A superlattice of the Na-montmorillonite model consisting of eight single-crystal cells (4a × 2b × c) was developed based on the initial single-crystal cell, as is shown in Fig. 3, where the length of the superlattice of Na-montmorillonite was 2.092 nm × 1.812 nm along the a and b directions, respectively. The initial lengths along the c-direction were set as 12.5, 15.5 or 18.5 Å. The models were composed of two clay layers, each containing four-unit cells in the x-direction and two-unit cells in the y-direction. The resulting overall plane dimensions of each model were 20.92 Å × 18.12 Å.

Fig. 3. Superlattice of Na-montmorillonite: (a) ball and stick; (b) polyhedron.

It should be noted that this superlattice of Na-montmorillonite exhibits negative charge, which is in contrast to natural montmorillonite. In addition, isomorphous substitution can be observed in the tetrahedral sheets and octahedral sheets of Na-montmorillonite (Kararborni et al., 1996). Hence, the superlattice of Na-montmorillonite was made to exhibit negative charge through lattice substitution. For one layer of Na-montmorillonite, all isomorphous substitutions were carried out by replacing two Si⁴⁺ with Al³⁺ in one tetrahedral sheet and four Al³⁺ with Mg²⁺ in one octahedral sheet, as is shown in Fig. 4. The resulting negative charge is compensated for with six sodium ions between the clay sheets; the positions of these substitutions are shown in Fig. 4. The chemical formula of the Na-montmorillonite model is Na_0.75(Si_7.75Al_0.25)(Al_3.5Mg_0.5)O₂₀(OH)₄. The superlattice of Na-montmorillonite after substitution is shown in Fig. 5.

Fig. 4. Top view of the position of Na-montmorillonite (a) before and (b) after replacement.

Fig. 5. Superlattice of Na-montmorillonite after replacement: (a) ball and stick; (b) polyhedron.

The water molecule is imposed on the established superlattice of the Na-montmorillonite model using an adsorption locator to evaluate the hydration degree of Na-montmorillonite. The water molecule utilizes the extended simple point charge (SPC/E) model, and it is pre-designed to be rigid. The geometric information and parameters of the SPC/E model are shown in Fig. 6 and Table 3, respectively. The simulated models of Na-montmorillonite with one to three water layers after absorbing water are shown in Fig. 7.

Fig. 6. Water molecule. The numbers shown in this figure are explained in Table 3.

Fig. 7. Na-montmorillonite structure with water molecules: (a) 30 water molecules; (b) 58 water molecules; (c) 87 water molecules.

Table 3. Parameters of the SPC/E model for water molecules.

Structural optimization and validation of the Na-montmorillonite model

The superlattice of Na-montmorillonite was optimized, and the relative parameters are shown in Table 4. Structural optimization was performed via geometry optimization in the Forcite module in Material Studio. The smart energy minimization algorithm was applied to perform the structural optimization using the parameters referred to by Li (2016). The hydration and swelling analysis were conducted after a second structural optimization of the Na-montmorillonite with 5–96 water molecules. The root mean square of the force between atoms (RMS Force), together with the root mean square of the displacement of atoms (RMS Displacement) were given as 0.005 kal (mol Å^–1)^–1 and 0.00005 Å, respectively.

Table 4. Structural optimization parameters.

RMS = root mean square.

As shown in Figs 8–10, the hydration layer formed gradually as the number of water molecules increased. Material Studio will search for the mathematical minimum-energy conformation of the Na-montmorillonite structure. After optimization, with respect to the Na-montmorillonite absorbing 30 water molecules, one layer of water molecules was formed, as highlighted by the blue line in Fig. 8a, whilst two water layers were formed with respect to the Na-montmorillonite absorbing 31 water molecules, as shown by the blue lines in Fig. 8b. It is interesting to note that only one clay sheet was shown at random after optimization, as the structures of the clay sheet are the same at both ends, which not only reflects the Na-montmorillonite structure but also gives the minimum-energy conformation of the Na-montmorillonite. Figures 9a and 10a demonstrate that the Na-montmorillonites absorbing 58 and 87 water molecules form two and three layers of water molecules, respectively, whilst three and four layers of water molecules are formed for Na-montmorillonites absorbing 59 and 88 water molecules, respectively, as shown in Figs 9b and 10b.

Fig. 8. Na-montmorillonite adsorption structure: (a) optimized structure with 30 water molecules; (b) optimized structure with 31 water molecules.

Fig. 9. Na-montmorillonite adsorption structure: (a) optimized structure with 58 water molecules; (b) optimized structure with 59 water molecules.

Fig. 10. Na-montmorillonite adsorption structure: (a) optimized structure with 87 water molecules; (b) optimized structure with 88 water molecules.

As shown in Fig. 11, before and after optimization, the energy of the cell decreased gradually with increasing numbers of water molecules. The energy after optimization was lower than that of the original structure to a significant degree.

Fig. 11. Variation of total energy with the number of water molecules (a) before and (b) after Na-montmorillonite optimization.

Hydration swelling analysis was performed to validate the Na-montmorillonite model. The relation of the variation in lattice parameter c to the number of water molecules was analysed. As shown in Fig. 12, the lattice parameter of Na-montmorillonite increases with increasing numbers of water molecules, which corroborates the findings of Boek et al. (1995a) and Xu et al. (2016).

Fig. 12. Variation of lattice parameter c with the number of water molecules.

Loading conditions

Uniaxial loading was applied to the superlattice Na-montmorillonite model by implementing the appropriate script into Material Studio. The script was written in Perl language. The initial stress condition was achieved with reference to five rounds of zero-stress balancing. There were 100 substeps for each zero-stress balancing. There were 500 balancing substeps and 500 generation substeps. The COMPASS force field and constant-pressure, constant-temperature synthesis were applied throughout. The barostat utilized the Souza-Martins algorithm, where the box time constant of the barostat was set to 10. The Nosé–Hoover thermostat algorithm was applied. As shown in Fig. 13, the applied uniaxial compressive stresses σ_zz were 0.10, 0.12, 0.14, 0.16, 0.18 and 0.20 GPa in the vertical direction and zero in all other directions, whilst the applied unidirectional shearing stresses τ _xy were 0.10, 0.12, 0.14, 0.16 and 0.18 GPa and zero in all other directions. The successful application of these two kinds of stresses can be proven via the successful calculation of the corresponding strain tensor. The strain tensor was calculated as follows:

(4)$${\rm \varepsilon }_{ij} = \displaystyle{1 \over 2}( {h\cdot h_0^{{-}1} + h_0^{{-}T} \cdot h^T} ) -1$$

where h is the (strained) cell matrix, h ₀ is the reference cell matrix, ɛ is the strain tensor and T denotes transpose.

Fig. 13. Stress loading illustration: (a) uniaxial compression; (b) unidirectional shearing. The red arrows denote the directions in which the stresses are applied.

Three cases with respect to the effects of various lattice parameters and degrees of hydration on the micro-mechanical properties of Na-montmorillonite under uniaxial compressive stress are shown in Table 5. In addition, two scenarios for Na-montmorillonite under unidirectional shearing are illustrated in Table 6.

Table 5. Numbers of water molecules contained in Na-montmorillonite at various initial lattice parameter c values in the uniaxial compression simulations.

Table 6. Numbers of water molecules contained in Na-montmorillonite at various initial lattice parameter c values in the unidirectional shearing stress simulations.

The initial velocity at a specified temperature is often generated in a random manner in MD or a random number is introduced in the thermostat and barostat. Therefore, MD simulations are chaotic, which leads to varying results when the same simulation is conducted again on the same computer. As the initial velocity being generated in this random way did not have much impact on the simulation results of the mechanical properties of Na-montmorillonite with various degrees of hydration at various initial c values, the same model was used when running 10 identical simulations for each case under the same parameters, and the final average value of the simulation results was then analysed.

Stress–strain relationship for uniaxial compression

Effects of lattice parameter c on the stress–strain relationship of Na-montmorillonite under uniaxial compression conditions

To investigate the effect of lattice parameter c on the mechanical responses of Na-montmorillonite, the same hydration degree was set with 63 water molecules for three different initial lattice parameters (c = 12.5, 15.5 or 18.5 Å). As shown in Fig. 14, the stress increased in general with increasing strain for all conditions. For c = 12.5 Å, the stress was linearly proportional to the strain before the stress reached 0.13 GPa, then the strain increment tended to decrease with increasing stress, before finally reaching a critical state of strain at 0.15 GPa. When the lattice parameter c increased from 12.5 to 18.5 Å, the strain was significantly greater under the same stress conditions. For c = 18.5 Å, the stress–strain curve was almost linear throughout. It was assumed that greater uniaxial compression stress would be needed to determine whether a critical strain state would be achieved. This phenomenon would indicate that when the lattice parameter c is greater, the Na-montmorillonite cell would experience greater strain under the same stress conditions. A smaller c value would indicate that the critical state would be achieved rapidly.

Fig. 14. Average uniaxial compressive stress–strain relationship for Na-montmorillonite at the same degree of hydration and with various initial layers.

Effects of hydration degree on the microscopic stress–strain relationship of Na-montmorillonite under uniaxial compression conditions

The uniaxial compressive stress–strain relationships for Na-montmorillonite at various hydration levels with respect to various lattice parameter c value are shown in Fig. 15. Figure 15a,c,e illustrates the stress–strain responses of Na-montmorillonites with 44, 74 and 142 water molecules, respectively. As the MD simulations may provide random results, 10 simulation tests were performed whilst keeping all conditions the same. The average values are plotted in Fig. 15b,d,f, indicating the effects of hydration degree on the stress–strain responses. The ultimate strain approached a critical state for all cases when c = 12.5 Å, as is shown in Fig. 15a,b; under the same stress conditions, the strain increased with decreasing numbers of water molecules. The Na-montmorillonite model became softer as the number of water molecules increased. The water molecules were compressed and condensed at this stage. The more rigid the water molecules, the more resistant the Na-montmorillonite is to compressive deformation. Similar trends can be seen in Fig. 15d,f, showing that the strain decreased with increasing degree of hydration when the stress conditions were kept the same. The stiffness of Na-montmorillonite increased as the number of water molecules increased.

Fig. 15. Uniaxial compressive stress–strain relationships for Na-montmorillonite at various hydration levels: (a,b) c = 12.5 Å; (c,d) c = 15.5 Å; (e,f) c = 18.5 Å.

Figure 16 shows the stress–strain responses for saturated Na-montmorillonite under various lattice parameter c values (12.5, 15.5 or 18.5 Å). When comparing the average stress–strain relationship for initial lattice parameter c values of 12.5, 15.5 or 18.5 Å of saturated Na-montmorillonite structures at 0.1–0.2 GPa, the stress–strain relationship demonstrated a tendency towards less strain with increasing numbers of water molecules whilst under the same stress conditions until the stress reached ~0.145 GPa. At greater stresses, Na-montmorillonite that contained 63 water molecules approached a critical state, where the strain remained almost constant with increasing stress. When the stress reached 0.18 GPa, the strain of Na-montmorillonite containing 63 water molecules was very close to that of Na-montmorillonite containing 163 water molecules, and the strain of Na-montmorillonite containing 63 water molecules was less than the strain values of Na-montmorillonites containing other numbers of water molecules. The strain of Na-montmorillonite containing 163 water molecules was always less than that of Na-montmorillonite containing 113 water molecules throughout these studies. It can also be seen that the stress–strain relationship changed from non-linear to linear with increasing numbers of water molecules.

Fig. 16. Average uniaxial compressive stress–strain relationship for saturated Na-montmorillonite at various initial lattice parameter c values.

Stress–strain relationship for unidirectional shearing

Effects of hydration degree on the microscopic stress–strain relationship of Na-montmorillonite under unidirectional shearing conditions

Due to the random nature of MD simulation, 10 simulation tests under the same conditions were performed. The average results of these tests for Na-montmorillonite at various hydration levels with respect to various lattice parameter c values are shown in Fig. 17. The strain generally increased with increasing stress at the same lattice parameter c value. It was observed that the unidirectional strain followed the same tendency as that of uniaxial compression mentioned above: the unidirectional strain declined gradually with increasing numbers of water molecules. This phenomenon illustrates that the rigid interlayer water can help Na-montmorillonite to resist shearing.

Fig. 17. Unidirectional shearing stress–strain relationships for Na-montmorillonite at various hydration levels: (a) c = 15.5 Å; (b) c = 18.5 Å.

Effects of hydration degree on the microscopic stress–strain relationship of saturated Na-montmorillonite under unidirectional shearing conditions

Figure 18 illustrates the unidirectional shearing stress–strain relationships for saturated Na-montmorillonite at various hydration levels with respect to various lattice parameter c values. Figure 18 demonstrates that the greater the water content, the less the shear strain between 0.10 and 0.18 GPa unidirectional shear stress. Figure 18 also shows that the strain of Na-montmorillonite at lower lattice parameter c values was less than at greater lattice parameter c values within a stress range of 0.10–0.18 GPa. Therefore, it is predicted that the shear strain of Na-montmorillonite containing 113 water molecules would be less than that of Na-montmorillonite containing 163 water molecules at the same unidirectional shear stress and over a greater range of unidirectional shear stress.

Fig. 18. Unidirectional shear stress–strain relationship for saturated Na-montmorillonite at various lattice parameter c values.

Conclusions

This paper presented a microscale model for the investigation of the uniaxial compressive and unidirectional shearing behaviours of Na-montmorillonite via MD simulations. The structures of Na-montmorillonite were established based on the Build Crystal module, and these were validated via hydration swelling analysis in terms of the lattice parameter c varying with the number of water molecules. Uniaxial compression and unidirectional shearing were then applied to the Na-montmorillonite. The stress–strain relationships obtained from uniaxial compression and unidirectional shearing were investigated, and the following main conclusions can be drawn from the simulations:

1. (1) The stress increased with increasing strain under all conditions when subjected to uniaxial compression. The initial lattice parameter c value impacted the critical state of the stress–strain curve significantly, whilst the hydration degree demonstrated negligible effects. A greater strain can be obtained when the lattice parameter c was greater under the same stress conditions. Increasing the hydration degree in terms of increasing the number of water molecules decreased the corresponding strain under the same stress conditions. The stiffness of Na-montmorillonite increased as the number of water molecules increased. With respect to saturated Na-montmorillonite, the stress–strain relationship demonstrated a more non-linear curve with fewer water molecules.

2. (2) The strain increased with increasing stress when subjected to unidirectional shearing at the same lattice parameter c values; however, the strain reduced with increasing hydration degree, which indicated that the interlayer water can help the Na-montmorillonite to resist shearing. The critical state was not as clear for unidirectional shearing as it was for uniaxial compression. For the saturated Na-montmorillonites with 113 and 163 water molecules, corresponding to the lattice parameter c values of 15.5 and 18.5 Å, respectively, the strain decreased gradually with increasing numbers of water molecules. It was observed that the saturated Na-montmorillonite with greater numbers of water molecules was more rigid than that with fewer water molecules.