Three new structural models of montmorillonite with differently distributed Al/Si and Mg/Al substitutions in the tetrahedral and octahedral clay layers are systematically developed and studied by means of MD simulations to quantify the possible effects of such substitutional disorder on the swelling behavior, the interlayer structure, and mobility of aqueous species. A very wide range of water content, from 0 to 700 mg water /g clay is explored to derive the swelling properties of Cs−montmorillonite. The determined layer spacing does not differ much depending on the clay model. However, at low water contents up to 1-layer hydrate (∼100 mg water /g clay ) the variation of specific locations of the tetrahedral and octahedral substitutions in the two TOT clay layers slightly but noticeably affects the total hydration energy of the system. Using atom−atom radial distribution functions and the respective atomic coordination numbers we have identified for the three clay models not only the previously observed binding sites for Cs + on the clay surface but also new ones that are correlated with the position of tetrahedral substitution in the structure. The mobility of Cs + ions and H 2 O diffusion coefficients, as expected, gradually increase both with increasing water content and with increasing distance from the clay surface, but they still remain 2 to 4 times lower than the corresponding bulk values. Only small differences were observed between the three Cs− montmorillonite models, but these differences are predicted to increase in the case of higher charge density of the clay layers and/or interlayer cations.