In this study, we construct the MMT model from the structure of pyrophyllite unit cell with the crystal constant a = 0.5160 nm, b = 0.8966 nm, c = 0.9347 nm, and α = 91.18°, β = 100.46°, and γ = 89.64° (Lee and Guggenheim, 1981; Skipper et al., 1995). The pyrophyllite unit cell ( Al 4 O 4 ( OH ) 4 ) ( Si 8 O 12 ) is then replicated (8 × 6 × 2) along the a, b, and c crystallographic directions, respectively to form a supercell of two TOT clay layers. The sizes of the replicated cell are Lx = 4.1280 nm, Ly = 5.3796 nm and Lz = 1.8694 nm. The constructed supercell not only decreases the subsequent excessive ordering of isomorphous substitution but also ensures that the atom in MMTs does not interact with its image in the next cell, when three-dimensional periodic boundary conditions are applied.

Based on the formed supercell, the isomorphous substitutions of octahedral Mg atom for Al atom and tetrahedral Al atom for Si atom with various compositions are introduced to establish different MMT models needed in this paper. Since one isomorphous substitution produces one negative charge (referring to the Supporting Information for details), which is referred to as the layer charge (Karnland, 2010), we introduce 48, 60, 72, 96, and 108 isomorphous substitutions to obtain the MMT models with the layer charge density (defined as the ratio between the number of layer charge to the number of unit cell, which is 96) of −0.5, −0.625, −0.75, −1.0, and −1.125 e/uc, respectively. For the MMT models with each layer charge density, three different isomorphous substitution positions, i.e., complete tetrahedral substitutions in the T sheet, entire octahedral substitutions in the O sheet and the ratio of the tetrahedral to octahedral substitutions is 1:1, are designed and the corresponding MMT models are denoted as M1, M2, and M3, respectively. All isomorphous substitutions obey Lowenstein’s substitution rule (i.e., the substitution sites cannot be adjacent to each other). We finally add the same number of Na cations as the layer charge between the layers of the MMT models to make sure that each MMT system is electrically neutral. Eventually, fifteen MMT models are established. Figure 1 shows the schematic view of the three MMT models and the snapshot of the molecular dynamic simulation supercell for the M3 Model with the layer charge of 0.75 e/uc.

In our molecular dynamic simulation, the Parinello−Rahman (PR) strain fluctuation formula (Ray and Rahman, 1985; Cui et al., 2007; Carrier et al., 2014) is used to calculate the anisotropic elastic parameters in the constant-pressure and constant-temperature (NPT) ensemble. The PR strain fluctuation formula is given as C i j k l = K B T 〈 V 〉 [ 〈 ε i j ε k l 〉 − 〈 ε i j 〉 〈 ε k l 〉 ] − 1 (1) where K B is the Boltzmann constant, equal to 1.38065 × 10^–23 J/K, T is the temperature in Kelvin, V is the volume in cubic meters, ε i j is the strain tensor, and   〈 ▪ 〉 denotes the ensemble average. The subscripts i, j, k and l run from one to three and represent the three dimensions in the Cartesian coordinates. The indexes one and two correspond to the x- and y-axes, respectively, which are parallel to the MMT layer, and the index three corresponds to the z-axis, which is perpendicular to the MMT layer. The strain tensor ε i j can be determined by ε i j = 1 2 [ 〈 h 〉 i k − T h k l T h l m h m j − 1 − δ i j ] (2) where δ i j is the Kronecker tensor, and h is the scaling matrix, which can be obtained through molecular dynamic simulation.

In this work, the LAMMPS code (Plimpton, 1995) is employed to perform the molecular dynamic simulation to obtain the scaling matrix h. In the simulation, a force field is required to describe the interactions of the atoms in the MMT, and we use the CLAYFF force field (Cygan et al., 2004) for this purpose (the interaction parameters and charge of the relevant atoms are shown in the Supporting Information). The CLAYFF force field consists solely of the electrostatic terms, the Lennard-Jones terms and the O-H bond constrained by a harmonic bond stretch potential and is widely used in the simulations of clay minerals (Ebrahimi et al., 2012; Carrier et al., 2014). In our molecular dynamic simulations, we set the time step to be 1 fs in the NPT ensemble under stratigraphic conditions (i.e., p = 10 MPa, and T = 323 K). The Langevin dynamics and Nosé−Hoover barostat are used to control the pressure and temperature, respectively. The cutoff distance of the nonbonded interaction is set to be 10.0 Å, and the long-range electrostatic interactions are calculated employing the particle-particle particle mesh (PPPM) method with an accuracy of 10^–4. The system is equilibrated for 2 ns. The box size (lx, ly, lz), tilt factors (xy, xz, yz) and volume of box (V) are then sampled every ten steps for 8 ns. After that, the sample data (lx, ly, lz, xy, xz, yz) are sorted to form the scaling matrix h in the form of h = [ l x     x y       x z 0                 l y         y z 0                 0               l z ] (3)

Then the resulting scaling matrix h is substituted into Eq. 2 to obtain the strain tensor ε i j , which are further integrated into Eq. 1 with V and T to determine the elastic tensors C i j k l of the MMTs.

When we show the elastic tensors, the Voigt notation is used to represent the calculated fourth-order stiffness tensor C i j k l by the second-order tensors C i j . The indexes change as follows: 11 → 1 , 22 → 2 , 33 → 3 , 23 → 4 , 13 → 5 , and 12 → 6 . The indexes 1, 2, and three correspond to the compressions in the x, y, and z directions, respectively, and the indexes 4, 5, and six correspond to the shears in the yz-, xz-, and xy-planes, respectively. Although 36 coefficients of the stiffness tensor are computed, we only report the nine coefficients C ₁₁, C ₂₂, C ₁₂, C ₃₃, C ₁₃, C ₂₃, C ₄₄, C ₅₅, and C ₆₆ because of the transverse isotropic characteristics of MMT (Ebrahimi et al., 2012; Carrier et al., 2014). According to the deformation modes involving the crystalline structure of the clay layer and interlayer space, the nine data are divided into three groups: the in-plane compressional coefficients C ₁₁, C ₂₂, and C ₁₂, the out-of-plane compressional coefficients C ₃₃, C ₁₃, C ₂₃, and the shear coefficients C ₄₄, C ₅₅, and C ₆₆.

The variations of the in-plane compressional coefficients C ₁₁, C ₂₂, and C ₁₂ with different layer charge densities for the three models with different layer charge distributions are shown in Figure 2. It is interesting that all the three in-plane compressional coefficients decrease approximately linearly with the increasing layer charge density in the three models. The reduction in the three coefficients indicates that the compressional resistance of the three models decreases with the increasing layer charge density. Specifically, the decreasing C ₁₁ and C ₂₂ represents the reduction in the compressional resistance in the x-axis and y-axis, respectively, and the declining C ₁₂ stands for the decrease of the compressional resistance in the y-direction when a given positive pressure is applied to the x-axis.

Since C ₁₁ and C ₂₂ reflect the strength of the intermolecular interactions of the MMT molecules in the x-axis and y-axis, respectively, and C ₁₂ reflects the strength of the intermolecular interactions in the xy-plane, the decrease of C ₁₁, C _22, and C ₁₂ can be explained by the weakening of the intermolecular interactions of the MMT in these directions. The intermolecular interaction is inversely proportional to the atomic distance (Cygan et al., 2004; Davarcioglu, 2011), and thus the decrease of the three coefficients with the increasing layer charge density may be related to the increase in the atomic distance of the MMT. Figure 3 shows the atomic distances (distance between the replaced atom or replacement atom and the adjacent oxygen atom) before and after the substitutions, and the results clearly demonstrate that the atomic distances in the MMT layer (composed of three sheets) after the substitutions are increasing, no matter whether the tetrahedral substitutions or the octahedral substitutions are introduced into the MMT. The increase in the atomic distance in the MMT layer may lead to the weakening of the intermolecular interactions within the MMT layer in all directions.

To further support the explained variation of the interaction strength, we also calculate the interaction energy (the difference between the molecules’ combined energy and all of their isolated energies) between the replaced atom or replacement atom and the adjacent oxygen atom (i.e., the forming atom pair) as shown in Figure 4. It can be seen that the interaction energy between the atom pair increases after the tetrahedral substitutions or the octahedral substitutions, directly indicating that the interaction within the MMT layer is weakened (Šponer et al., 1999; Lovelock, 2017). Considering the in-plane compressional coefficients are related to the intermolecular interactions in the x- and y-axes of the MMT layer (Carrier et al., 2014), the increase in the atomic distances will reduce these intermolecular interactions in these directions and hence the in-line compressional elastic coefficient C ₁₁, C ₂₂, and C ₁₂.

It is worth noting that although C ₁₁, C ₂₂, and C ₁₂ of all the models are decreasing, we find the values and the variations of C ₁₁, C ₂₂, and C ₁₂ in the three models are different. First, it is evident that C ₁₁, C ₂₂, and C ₁₂ in the M2 Model are larger than those in the other two models when the layer charge densities are the same. This means that C ₁₁, C ₂₂, and C ₁₂ of the MMT with the layer charges distributed in the O sheet are greater than those of the MMT with the layer charges distributed in the two T sheets and in both the O sheet and the T sheets. This result can be closely related to the structure of the MMT. As shown in Figure 1D, the MMT consists of a sandwich of one O sheet between two T sheets, called a TOT structure. Since the O sheet is constrained by the two T sheets and there are interactions among the sheets, the increase in the atomic distance and in the interaction energy between the atomic pairs in the O sheet after the octahedral substitution is weakened. As a result, the intermolecular interactions in the xy-plane are less reduced, leading to a greater C ₁₁, C ₂₂, and C ₁₂ in the M2 model.

In addition to the difference in the values of the three in-plane compressional coefficients, it is also interesting to note that the changing rates of the compressional coefficient with layer charge densities in the three models are different. Table 2 tabulates the slopes of C ₁₁, C ₂₂, and C ₁₂ after linear fitting for the three models. By comparing the slopes of the same parameter in the three models, we find that the M1 model has the largest slopes of C ₁₁, C ₂₂, and C ₁₂, followed by the M3 model and the M2 model. This result implies that the layer charge density has a greater effect on C ₁₁, C ₂₂, and C ₁₂ of MMT with the layer charges distributed in the T sheet. The difference in the slopes of the in-plane compressional coefficients in the three models can again be explained through the structure of the MMT, where the TOT architecture will confine the expansion of the O sheet, and hence the octahedral substitution has less impaction on the elastic properties, resulting in a gentler slope of the in-plane compressional coefficients with varying layer charge density.

Figure 5 shows the variations of the out-of-plane compressional coefficients C ₃₃, C ₁₃, and C ₂₃ with the increase of layer charge density for the three models with different layer charge distributions. Unlike the in-plane compressional coefficients that show consistent downward trends in the three models, the out-of-plane compressional coefficients show very different variations with layer charge density in the three models. Therefore, to more conveniently understand the influences of the layer charge distribution on the out-of-plane compressional coefficients, we choose to analyze the same elastic coefficient in the different models.

We first compare the variations of C ₃₃ with layer charge density in the three models that represents the compressional resistance of the MMT in the z-axis. As shown in Figure 5A, C ₃₃ of the M1 Model increases slightly and then decreases linearly with the increase of layer charge density. The initial increasing and then decreasing C ₃₃ indicates that the increasing layer charge density first increases and then reduces the compressional resistance of the M1 Model in the z-axis. In a different manner, C ₃₃ of the M2 Model consistently increases with increasing layer charge density, suggesting that the compressional resistance of the M2 Model in the z-axis is increasing with the increasing layer charge density. For the M3 Model, C ₃₃ doesn’t vary much as the layer charge density increases, implying that the layer charge density is not affecting significantly the compressional resistance of the M3 Model in the z-axis.

Because C ₃₃ reflects the strength of the intermolecular interactions of the MMT in the z-axis, the different trends of C ₃₃ with the increase of layer charge density in the different models can be explained by the difference in the z-axis intermolecular interaction strength of the MMT caused by the layer charge distributions. The intermolecular interactions of the MMT in the z-axis involve the intermolecular interactions in the z-axis both within the MMT layer (each layer is composed of three sheets) and between the neighboring MMT layers because the MMT contains the MMT layers and the interlayers in the z-axis as shown in Figure 1D. We have employed atomic distance inversely to represent the strength of the intermolecular interactions in the previous section, however, the atomic distance can only reflect the intermolecular strength of the MMT layer in the x- and y-axes, but will not fully represent the intermolecular strength in the z-axis. This is because the z-axis intermolecular interaction includes not only the z-axis intermolecular interaction within each MMT layer but also that between the MMT layers (i.e., the interlayers). To characterize the strength of the intermolecular interactions of the MMT in the z-axis, we will use the interaction energy of the MMT layer and the interlayer. However, the interaction energy will not directly quantify the interaction strength in a particular direction because it contains the information of the intermolecular interactions in all directions. Therefore, we combine the interaction energy with the calculated atomic distance (as shown in Figure 3) that reflects the intermolecular interactions in the x- and y-axes for the better characterization of the intermolecular strength in the z-axis.

Figure 6 shows the simulated interaction energy in the MMT layer and the interlayer as a function of layer charge density for the three models. In the three models, the interaction energy in the interlayer all decreases with the increasing layer charge density, indicating enhanced interlayer interactions (Shiu and Tsai, 2014; Wang et al., 2018). The decreasing interlayer interaction energy with increasing layer charge density can be explained by the fact that with increasing layer charge density there are more compensating Na cations in the interlayer that bind stronger the MMT layers (Zhong et al., 2021), resulting in lower interlayer interaction energy. On the contrary, the interaction energy in the MMT layer varies differently with the increase of layer charge density in the three models. The layer interaction energy of the MMT increases remarkably with the increasing layer charge density in the M1 model, and the increase becomes less significant in the M3 model, indicating the strength of the intermolecular interactions within the MMT layer are reducing significantly and gently in the M1 Model and M3 Model, respectively. In a different way, the interaction energy of the MMT layer reduces with the increasing layer charge density in the M2 Model, suggesting an increasing interaction strength of the intermolecular within the MMT layer with the increasing layer charge density.

The decreasing interaction strength of the MMT layer (the increasing layer interaction energy) in the M1 Model can be understandable, because the isomorphous substitutions in the T sheet will increase the atomic distance (as shown in Figure 3A) and hence reduce the interaction strength, and therefore more isomorphous substitutions with greater layer charge density will result in the more significant reduction in the interaction strength. However, the increasing atomic distance due to the isomorphous substitutions cannot explain the obtained increasing interaction strength (decreasing layer interaction energy) in the M2 Model. Although the isomorphous substitutions in the O sheet (M2 Model) will increase the atomic distance, the increase occurs only in the O sheet, and with the increasing atomic distance in the O sheet, the distance of the atoms between the O and T sheets is reduced. As a result, the interaction strength between the O sheet and the T sheets gets enhanced, leading to an increasing layer interaction strength. Accordingly, the obtained increasing layer interaction strength with layer charge density in the M2 Model implies that the increasing interaction strength between the O sheet and the T sheets overwhelms the decreasing interaction strength in the O sheet. Similarly, the gently decreasing interaction strength of the MMT layer with layer charge density in the M3 Model indicates that the effects of the increasing interaction strength between the O sheet and the T sheets are less significant than the decreasing interaction strength in the O and T sheets.

After analyzing the impact of the layer charge density on the two interactions and on the interactions in the x- and y-axes (which can be represented by the in-plane compressional coefficients presented above) in the three models, the variations of the calculated C ₃₃ in the three models that reflect the interaction strength in the z-axis can be more understandable. For the M2 Model, the interaction strength in the x- and y-axes decreases with increasing layer charge density, however, the total interaction strength (i.e., the interaction strength in both the MMT layer and the interlayer) increases with layer charge density, and accordingly the z-axis interaction strength (i.e., C ₃₃) must increase with the increasing layer charge density. Unlike C ₃₃ in the M2 Models, C ₃₃ in the M1 Model and M3 Model can not be directly explained by the total interaction strength and the interaction strength in the x- and y-axes, but can reflect the competition between the interaction strength of the layer and interlayer in the z-axis (i.e., the difference between the layer and interlayer interaction strength and their contribution from the x- and y-axes, respectively).

For the M1 Model, because the total interlayer interaction strength increases while the interaction strength in the x- and y-axes decreases with increasing layer charge density, the z-axis interlayer interaction strength must increase with the increasing layer charge density. On the other hand, the variation of the layer interaction strength in the z-axis with the increasing layer charge density cannot be determined because both the layer interaction strength and the interaction strength in the x- and y-axes decrease with layer charge density. When the layer charge density is low, it can be difficult to deduce the variation of the layer interaction strength in the z-axis from the increasing C ₃₃, because the interlayer interaction strength in the z-axis also increases with increasing layer charge density. However, when the layer charge density increases from 0.625 e/uc, the decreasing C ₃₃ indicates the layer interaction strength in the z-axis also decreases and the decreasing of the layer interaction strength in the z-axis is playing a more significant role than the increasing interlayer interaction strength in the z-axis.

Since the trends of the interlayer interaction strength and the layer interaction strength in the M3 Model are similar to those in the M1 Model, it is reasonable that the z-axis interlayer interaction strength increases (because the interaction strength in the x- and y-axes decreases) and the layer interaction strength in the z-axis can not be determined (because both the layer interaction strength and the interaction strength in the x- and y-axes decrease) in the M3 Model. However, according to the fact that the C ₃₃ of the M3 Model doesn’t vary much with the increase of the layer charge density and the z-axis interlayer interaction strength increases, we can infer that the layer interaction strength in the z-axis decreases and the decreasing interlayer interaction in the z-axis balances with the increasing layer interaction in the z-axis.

Having analyzed the variations of C ₃₃ with layer charge density and provided a plausible explanation through the interaction energy, we proceed to present and interpret the variations of C ₁₃ and C ₂₃ with layer charge density in the three models with different layer charge distributions. Interestingly, although the values are much smaller, C ₁₃ and C ₂₃ show very similar trends with C ₃₃ in each model, and this indicates that the interaction strength in the z-axis that determines C ₃₃, may also be the main factor influencing C ₁₃ and C ₂₃. Therefore, it can be reasonable that the difference in variations of the out-of-plane compressional coefficients with layer charge density in the three models with different layer charge distributions is controlled mainly by the layer interaction strength in the z-axis.

Figure 7 shows the variations of the in-plane shear coefficient C ₆₆ and the out-of-plane shear coefficients C ₄₄ and C ₅₅ with the increasing layer charge density. We can see that C ₆₆ of the three models keeps almost constant at around 80 GPa as the layer charge density increases, indicating the shear deformation resistance in the xy-plane is independent of the density and distribution of the layer charges. Unlike the flat in-plane shear coefficient C ₆₆, the out-of-plane shear parameters C ₄₄ and C ₅₅ are fluctuant with the increasing layer charge density in all the three models, and their values in each model are approximately the same. This suggests the shear deformation resistance in the yz-plane and xz-plane is not significantly affected by the layer charge density and the layer charge distribution. Since the shear coefficients are mainly reflecting the torsion interaction of the bonding atoms at the molecular level (Sridharan, 2014; Yang and Xu, 2021), the unchanging shear coefficients indicate that the torsion interaction of the bonding atoms is independent of the varying layer charge density and layer charge distribution.

Although the layer charge density and layer charge distribution are not impacting greatly the shear coefficients, it is found that the out-of-plane shear elastic coefficients C ₄₄ and C ₅₅ are considerably smaller than the in-plane shear elastic constant C ₆₆. This result indicates that the shear deformation is more likely to occur in the yz-plane and xz-plane than in the xy-plane of the MMTs. This may be related to the interlayer structure of the MMTs as shown in Figure 1D. Because the interlayer structure is perpendicular to the z-axis and the interlayer interactions are much weaker than the intermolecular interactions within the MMT layer (Zhang et al., 2018; Zhao et al., 2021), the yz- and xz-planes are easier to form slip deformation than the xy-plane when a shear force is applied, resulting in the smaller shear coefficients C ₄₄ and C ₅₅ than C ₆₆.