Laterite is a reddish-brown clayey soil enriched with iron, aluminium oxides and clay minerals. It is formed through the lateritization of parent rocks in humid tropical climates (Rashid et al., Reference Rashid, Tabatabaei, Horpibulsuk, Mohd Yunus and Hassan2019). Owing to its distinctive formation environment, laterite exhibits distinct engineering properties compared to conventional soils. While it demonstrates poor physical characteristics (e.g. high liquid and plastic limits, elevated porosity, high water content and low density), it possesses favourable mechanical properties, including high strength and low compressibility. This gives rise to a seemingly paradoxical combination in which its physical indicators contradict its mechanical parameters. In recent years, the distinctive properties of laterite have made it a valuable industrial material, finding widespread application in construction, environmental engineering and materials science (Wang et al., Reference Wang, Meng, Jiang, Sarwar, Fu and Yang2022). Its low permeability renders it an economical and effective impervious material, being widely used in landfill capping, tailings dam rehabilitation and as clay liners in earth-rock dams (He et al., Reference He, Zhang, Chen and Long2022). Additionally, laterite’s high adsorption capacity and strong ion-exchange capabilities enable its use as a soil amendment for heavy-metal remediation (Kong & Luo, Reference Kong and Luo1993; Mengue et al., Reference Mengue, Mroueh, Lancelot and Eko2017). However, the presence of hydrophilic minerals in laterite makes it susceptible to significant volume changes in response to wet–dry cycles, leading to pronounced swelling and shrinkage behaviour (Zeng et al., Reference Zeng, Pan, Wu, Zhang and Liang2022). These characteristics can induce various geotechnical issues, including slope instability, pavement cracking and the failure of waterproof barriers (Jiao et al., Reference Jiao, Liu, Ji, Zhang, He and Liao2022).
Research has confirmed that colloidal oxides are the primary factor underlying of the unique geotechnical properties observed in laterite (Zhang & Kong, Reference Zhang and Kong2014; Tchakouté et al., Reference Tchakouté, Bewa, Fotio, Dieuhou, Kamseu and Rüscher2021). When these oxides are removed, the physical, chemical and mechanical characteristics of laterite become comparable to those of ordinary cohesive soils. Although extensive studies have explored the macroscopic mechanical behaviour of laterite (Lu et al., Reference Lu, Dong, Wang, Ding, Wang, Hui and Wang2018; Ehrlich et al., Reference Ehrlich, Almeida and Curcio2019), the microscale interactions between colloidal oxides and clay minerals (particularly their role in controlling the microstructure and mechanical properties of laterite aggregates) remain poorly understood. A key challenge lies in deciphering the aggregate structure formed by colloidal oxides and clay particles, as this knowledge is critical for controlling laterite’s distinctive physicomechanical behaviour. However, the chemical complexity of clay minerals complicates the characterization and tracking of their aggregation processes, leaving significant gaps in our understanding of microscale dynamics. Recent advances in nanoscale geochemical characterization and simulation techniques, such as molecular dynamics (MD; Meng et al., Reference Meng, Yin, Li, Zhong, Sheng and Nie2019; Yang et al., Reference Yang, Han, Hu and He2019a; Luo et al., Reference Luo, Ou, Chen, Zhang, Xia, Zhu and Zhou2022; Liu et al., Reference Liu, Chen, Li and Zhao2023; Xu et al., Reference Xu, Liu, Deng, Wang, Xue and Wang2023), now provide powerful tools to investigate these interactions at the microscale and nanoscale.
The accuracy of MD simulations fundamentally depends on the reliability and comprehensiveness of the employed force fields. Although successive force fields developed specifically for clay minerals, including the Matsuoka–Clementi–Yoshimine (MCY), Consistent Force Field (CFF) and Consistent Valence Force Field (CVFF), have advanced research in this area, these models typically approximate clay mineral lattices as rigid slabs and do not adequately account for isomorphic substitutions within the crystal structure (Dauber‐Osguthorpe et al., Reference Dauber‐Osguthorpe, Roberts, Osguthorpe, Wolff, Genest and Hagler1988; Yotsuji et al., Reference Yotsuji, Tachi, Sakuma and Kawamura2021). The Clay Force Field (ClayFF), which was parameterized specifically for clay mineral systems (Cygan et al., Reference Cygan, Liang and Kalinichev2004), has been successfully applied to investigate the adsorption of metal ions and small molecules in clay interlayers (Li et al., Reference Li, Liu, Tang and Zhang2020), lattice expansion induced by molecular intercalation (Zhang et al., Reference Zhang, Liu, Gao, Li, Liu and Li2016) and the coupled effects of charge distribution, interlayer cations and temperature on clay-swelling behaviour (Teich-McGoldrick et al., 2015).
MD simulation has been widely employed in clay research due to its ability to reveal atomic-scale interaction mechanisms. This approach has not only successfully characterized the high crystalline water content and restricted diffusion properties within clay interlayers, but also elucidated the bond rupture mechanism underlying the yielding and failure of clay crystals (Yang et al., Reference Yang, He, Lu and Gong2019b; Zhu et al., Reference Zhu, Shen, Shao and He2021). However, most existing simulations are based on idealized, infinitely extended clay models, which fail to accurately capture the structural complexity and environmental responsiveness of natural clay aggregates (Xu & Zheng, Reference Xu and Zheng2018; Cygan et al., Reference Cygan, Greathouse and Kalinichev2021).
To address these limitations, recent studies have systematically investigated clay interfacial interactions by incorporating factors such as nanoparticle morphology and chemical composition. For instance, Cui et al. (Reference Cui, Zhaojie, Yang and Wu2015) observed the formation of solid-like ordered structures at solid–fluid interfaces, while advanced methodologies combining the weighted histogram analysis method (WHAM) with free-energy perturbation have enabled quantitative assessment of intermolecular interactions at the atomic scale (Luzhkov, Reference Luzhkov2008; Zhu et al., Reference Zhu, Whittle and Pellenq2022). These techniques provide a crucial foundation for investigating isomorphic substitution, a defining structural characteristic of clay minerals.
Isomorphic substitution represents a fundamental characteristic of clay minerals, occurring when silicon in the silicate tetrahedra (T) and/or aluminium in the octahedra (O) of phyllosilicate structures is replaced by alternative metal cations. In the tetrahedral layer, Si4+ is frequently substituted by Al3+ or Fe3+, whereas octahedral Al3+ can be replaced by various cations, including Mg2+, Fe3+, Fe2+, Cr3+, Ni2+ and Zn2+ (Bergaya et al., Reference Bergaya, Theng and Lagaly2006). This substitution mechanism typically generates a net negative charge in silicate clay minerals, creating intense electric fields ranging from 108 to 109 V m–1 near mineral surfaces (Liu et al., Reference Liu, Li, Li, Xie, Ni and Wu2014). Extensive research has examined ion-adsorption phenomena on clay mineral surfaces, particularly focusing on ionic specificity effects arising from isomorphic substitution. These studies have established distinct adsorption preference sequences: K+ > Na+ > Ca2+ > Cs+ > Ba2+ for cations (Underwood et al., Reference Underwood, Erastova and Greenwell2016) and OH− > F− > Cl− > I− for anions (Li, Reference Li2018). Molecular simulation approaches have provided valuable insights into these systems, enabling the investigation of the elastic and fracture properties of model clay minerals (Hantal et al., Reference Hantal, Brochard, Laubie, Ebrahimi, Pellenq, Ulm and Coasne2014) and the determination of nanoscale elastic constant matrices for clay–water systems in interlayer spaces (Ebrahimi et al., Reference Ebrahimi, Pellenq and Whitttle2012).
The aforementioned simulations have proven instrumental in deciphering microscale particle interaction mechanisms, providing critical insights into the interplay between clay minerals and colloidal oxide particles in natural clayey soils. To tackle persistent engineering challenges (including slope instability, differential foundation settlement, pavement cracking and waterproof barrier failure resulting from laterite’s distinctive swell–shrink behaviour, fissuring tendency and moisture sensitivity), this investigation employs MD simulations that integrate umbrella sampling and the WHAM. This approach specifically examines how isomorphic substitution influences illite–colloidal alumina interactions in laterite systems. Through a systematic analysis of potential of mean force (PMF) and differential charge density distributions across various substitution scenarios, we elucidate the fundamental interaction mechanisms between illite and colloidal alumina. The findings establish a scientific foundation for modulating laterite’s engineering characteristics and support the development of novel laterite-based materials (Nandi et al., Reference Nandi, Zaccaron, Raupp-Pereira, Arcaro, Bernardin and Montedo2023; Zhu et al., Reference Zhu, Zhang and Yong2023), offering both substantial theoretical significance and promising practical applications.
Simulation method
Model construction
Accurate characterization of mineral constituents, including their content and crystal structure, is essential for employing MD simulations to investigate the unique physicomechanical properties of natural laterite. In preliminary characterization studies (Jia et al., Reference Jia, Wu, Ren and Lin2022), X-ray diffraction (XRD) analysis revealed illite as the dominant clay mineral in laterite, comprising 22.0 wt.% of the total mass, with colloidal alumina representing 19.8 wt.%. The XRD results further identified quartz and feldspar as weathering residues, which can be considered inert components relative to the reactive clay minerals and colloidal oxides. Extensive research has established that colloidal oxides are primarily responsible for the distinctive geotechnical behaviour of laterite (Zhang & Kong, Reference Zhang and Kong2014; Tchakouté et al., Reference Tchakouté, Bewa, Fotio, Dieuhou, Kamseu and Rüscher2021). To elucidate the fundamental interaction mechanisms between these colloidal oxides and clay minerals, we developed a MD model incorporating both illite and colloidal alumina, with particular attention given to the effects of isomorphic substitution within illite.
The fundamental structural units of clay minerals consist of silicon–oxygen tetrahedra (T) and aluminium–oxygen octahedra (O). Illite exhibits a characteristic TOT layered structure, comprising two silica tetrahedral sheets sandwiching one alumina octahedral sheet. Based on experimental data from Drits et al. (Reference Drits, Zviagina, McCarty and Salyn2010), we adopted the chemical formula KAl2(SiAl)4O10(OH)2 for illite in this study. The crystal structure was carefully constructed through the following steps: first, a rectangular illite model was generated using Materials Studio software, which was subsequently trimmed along two [010] and four [110] crystallographic directions (Churakov, Reference Churakov2006) to achieve a near-hexagonal morphology (Fig. 1). The atomic configuration at the model edges was processed following a method consistent with that of Jia et al. (Reference Jia, Wu, Lin and Jiang2023): dangling oxygen bonds in the tetrahedral layer were saturated with hydrogen atoms to form silanol (Si–OH) groups, while aluminium sites at the octahedral layer edges were left coordinatively unsaturated to mimic natural clay-edge characteristics. This selective saturation strategy effectively preserves both the charge balance and reactivity of the edge structure. Natural illite commonly undergoes isomorphic substitution following Löwenstein’s rule (Loewenstein, Reference Loewenstein1954), where: (1) tetrahedral Si4+ can be replaced by Al3+; and (2) octahedral Al3+ can be substituted by Mg2+ or Fe2+. When one out of every four Si4+ in the tetrahedral sheet is replaced by Al3+, the tetrahedral sheet reaches its maximum substitution capacity (Fig. 1b). Similarly, when the octahedral sheet is substituted by Mg2+ or Fe2+ at intervals of one Al3+, the octahedral sheet achieves its maximum substitution level (Fig. 1c). Based on these maximum substitution levels, we constructed six distinct illite models representing different substitution scenarios. These modified illite flakes were then combined with colloidal alumina particles in a simulation box (Fig. 1d), yielding final crystal parameters of a = 39.13 Å, b = 43.34 Å and c = 44.84 Å, with α = β = γ = 90°. The negative charges induced by isomorphic substitution were neutralized by K+ counterions adsorbed at the surface.
Model construction: (a) front view of the illite along the {010} plane, with K+ located near the {00-1} surface; (b) and (c) represent isomorphic substitution in the silicate tetrahedral and aluminous octahedral layers, respectively; (d) is the combined model of illite with colloidal alumina particles.
Figure 2 presents the six isomorphic substitution combinations of the tetrahedral silicate (T) and octahedral aluminate (O) layers within illite clay flakes. (1) The tetrahedral and octahedral sheets are fully saturated with Si and Al atoms, respectively, without any occurrence of isomorphic substitution (No-IS). (2) In the tetrahedral sheet, Si atoms remain unsubstituted, while every alternate Al3+ in the octahedral sheet is replaced by Mg2+ (O-Mg). (3) Similarly, Si atoms in the tetrahedral sheet are unsubstituted, and every alternate Al3+ in the octahedral sheet is replaced by Fe2+ (O-Fe). (4) One Si4+ out of every four in the tetrahedral sheet is substituted by Al3+, with no substitution occurring in the octahedral Al atoms (T-Al). (5) One Si4+ out of every four in the tetrahedral sheet is substituted by Al3+, and every alternate Al3+ in the octahedral sheet is replaced by Mg2+ (T-Al&O-Mg). (6) One Si4+ out of every four in the tetrahedral sheet is substituted by Al3+, and every alternate Al3+ in the octahedral sheet is replaced by Fe2+ (T-Al&O-Fe). Based on these six isomorphic substitution scenarios for the illite clay flakes, we investigated the interactions between the {00-1} and {001} surfaces with colloidal alumina.
Six isomorphic substitution scenarios in the tetrahedral (T) and octahedral (O) layers of illite.
Simulation details
MD simulations were conducted using the LAMMPS software package (Plimpton, Reference Plimpton1995). The ClayFF was employed to describe the interactions between atoms of illite and colloidal alumina, with computational parameters provided by Cygan et al. (Reference Cygan, Liang and Kalinichev2004). To simulate the properties of an infinite system and avoid surface effects, this study employed periodic boundary conditions in all directions. Initially, the system was simulated for 1 ns under the isothermal–isobaric (NPT) ensemble at 300 K and 1 bar to achieve steady state at the specified temperature and pressure, thus obtaining the stable configurations formed by the interaction between illite and colloidal alumina particles. The interaction sites of colloidal alumina on the {00-1} and {001} surfaces of illite clay flakes were analysed. During the simulation, the particle–particle particle–mesh (PPPM) solver was used to handle Coulombic interactions with a cutoff radius of 12.5 Å and a computational precision of 1 × 10–4. The Nose–Hoover thermostat and barostat were used to control temperature and pressure, with a timestep of 1 fs. Finally, the radial distribution functions (RDFs) were computed using trajectories collected after the system reached equilibrium.
To quantitatively elucidate the interfacial interactions between illite and colloidal alumina at the microscopic level from a thermodynamic perspective, this study employed an integrated computational approach combining MD simulations, umbrella sampling and the WHAM. This strategy was adopted based on the following rationale: while conventional MD simulations can provide atomic-scale dynamic processes and equilibrium configurations, they often fail to adequately sample the transition across the free-energy barrier separating the dissociated and bound states. Therefore, umbrella sampling was implemented using the collective variables module in LAMMPS, wherein harmonic potentials with force constants of 300 and 600 kcal (mol Å–2)–1 were applied along the reaction coordinate to ensure sufficient overlap between adjacent sampling windows. Each simulation window was run for 2 ns to guarantee result reliability. It should be noted that the distance between the centres of mass of illite and colloidal alumina was selected as the reaction coordinate in this study. The initial separation distance of 16 Å between illite and colloidal alumina was systematically reduced in two stages: first, at a rate of 0.05 Å ns–1 from 16 to 10 Å, followed by a slower approach at 0.025 Å ns–1 for the remaining distance. Throughout this procedure, the illite clay sheets were held immobile. Finally, the WHAM (Kumar et al., Reference Kumar, Rosenberg, Bouzida, Swendsen and Kollman1992) was implemented to process n output files, evaluate distance sampling distributions and reconstruct the complete PMF profile. Additionally, differential charge density distributions were computed using Materials Studio software to systematically examine illite–colloidal alumina interactions across different isomorphic substitution scenarios, providing fundamental insights into their interfacial bonding mechanisms.
Results
Interaction morphologies between illite and colloidal alumina
Figure 3 presents the equilibrium interaction configuration between illite clay sheets and alumina colloids after relaxation, revealing pronounced adsorption capabilities on both the {00-1} and {001} crystallographic surfaces of illite for colloidal alumina particles. It should be noted that the negative surface charge of illite is neutralized by K+ ions, which form strong inner-sphere surface complexes preferentially above the siloxane cavities within the clay matrix (Morisada et al., Reference Morisada, Muranishi, Shinto and Higashitani2005). Equilibrium K+ concentrations vary significantly on {00-1} surfaces under different isomorphic substitution conditions. When Al atoms in the octahedral layers are substituted by Mg2+ or Fe2+ (O-Mg and O-Fe), the impact on the surface negative charge of illite is relatively minor. Conversely, substitution of Si atoms in the tetrahedral layers by Al3+ (T-Al) results in a marked increase in the surface negative charge of illite. This phenomenon may arise because the octahedral layers, positioned in the centre of the clay flake, allow for an appropriate redistribution of surrounding atomic charges to balance the charge disparities after substitution. In contrast, isomorphic substitution occurring within the tetrahedral layers, which are located at the surface of illite, directly affects the surface negative charge.
Interaction morphologies of illite {00-1} and {001} surfaces with colloidal alumina under various isomorphic substitution conditions.
Unlike traditional studies that treat clay minerals as homogeneously charged entities through simplified approaches (Cai et al., Reference Cai, Wu, Qin and Yan2025), this work quantitatively evaluates the differential responses of distinct crystallographic planes ({001} vs {00-1}) to charge variations by constructing six substitution models, thereby establishing a structure–property correlation between crystallochemical characteristics and interfacial energy. The presence of K+ ions on the {00-1} surface of illite creates distinctly different interaction patterns with colloidal alumina compared to the {001} surface. This differential behaviour stems from electrostatic complementarity: while K+ carries a positive charge, the oxygen and aluminium atoms in colloidal alumina possess negative and positive charges, respectively. The heterogeneous charge distribution across colloidal alumina surfaces drives negatively charged oxygen atoms to preferentially orientate towards positively charged K+ sites. As a result, colloidal alumina particles demonstrate enhanced adsorption stability in K+-enriched regions of the {00-1} surface, with their constituent oxygen atoms undergoing conformational adjustments to minimize the K+–O distances.
Figure 4 presents the RDFs of K+ with O and Al atoms in colloidal alumina on the illite {00-1} surface. It can be seen from the RDF curves that the distance between the K atom and the O atom is closer than that between the K atom and the Al atom. Isomorphic substitution substantially influences these interactions, as evidenced by peak position variations: for the non-substituted case (No-IS), the K–O peak occurs at 2.69 Å, while substitutions (O-Mg, T-Al&O-Fe and T-Al&O-Mg) shift this peak inward to 2.56 Å, indicating enhanced K+–O attraction. Notably, T-Al and O-Fe substitutions result in weaker interactions (K–O peaks at 3.06 and 2.94 Å, respectively) compared to No-IS. Regarding K–Al interactions, No-IS shows its first peak at 3.69 Å, while substitutions generally strengthen this interaction – most remarkably for T-Al&O-Fe (peak position at 3.06 Å), followed by T-Al&O-Mg. Interestingly, among octahedral substitutions alone, O-Mg induces significantly stronger K+–colloid interactions than O-Fe.
RDFs of K+ with O and Al atoms of colloidal alumina on the {00-1} surface of illite. g(r) = relative probability density of finding another particle at a distance r from a reference particle.
Further analysis of isomorphic substitution effects was conducted through RDF analyses between substitution atoms and colloidal alumina components. Figure 5a,b presents the RDFs for octahedral Fe(O-Fe) and Mg(O-Mg) substitutions, respectively, revealing distinct interaction patterns: O-Fe shows strong Fe–O bonding (shorter distances) with weaker Fe–Al interactions, while O-Mg demonstrates the opposite trend, with preferential Mg–Al bonding. This phenomenon may be attributed to the distinct electronic and structural characteristics of the substituting cations. While Fe atoms can exhibit either divalent (Fe2+) or trivalent (Fe3+) states depending on their electron configuration, Mg atoms predominantly maintain a stable divalent (Mg2+) state. The combined substitution cases (T-Al&O-Fe and T-Al&O-Mg) are shown in Fig. 5c,d, with the first-shell coordination distances of tetrahedral Al substituents with respect to Al/O atoms in colloidal alumina being 2.94, 3.56 Å and 3.56, 4.06 Å, respectively. These results confirm that T-Al&O-Fe induces the strongest colloidal alumina adsorption, followed by T-Al&O-Mg, which is consistent with the previously observed K+-mediated adsorption trends. The systematic RDF analysis provides atomic-scale evidence for the cation-specific effects of isomorphic substitution on colloidal alumina interactions.
In contrast to the K+-rich {00-1} surface, the {001} surface of illite lacks K+ ions, enabling direct adsorption of colloidal alumina onto the clay matrix. As shown in Fig. 3, in cases of O-Mg and T-Al, an oxygen atom from the colloidal alumina is preferentially drawn to the sites of substitution (Fig. 3a,b for the {001} surface). For T-Al&O-Fe and T-Al&O-Mg, the colloidal alumina is drawn towards the areas immediately above the most densely substituted regions (Fig. 3c,d for the {001} surface), where the proximity between the Al and O atoms of the colloidal alumina and the substituting atoms is notably reduced compared to other cases. In instances of O-Fe and No-IS, the colloidal alumina gravitates towards the cavities above the tetrahedral silicate structures (Fig. 3e, f for the {001} surface).
Figure 6 presents the RDFs between substitution atoms on the {001} surface of illite and atoms in colloidal alumina. When octahedral Al atoms are replaced by Fe or Mg (i.e. O-Fe and O-Mg), the illite exhibits distinct interaction patterns with colloidal alumina: for the O-Fe system, illite preferentially interacts with oxygen atoms of colloidal alumina, as evidenced by the first Fe–O RDF peak at 6.06 Å (vs Fe-Al at 6.31 Å). For the O-Mg system, the interaction shifts to aluminium atoms, with the first Mg–Al RDF peak at 5.19 Å (vs Mg–O at 5.56 Å), indicating stronger adsorption resulting from Mg substitution. In T-Al&O-Fe and T-Al&O-Mg systems (Fig. 6c,d), the Al atoms of colloidal alumina preferentially interact with illite, resulting in significantly more stable adsorption compared to octahedral substitution alone.
RDFs of substitution atoms on the illite {00-1} surface and atoms from colloidal alumina. g(r) = relative probability density of finding another particle at a distance r from a reference particle.
RDFs of substitution atoms on the illite {001} surface and atoms from colloidal alumina. g(r) = relative probability density of finding another particle at a distance r from a reference particle.
The adsorption behaviour of colloidal alumina on illite’s {00-1} and {001} surfaces under various substitution conditions reveals that clay mineral–colloidal oxide interactions are governed by three key factors: (1) the nature of charge-balancing cations; (2) the chemical identity of substituting atoms; and (3) their crystallographic substitution sites. These parameters collectively modulate the clay surface charge distribution, thereby regulating the colloidal oxide adsorption efficiency.
PMFs between illite and colloidal alumina
The PMF calculations characterize the interaction energetics between colloidal alumina and illite clay sheets during their approach (Fig. 7). Analysis of the PMF profiles reveals three distinct interaction regimes: (1) weak repulsion at large separation distances; (2) progressively stronger attraction with decreasing interparticle distance that reaches maximum binding energy; followed by (3) rapid development of repulsive forces, forming a characteristic energy well. Notably, the {00-1} surface exhibits significantly stronger attraction to colloidal alumina than the {001} surface, which is attributed to enhanced interactions mediated by charge-balancing K+ cations that form electrostatic pairs with negatively charged oxygen atoms (O2–) in colloidal alumina. Isomorphic substitution systematically modulates these interactions, with the {00-1} surface showing reduced attraction post-substitution, whereas the {001} surface displays increased binding affinity. Specifically, octahedral Mg substitution induces greater PMF changes than Fe substitution on the {00-1} surface (Fig. 7b,c), indicating weaker attraction. This trend persists when combined with tetrahedral Al substitution (Fig. 7e,f), where Mg substitution again results in weaker attraction compared to Fe substitution. These results conclusively demonstrate that octahedral Mg substitution more significantly weakens {00-1} surface interactions with colloidal alumina than Fe substitution, providing atomic-scale insights into cation-specific effects on clay–colloid interactions.
PMFs of interaction between illite and alumina: (a) no isomorphic substitution occurred (No-IS); (b) and (c) represent Al atoms in the octahedral sheet being replaced by Mg (O-Mg) and Fe (O-Fe), respectively; (d) Si atoms in the tetrahedral sheet are substituted by Al atoms (T-Al); (e) and (f) represent replacement of Si atoms in the tetrahedral sites by Al atoms, and Al atoms in the octahedral sheet are replaced by Mg (T-Al&O-Mg) and Fe (T-Al&O-Fe) atoms, respectively.
Figure 8 provides a schematic representation of characteristic PMF features, where points a and c correspond to maximum repulsion and maximum attraction, respectively, whereas points b and d indicate positions with zero interaction force between particles. To clarify the effects of isomorphic substitution on illite–alumina interactions in laterite, Table 1 summarizes numerical values at key positions on the PMF curves. For non-substituted systems, the {00-1} surface of illite exhibits maximum repulsion of 8.45 kcal mol–1 at 8.65 Å separation from alumina particles and maximum attraction of –51.16 kcal mol–1 at 1.75 Å. Octahedral Mg substitution yields the weakest attraction among all substitution cases, showing maximum attraction of –25.99 kcal mol–1 at 4.11 Å, followed by combined tetrahedral Al and octahedral Mg substitution (T-Al&O-Mg). Notably, combined tetrahedral Al and octahedral Fe substitution (T-Al&O-Fe) produces maximum repulsion of 13.48 kcal mol–1 at 6.06 Å. Generally, isomorphic substitution enhances attraction between the {001} surface and colloidal alumina to varying degrees depending on substitution type and location (Table 1). The non-substituted {001} surface shows maximum attraction of only –17.22 kcal mol–1, whereas T-Al&O-Fe substitution increases this value to –29.27 kcal mol–1.
Schematic diagram of typical positions in a PMF curve: points a and c denote the locations of maximum repulsive and attractive forces, respectively, and points b and d represent positions where the interaction force is zero.
Typical positional characteristics of PMF curves under various substitution conditions.
The simulation results demonstrate that the presence of K+ ions induces significantly more complex interactions between the {00-1} surface of illite clay sheets and colloidal alumina compared to the {001} surface. As shown in Fig. 9, the PMF curves for O-Mg and T-Al substituted systems exhibit two local energy minima (points A and C) separated by a transition state (point B). In the O-Mg and T-Al substitution simulations, the attraction between illite and colloidal alumina reaches local maxima at 1.25 and 1.24 ns, respectively (point A in Fig. 9), decreases to the transition state (point B) at 1.55 and 1.50 ns,, respectively, and subsequently increases to the second local minimum (point C). The corresponding interaction configurations at points A, B and C clearly illustrate the dynamic evolution of illite–alumina interactions. During the simulation, positively charged K+ ions preferentially form electrostatic bonds with negatively charged O2– ions in alumina, resulting in a relatively stable coordination between one terminal oxygen atom and K+ at point A. Simultaneously, interactions between other oxygen atoms and K+ create unbalanced electrostatic forces that induce rotational displacement of colloidal alumina, forming the transition state at point B. Furthermore, the negatively charged sites generated by isomorphic substitution in illite strongly attract positively charged Al atoms, driving further reorientation of colloidal alumina to minimize Al–substituent distances and establish the more stable configuration at point C. This multistage interaction process, governed by the combined effects of K+ cations and substitution-induced negative charges, explains the characteristic PMF profile, with its distinct local minima and transition states. The appearance of these metastable states (points A and B) in the PMF curves fundamentally accounts for the high sensitivity and strong structural stability observed in natural laterite clays.
PMFs of interaction between the {00-1} surface of illite and colloidal alumina during (a) O-Mg and (b) T-Al substitution: points A and C, respectively, represent the locally stable states, and point B denotes the transition state. The blue dashed circles and red curved arrows are used to highlight the position and movement direction of the colloidal alumina.
Differential charge density distributions
To systematically investigate the effects of isomorphic substitution on illite–alumina interactions, we calculated differential charge density distributions using Materials Studio software under various substitution conditions, obtaining cross-sectional charge density profiles (Fig. 10). The analysis reveals highly non-uniform charge distributions on both illite clay sheets and colloidal alumina particles. In non-substituted systems, Al atoms in octahedral sites and Si atoms in tetrahedral sites predominantly lose electrons while oxygen atoms gain electrons (Fig. 10b,c). Similarly, in colloidal alumina, Al atoms exhibit electron depletion and O atoms show electron accumulation. Octahedral Mg substitution results in reduced electron loss at substitution sites compared to non-substituted cases (Fig. 10e), indicating a slight decrease in overall charge density that explains the weakest attraction observed in Mg-substituted systems. In contrast, octahedral Fe substitution leads to electron accumulation around Fe centres (Fig. 10g), creating enhanced charge heterogeneity that facilitates exceptionally strong interfacial bonding, evidenced by the minimal equilibrium separation distance of 1.13638 Å at maximum attraction. Tetrahedral Al substitution shows moderately reduced electron loss at substitution sites, although this effect is less pronounced than in octahedral Mg substitution cases.
Differential charge density maps of illite clay flakes and colloidal alumina, illustrating: (a) a schematic representation of the differential charge density sections; and the A, B and C sections for (b) the non-substituted case, (e) octahedral Mg substitution, (f) tetrahedral Al substitution and (g) octahedral Fe substitution.
The variation in charge density of illite clay sheets induced by isomorphic substitution plays a crucial role in governing colloidal alumina adsorption. Specifically, Mg and Al substitutions decrease the charge density of illite, whereas Fe substitution increases it. The presence of charge-balancing K+ cations on the {00-1} surface introduces additional complexity to illite–alumina interactions, as positively charged K+ forms electrostatic attractions with oxygen atoms in colloidal alumina, whereas the negatively charged substitution sites preferentially interact with the aluminium atoms of alumina. In contrast, the {001} surface, lacking K+ cations, exhibits more direct interactions, where the negatively charged substitution sites immediately engage with colloidal alumina. Consequently, when isomorphic substitutions occur in both tetrahedral and octahedral sites, the interaction strength between illite and colloidal alumina increases proportionally with the increased charge density.
Discussion
Although the MD simulations presented in this study yield significant insights into the role of isomorphic substitution at the atomic scale, it is crucial to acknowledge the inherent limitations of the theoretical framework employed.
The primary limitation stems from the use of classical, non-reactive force fields (e.g. ClayFF). These force fields rely on fixed partial charges and predefined interaction potentials, which do not permit chemical reactions, bond breaking/forming or dynamic changes in charge states to occur during the simulation. Consequently, the process of isomorphic substitution itself is not directly observed; rather, its thermodynamic feasibility and structural consequences are inferred from the computed PMFs and equilibrium configurations. For example, although iron can exhibit valence transitions between Fe2+ and Fe3+ in natural environments, it was modelled with a fixed Fe2+ valence state in this study, reflecting the dominant oxidation state under the simulated conditions, as the classical force field cannot simulate electron transfer processes.
Furthermore, the simulation models constitute a necessary simplification of the highly heterogeneous lateritic soil system. This study focuses on a well-defined binary system (illite and colloidal alumina) in an aqueous environment, isolating the phenomenon of interest but inevitably omitting the influence of other natural soil constituents, such as organic matter, diverse mineral phases and complex porewater chemistry. The use of simplified model structures, while providing molecular-level clarity, may not fully capture the complexity of natural clay edges and oxide surfaces.
Despite these limitations, the power of the MD approach lies in its unique ability to test specific mechanistic hypotheses and provide thermodynamic data that are extremely challenging to obtain experimentally. The PMF profiles and structural analysis offer a quantitative foundation for interpreting macroscopic engineering properties. The findings of this work – particularly the identification of favourable interaction sites and the quantitative impacts of different substitution types on binding affinity – generate specific, testable predictions. These predictions can guide future experimental research – for instance, by targeting the predicted local atomic environments around substitution sites using advanced spectroscopic techniques such as solid-state nuclear magnetic resonance or X-ray absorption spectroscopy for validation. Looking forward, the insights gained here could also motivate the application of more advanced simulation methods, such as reactive force fields or quantum mechanical/molecular mechanical approaches, which could potentially model the dynamic process of isomorphic substitution explicitly. Although the MD simulations presented in this study yield significant insights into the role of isomorphic substitution at the atomic scale, it is crucial to acknowledge the inherent limitations of the theoretical framework employed.
Conclusions
Laterite’s unique engineering properties (including its high strength, low compressibility, low permeability and strong ion-exchange capacity) have enabled its widespread application across the construction, environmental remediation, energy and chemical industries. Although colloidal oxides have been identified as the fundamental cause of these distinctive geotechnical characteristics, the microscopic mechanisms governing their interactions with clay minerals remain incompletely understood. This knowledge gap particularly concerns how isomorphic substitution and interlayer cations collectively influence colloidal oxides’ control over laterite’s microstructural organization and mechanical behaviour. Building upon the XRD characterization of natural laterite, this study systematically investigates isomorphic substitution effects on illite–alumina interactions through MD simulations incorporating umbrella sampling and the WHAM. Our approach elucidates stable configurations, PMF profiles and charge density distributions, ultimately revealing the atomic-scale interaction mechanisms under various substitution conditions. Key findings include the following.
Interfacial stabilization mechanism
The heterogeneous charge distribution across illite–alumina interfaces creates a complex stabilization mechanism. On K+-rich {00-1} surfaces, electrostatic attractions arise between (1) K+ cations and oxygen anions (O2–) of alumina and (2) anionic sites induced by isomorphic substitution and the aluminium cations (Al3+) of alumina. These competing interactions generate two metastable states, separated by a rotational transition state. The distinct PMF signatures of these states explain the unique sensitivity and structural stability of laterite.
Critical role of substitution type
The identity and location of substituent atoms critically determine the interaction strength: (1) octahedral Mg and tetrahedral Al substitutions reduce electron loss at the substitution sites, decreasing charge density and consequently weakening the attraction; and (2) octahedral Fe substitution, conversely, promotes electron accumulation around the Fe centres, enhancing charge heterogeneity and resulting in exceptionally tight binding (equilibrium distance of 1.14 Å at maximum attraction).
Surface-specific responses
The sensitivity of the PMF to charge redistribution exhibits distinct, surface-dependent trends: (1) on {00-1} surfaces, substitutions attenuate attraction, an effect most pronounced for Mg; (2) on {001} surfaces, binding strength increases proportionally with charge density enhancement, maximized by Fe substitution; and (3) this dichotomy stems from substituent valence states: Mg2+/Al3+ maintain fixed oxidation states, whereas the behaviour of iron, which can exist as both Fe2+ and Fe3+ in natural systems, was modelled as Fe2+ in this study, consistent with the dominant valence state under the simulated conditions.
CRediT authorship contribution
Jing Yan: Methodology, Software, Writing – original draft. Daoyong Wu: Supervision, Conceptualization, Investigation, Writing – review & editing. Zipeng Qin: Data curation, Investigation. Jie Wang: Formal analysis, Software. Weizhong He: Resources, Formal analysis. Qingliang Zhu: Resources, Investigation. Yan Tian: Software, Validation.
Data availability
No data were used for the research described in the article.
Acknowledgements
None.
Financial support
This research was supported by the National Natural Science Foundation of China (Grant No. 42002280), the Science and Technology Foundation of Guizhou Province (No. ZK[2022]Key018) and Zhejiang Provincial Natural Science Foundation of China under Grant No. LTGG24E090002.
Competing interests
The authors declare none.
