# Atomistic Computer Simulations of Uranyl Adsorption on Hydrated Illite and Smectite Surfaces

Anna D Krot, Irina E Vlasova, Evgeny V Tararushkin, Andrey G. Kalinichev

► To cite this version:

Anna D Krot, Irina E Vlasova, Evgeny V Tararushkin, Andrey G. Kalinichev. Atomistic Computer Simulations of Uranyl Adsorption on Hydrated Illite and Smectite Surfaces. Minerals, 2024, 14 (1), pp.109. 10.3390/min14010109. [hal-04412668](https://hal.science/hal-04412668v1)

HAL Id: hal-04412668

https://hal.science/hal-04412668v1

Submitted on 19 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Atomistic Computer Simulations of Uranyl Adsorption on Hydrated Illite and Smectite Surfaces

Anna D. Krot<fnref n="1" />, Irina E. Vlasova<fnref n="1" />, Evgeny V. Tararushkin<fnref n="1" />,<fnref n="2" /> and Andrey G. Kalinichev<fnref n="3" />,*

**Abstract:** A quantitative understanding of the molecular-scale mechanisms of radionuclide sorption on different clay minerals is crucial for the development and safe implementation of geological nuclear waste disposal technologies. We apply classical molecular dynamics (MD) computer simulations to study the adsorption of uranyl on the external basal surfaces of two typical clay models. In the illite model, negative charge is primarily localized in the tetrahedral sheets, while in the lower-charge smectite model, the isomorphic substitutions are introduced in the octahedral sheet. The comparison of atomic density distributions at the clay surfaces and adsorption-free energies profiles as a function of distance from these surfaces demonstrates that overall U behavior at the basal clay surface is quite similar for illite and smectite. Uranyl is sorbed as a mixture of outer-sphere aqua complexes [UO₂(H₂O)₅]<fnref n="2" />⁺ and hydrolyzed aqua complexes [UO₂(H₂O)₄₅OH]⁺ on both surfaces. The structural and compositional differences between the models do not greatly affect the uranyl's nearest coordination environment and are mainly reflected in the specific localization and orientation of the uranyl ions at both surfaces and in the magnitude of the adsorption-free energies. The observed quantitative characteristics of uranyl interactions with illite and smectite surfaces will help to better understand U behavior during the sorption process on clay minerals for the entire range of mixed-layer illite–smectite structures. A comparison of two versions of the ClayFF force field in the simulations made it possible to more accurately and quantitatively evaluate some subtle features of the uranyl–clay interactions and to obtain a more precise composition of uranyl complex with the modified ClayFF force field (ClayFF-MOH).

**Keywords:** uranyl; clay; illite; smectite; basal surface; structural properties; adsorption-free energy; molecular dynamics

## 1. Introduction

Due to their low hydraulic permeability, high sorption capacity to various radionuclides, and long-term stability, clay host rocks are now universally considered one of the most important components of the complex natural and engineered barrier systems for radioactive waste disposal and storage in deep geological formations [1–4]. Moreover, clays are an ubiquitous constituent of many geological formations [5,6]. Therefore, it is highly important to quantitatively understand and predict the behavior of radionuclides in their interaction with clay minerals under the conditions of the geological disposal of spent nuclear fuel and radioactive waste [7–11]. Among radionuclides, uranium deserves special attention, since its amount and role in the nuclear industry is the highest [7].

Various aspects of uranium sorption, in the form of uranyl cation, UO₂<fnref n="2" />⁺, on clay minerals have already been extensively studied [12–25]. It was shown that clay minerals demonstrate a high sorption capacity towards uranyl in the wide range of pH ~4–8. In

<figure><img src="image_2.png" /></figure>

### check for updates

**Citation:** Krot, A.D.; Vlasova, I.E.; Tararushkin, E.V.; Kalinichev, A.G. Atomistic Computer Simulations of Uranyl Adsorption on Hydrated Illite and Smectite Surfaces. *Minerals* **2024**, *14*, 109. https://doi.org/10.3390/min14010109

Academic Editors: Katharina Müller and Norbert Jordan

Received: 9 July 2023
Revised: 15 January 2024
Accepted: 15 January 2024
Published: 19 January 2024

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ Department of Chemistry, Lomonosov Moscow State University, 119991 Moscow, Russia; anna.d.krot@gmail.com (A.D.K.); irinaeng@gmail.com (I.E.V.); etararushkin@hse.ru (E.V.T.)</footnote>

<footnote n="2">² International Laboratory for Supercomputer Atomistic Modelling and Multi-Scale Analysis, National Research University Higher School of Economics, 101000 Moscow, Russia</footnote>

<footnote n="3">³ Laboratoire SUBATECH (UMR 6457—IMT Atlantique, Nantes Université, CNRS/IN2P3), 44300 Nantes, France</footnote>

<footnote>* Correspondence: kalinich@subatech.in2p3.fr</footnote>
this acidity range, uranyl is partially hydrolyzed and forms mixed hydroxo-carbonate complexes and adsorbs via the formation of inner-sphere complexes at clay mineral edge sites. At a higher pH, sorption decreases as uranyl forms highly soluble $[UO_2(CO_3)_3]^{4-}$ anionic carbonate complexes in solution. A significant effect of $CO_2(g)$ and competing cations, especially $Ca^{2+}$, has been demonstrated in numerous investigations [12–16]. Several studies were aimed at revealing the local nanoscale structure of adsorbed uranyl ions. For this, X-ray absorption spectroscopy (XAS) and, in particular, its part beyond the absorption edge, extended X-ray absorption fine structure spectroscopy (EXAFS), is widely used. When the uranyl ion is adsorbed at the SiOH/AlOH edge sites of a clay mineral, the EXAFS results typically demonstrate a split of the equatorial coordination sphere of O atoms at distances of ~2.3 and 2.5 Å with coordination numbers (CNs) of two and three, respectively, and Si/Al atoms at distances of ~3.1, 3.3, and 3.6 Å. At a pH < 4 and an ionic strength in the range of ~0.01–1.0 M, uranyl cation in solution is present as a non-hydrolyzed aqueous ion $[UO_2(H_2O)_5]^{2+}$ and undergoes adsorption on clay minerals by an ion-exchange mechanism, forming outer-sphere complexes at external basal surfaces and in the interlayer space. When uranyl forms outer-sphere complexes with only $H_2O$ molecules in its local coordination, the EXAFS spectra reveal a single coordination sphere of equatorial O atoms at ~2.40 Å with CN of five to six and the absence of contributions from Si/Al coordination spheres [17–25].

In addition to numerous experimental studies, atomistic computer simulations are also now widely used to obtain detailed quantitative information on the adsorption of aqueous uranyl species on various clay surfaces at a fundamental molecular scale [26–35]. Classical atomistic simulations can mostly be used to characterize the adsorption at the external basal surfaces and interlayer spaces of clay minerals. To analyze inner-sphere complexes at the edge sites, the chemical bond formation should be accurately taken into account, which requires a quantum chemical approach (e.g., density functional theory, DFT) [36–41]. However, although the interactions of $UO_2^{2+}$ can only be partially evaluated in classical atomistic simulations, they are still highly beneficial because they allow us to study much larger model systems for much longer time intervals than it is currently possible with quantum approaches [42–44].

ClayFF [42] is a classical molecular modeling force field that is currently one of the most widely used and thoroughly tested force fields for atomistic computer simulations of clay minerals, their fluid interfaces, and other related materials under a very wide range of conditions [43,44]. Its original version [42] has been recently upgraded by adding several new terms that allow for a more accurate description of the Metal-O-H (MOH) angular bending motion in hydroxide crystals and on their edges [44]. Recent molecular simulations of brucite, $Mg(OH)_2$ [45], and 10-Å phase—an important member of the family of dense hydrous magnesium silicates (DHMS)—which play a key role in transporting and storing water in the Earth's mantle at subduction zones, $Mg_3Si_4O_{10}(OH)_2 \cdot nH_2O$ [46], have already clearly demonstrated the importance of the MOH terms for a more accurate description of mineral systems containing structural hydroxyl groups. In addition to a more accurate reproduction of the crystal properties of hydroxide minerals themselves, the introduction of the MOH terms also affects the structural and dynamic properties of interfacial aqueous solutions compared to the original, less accurate version of ClayFF [42,47]. Moreover, the introduction of the M-O-H angle bending terms in the ClayFF parametrization made it possible to simulate crystal structures at higher temperatures and pressures than with the original ClayFF version [45].

Illite and smectite crystal structures are formed by so-called 2:1 TOT layers, where one sheet of octahedrally (O) coordinated $AlO_6$ units is sandwiched from both sides by two sheets of tetrahedrally (T) coordinated $SiO_4$ units. The T and O sheets are bound together via bridging oxygen atoms of the $AlO_6$ and $SiO_4$ polyhedra. Such 2:1 TOT-layered structures can easily accommodate various isomorphic substitutions ($Al^{3+}$ for $Si^{4+}$ in the T-sheets, $Mg^{2+}$ or $Fe^{2+}$ for $Al^{3+}$ in the O-sheets, etc.), leading to a negative structural charge of the TOT layers. This charge is compensated by the presence of cations in the interlayer space of clay structures.
In smectite, the substitutions are mainly localized in the octahedral sheet, with a common crystallographic unit cell formula of $M_x[Al_{4-x}(Fe,Mg)_x]Si_8O_{20}(OH)_4 \cdot nH_2O$, where M is, typically, $Ca^{2+}$, $Na^+$, and $K^+$, and the extent of substitution (x) varies from ~0.2 to ~0.6. The relatively low structural charge and high hydration energy of the interlayer cations promote the hydration, expansion, and swelling of these layered structures in the presence of water. On the contrary, illite has a non-expandable clay structure, mainly with tetrahedral substitutions, such as $K_{x+y}(Al_{4-x}Mg_x)[Si_8-yAl_y]O_{20}(OH)_4$; $x + y$ is typically in the range from ~0.6 to ~1, i.e., significantly higher than in smectites [48]. $K^+$ ions, typically representing the charge-compensating interlayer cations in illite, remain nonhydrated.

In geochemical systems, the process of illitization takes place over time, and clay minerals are often present as mixed-layer illite-smectite materials [49–51]. Their more disordered crystal structures are often not very well resolved and their description is based on certain assumptions, such as the Markovian model [52], or the fundamental particle model [53,54]. The former model is based on a statistical approach and the clay mineral structure is represented as a sequence of illitic and smectitic layers, where the nature of each layer is dependent on several previous layers. According to the latter model, the mixed-layer clay structure consists of multi-layer particles of illitic composition, and the smectitic behavior is the consequence of water adsorption at the external surfaces of these particles.

Here we use classical molecular dynamics (MD) computer simulations to characterize in detail the interaction of hydrated uranyl ions with the clay surfaces of the two end-members of the mixed-layer clay materials. The main purpose of this work is to quantify the effects of the specific structural and compositional features of the external basal surfaces of mixed-layer illite-smectite clays on the structure and adsorption capacity of hydrated uranyl species at the nanoscale. Additionally, we evaluate how recent improvements in the ClayFF force field affect the results of these simulations.

## 2. Models and Methods

### 2.1. Structural Models

Our model of “illite” was constructed from the ideal pyrophyllite crystal structure ($a = 5.2021$ Å, $b = 8.9797$ Å, $c = 10.226$ Å, $\alpha = \gamma = 90^{\circ}$, $\beta = 101.57^{\circ}$) [55], while the “smectite” model was prepared from the crystal structure of montmorillonite ($a = 5.18$ Å, $b = 8.98$ Å, $c = 15.00$ Å, $\alpha = \beta = \gamma = 90^{\circ}$) [56]. In both cases, $Mg^{2+}$ and $Al^{3+}$ substitutions were randomly introduced into the octahedral or tetrahedral sheets of the TOT layers, respectively. The Loewenstein rule was obeyed to avoid two adjacent substitutions [57]. The resulting substituted crystallographic unit cells were then multiplied along the z-axis to build crystals consisting of 7 TOT layers. In the interlayer space of the illite model, dehydrated $K^+$ ions were introduced in the amount equal to the net negative charge of the substituted TOT layer. Similarly, the interlayer space of the smectite model was filled with $[Ca(H_2O)_6]^{2+}$ to compensate for the negative charge of the surface. The illite model was specifically constructed with the structural negative charge close to the ISCz-1 mixed-layer clay mineral ($Z_T = 0.81$ per formula unit (p.f.u.)), $Z_O = 0.48$ p.f.u.) [58], which is currently under our experimental investigation. From the known amount of fixed ($K^+$, $Na^+$) and exchangeable ($Ca^{2+}$, $Mg^{2+}$) cations and illite and smectite fractions in the mixed-layer mineral (0.3 and 0.7, respectively), the charges of illitic and smectitic layers were calculated to be 1.48 p.f.u. and 0.86 p.f.u., respectively. To estimate tetrahedral and octahedral charge in illitic and smectitic layers separately, a correlation plot from [54] was used and the following charges were specified: in illitic layers, $Z_O = 0.3$ p.f.u. and $Z_T = 1.18$ p.f.u; in smectitic layers, $Z_O = 0.76$ p.f.u. and $Z_T = 0.1$ p.f.u. These charges were approximated to simplify the model construction. Therefore, the illite model contained 5 octahedral substitutions and 8 tetrahedral substitutions in a single TOT sheet composed of 10 unit cells ($2 \times 5$ in x and y directions, respectively), while the smectite model was constructed with only octahedral charge localization and the TOT sheets contained 8 octahedral substitutions.

The resulting structural formulas of modeled clay minerals are $K_{1.3}(Al_{3.5}Mg_{0.5})[Si_{7.2}Al_{0.8}]O_{20}(OH)_4$ and $[Ca(H_2O)_6]_{0.4}(Al_{3.2}Mg_{0.8})[Si_8]O_{20}(OH)_4$ for the illite and smectite models,
respectively. In the illite model, the charge distribution was calculated to be equal to that of the external surface of mixed-layer illite-smectite ISCz-1 in terms of the fundamental particle theory [53], while in the smectite model only the total octahedral charge of the ISCz-1 composition was considered, as it should be in smectitic layers according to the statistical Markovian model [50], and therefore appeared to be slightly higher than it is common for smectite minerals.

To simulate the interface of illite/smectite with a UO₂(OH)⁺ aqueous solution, a system consisting of 4 × 10 × 7 illite/smectite unit cells was cleaved parallel to the (001) crystallographic plane in the middle of the interlayer space. After cleavage, half of the interlayer K⁺/Ca²⁺ ions located at the external surfaces were randomly retained by each of the two created basal surfaces. These surfaces were brought into contact with a layer of UO₂(OH)⁺ aqueous solution approximately 80 Å thick. After the insertion of the aqueous layer, the dimensions of the simulation supercells were approximately 36 × 52 × 144 Å³ and 36 × 52 × 156 Å³ for the illite and smectite models, respectively. In both cases, the number of H₂O molecules (4856) introduced into the nanopore corresponded to the density bulk liquid water under ambient conditions, 1 g/cm³. Eight [UO₂OH]⁺ complexes were initially uniformly distributed in the aqueous nanopores, resulting in a [UO₂²⁺] concentration of ~0.1 M. The solution ionic strength was not adjusted specifically and could be assumed to be equal to 0 M. A few OH⁻ anions were added to the aqueous phase to ensure the electrostatic neutrality of the simulated systems in both cases. In this work, partly hydrolyzed [UO₂OH]⁺ and aqua non-hydrolyzed complexes are understood as uranyl cation with the first coordination sphere (as defined by the first peak of the radial distribution function; see below) of H₂O/OH⁻, although water molecules in this coordination shell are dynamic (rapidly exchanging with other molecules in solution).

## 2.2. Force Field Parameters

All MD simulations were performed with the new modified version of the ClayFF force field (ClayFF-MOH, [44]) and with the original parametrization (ClayFF-orig, [42]) for comparison. ClayFF assigns partial charges to different structural oxygen atoms by taking into account, when necessary, their local coordination to the tetrahedral and/or octahedral substitutions in the crystal structure of minerals. The values of the additional parameters for the Metal-O-H (M-O-H) angle bending terms of the more recent ClayFF-MOH version were determined for the Mg-O-H and Al-O-H angles by fitting the results of classical MD simulations to the structural and dynamic results of DFT calculations for brucite, Mg(OH)₂, and gibbsite, Al(OH)₃, respectively, using a simple harmonic functional form [44]:

$$E_{\text{MOH}} = k(\theta - \theta_0)^2 \quad (1)$$

where $θ_{0,\text{MgOH}} = 120^{\circ}$, $k_{\text{MgOH}} = 25.104 \text{ kJ} \cdot \text{mol}^{-1} \cdot \text{rad}^{-2}$, $θ_{0,\text{AlOH}} = 116^{\circ}$, and $k_{\text{AlOH}} = 62.76 \text{ kJ} \cdot \text{mol}^{-1} \cdot \text{rad}^{-2}$. These terms are used to properly constrain the librational motions of the structural hydroxyl groups in our clay models. In addition, for a more accurate modeling of the stretching vibrations of the hydroxyl ions, the Morse potential was used, instead of the original harmonic potential [59]:

$$E_{\text{OH}} = D_0 \left[ 1 - e^{-\alpha(r-r_0)} \right]^2 \quad (2)$$

where $D_0 = 552.288 \text{ kJ} \cdot \text{mol}^{-1}$, $\alpha = 2.1815 \text{ Å}^{-1}$, and $r_0 = 0.945 \text{ Å}$.

The atomic partial charges and Lennard-Jones parameters for the uranyl cation, UO₂²⁺, were taken as for a non-polarized species. Although the polarizability of uranyl is important for complexation onto clay mineral surfaces, models with non-polarized UO₂²⁺ have proven to be reliable for obtaining close to the experimental data results [27,28,32,34,60]. H₂O molecules were modeled using the flexible extended simple point charge (SPC/E) potential [61], and the intermolecular interaction parameters for aqueous OH⁻ consistent with this water model were taken from Greathouse et al. [26].
## 2.3. *Simulation Details*

MD simulations were performed using the LAMMPS software package [62,63]. Periodic boundary conditions [64,65] in three dimensions were applied in all simulations. A cutoff radius of 12.5 Å was used for the short-range intermolecular forces, and the particle–particle particle–mesh method (PPPM) was used to account for the long-range electrostatic interactions [65–67]. Standard Lorentz–Berthelot mixing rules ($ε_{ij} = \sqrt{ε_{ii}ε_{jj}}$; $σ_{ij} = (σ_{ii} + σ_{jj})/2$) were applied to calculate the Lennard-Jones parameters of interaction between unlike atoms. The Newtonian equations of motion were numerically integrated with a time step of 1.0 fs via velocity Verlet algorithm [64,65,68]. The constructed models were initially equilibrated for 5.0 ns at a temperature of 298 K and pressure of 1 bar using the Nosé–Hoover thermostat [69,70]. After the NPT equilibration, the equilibrium NVT-simulation was performed for 20.0 ns using the Nosé–Hoover thermostat [69,70] to collect the dynamic trajectories of all atoms of the simulated system for further statistical analysis. The quantitative analysis of the structure of clay–solution interfaces was performed following the well-developed approaches that are commonly used for such kinds of interfacial simulations [47].

In addition, a series of constrained MD simulations were performed to obtain more information on the uranyl coordination at the basal surfaces of illite and smectite models and to quantitatively assess its adsorption-free energy as a function of the distance from the surface. The metadynamics technique implemented in the PLUMED code [71] was used for this part of the simulation project. In these simulations, the distance z between the center of mass of the uranyl ion and the center of mass of the nearest ditrigonal cavity without tetrahedral substitutions on the clay basal surfaces was constrained by a collective variable (CV). The CV was limited to a distance of 12 Å from the clay surface. To limit the movements of uranyl in the lateral x and y directions, we used an additional vertical boundary with a weak repulsive potential [72]. In these simulations, only the position of the uranyl ion was constrained without controlling its immediate environment in the first coordination sphere. Well-tempered metadynamics simulations were performed using the “gaussian hills” approach [71]. Their initial height and width were 0.4 kJ·mol⁻¹ and 0.5 Å, respectively, and they were introduced every 250th time step with a bias factor of 15. All constrained MD simulations were performed for 40.0 ns with the same time step of 1.0 fs.

# 3. Results and Discussion

## 3.1. *Atomic Density Profiles*

To assess the effects of the magnitude and localization of the structural negative charge in the clay layers on the local structure of the uranyl surface complexes, we first compared the one-dimensional atomic density profiles of solution species at the clay interface. The closest to the surface peaks of U density are located at 4.10 Å from the illite surface for both versions of the ClayFF force field (Figure 1). For the smectite surface, the peak position is slightly different between the ClayFF-orig and ClayFF-MOH models. For ClayFF-MOH, it is located a little farther from the surface by ~0.1 Å. However, such a difference is probably due to minor differences in the local U speciation and can be neglected. The peak at ~4.1 Å corresponds to the nearest outer-sphere uranyl complex at the basal clay surface, which is characterized by hydrogen bonding between the uranyl-coordinating H₂O molecule/-OH group and the surface. The halfwidth of the peak (at ~4.1 Å) is larger for the smectite model with both versions of ClayFF, pointing to a higher mobility of uranyl complexes on the smectite surface. The higher integral peak intensity indicates higher adsorption of uranyl on the illite surface, due to the stronger electrostatic interaction with the higher-charge tetrahedrally substituted surface. Although it is known from experimental studies that uranyl adsorption is higher on smectites, note that in the current work, we investigated only uranyl adsorption at the external basal surfaces via ion-exchange mechanism, which is predominant at a pH < 4. Overall, higher adsorption of uranyl with ClayFF-MOH for both illite and smectite surfaces is predicted according to the U atomic density profiles. However, the average distance of adsorbed uranyl ions from the surface is similar for both versions of ClayFF.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Atomic density profiles of U atoms of uranyl ions as functions of their distance from the basal clay surface modeled with two versions of the ClayFF force field.</figcaption></figure>

There is also a distinct maximum of U density located at ~7 Å, indicating the formation of the second uranyl outer-sphere complex, farther from the surface. The calculated two-dimensional surface atomic density distributions within the 8 Å layer of the interfacial solution show that uranyl complexes are preferably adsorbed at the center of the ditrigonal cavities of the tetrahedral sheet for illite and smectite surfaces (Supplementary Materials Figures S1 and S2). For illite, uranyl is preferably adsorbed at the center of the ditrigonal cavities near substitutions (Figure S2). Such adsorption sites on the basal surface of clay minerals are known to be the most favorable for cations [73–75].

Under the conditions studied (i.e., considering only adsorption via ion-exchange mechanism on the basal surfaces at pH < 4), the adsorption of uranyl on the illite surface is stronger than on the smectite surface, which is seen from the atomic density profiles of the uranyl oxygen (Ou) atoms (Figure 2). The Ou profile on the illite surface has only two peaks located at 2.9 Å and at 5.2 Å. The distance between these peaks at 2.9 Å and at 5.2 Å is 2.3 Å, which is less than the full equilibrium length of a uranyl cation (~1.8 × 2 = 3.6 Å) [76]. Thus, uranyl is mainly adsorbed with an inclination of its O-U-O axis at about 40° to the illite surface, as predicted with both versions of ClayFF. A similar value of 36° was obtained earlier in the simulation of uranyl adsorption on a muscovite surface with ClayFF-orig [32]. The tilting of the uranyl ion on the illite surface could be explained by a stronger interaction of U with O atoms from the tetrahedral substitutions (basal O). Recent ab initio simulations showed the presence of weak chemical bonds between U and basal O atoms due to the overlap of uranium f- and oxygen p- electronic states [36].

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Atomic density profiles of O<sub>u</sub> atoms of uranyl as functions of their distance from the basal clay surface modeled with two versions of the ClayFF force field.</figcaption></figure>
According to a previous MD simulation of UO₂²⁺ sorption on the muscovite surface [32], the atomic density profile of U shows two distinct peaks at 2.0 Å and 3.25 Å, while in the presence of competing K⁺ ions, only a 3.25 Å peak was observed. Our results show only one outer-sphere uranyl complex with a U peak at 4.1 Å (Figure 1). The peak at 2 Å is absent in our results due to the presence of competing K⁺ and Ca²⁺ cations, in line with the previous study [32]. The difference of ~1 Å between 4.1 Å in our results and 3.25 Å in previously reported results for similar systems points to a strong effect of solution composition (ionic strength, nature, and concentration of the counterions and UO₂²⁺ concentration), which is not exactly the same in our studies. It is also interesting to note that while the distance from U to the surface differs significantly, the inclination angle of 40° is preserved. Therefore, our results are similar to those reported previously for uranyl adsorption on muscovite, although the distances are somewhat different due to differences in clay mineral structure, solution composition, etc. Specifically, we should point out that the amount of tetrahedral substitutions in the muscovite model studied earlier [32] is more than two times higher than in our current illite model, which obviously influences the results and makes any direct quantitative comparison difficult and inconclusive. The shape of the curves is similar for both versions of ClayFF, indicating only minor changes in the interfacial complex orientation. From this observation, we can conclude that the modification of the force field does not significantly affect the position of the complex at the illite surface, which is not the case for smectite.

At the smectite surface, the orientational effect is less pronounced, but the curve shape is somewhat different between the ClayFF-orig and ClayFF-MOH models. With ClayFF-orig, two smooth peaks of the surface-nearest O_u atoms at 2–3 Å are observed, which indicates the absence of a preferential uranyl orientation at the smectite basal surface. With ClayFF-MOH, the distribution of O_u is also smooth, but a more intense maximum at 2.9 Å is observed, pointing to the increased preference for inclined uranyl ions orientation. The increased fraction of inclined uranyl ions may be a consequence of the increased near-surface OH⁻ concentration in the smectite ClayFF-MOH model (Figure 3). This enables the formation of more hydrolyzed [UO₂OH]⁺ complexes in the smectite ClayFF-MOH model.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Atomic density profiles of O<sub>hw</sub> atoms (O<sub>hw</sub> is O atom of OH<sup>-</sup> in solution) as functions of their distance from the basal clay surface modeled with two versions of the ClayFF force field.</figcaption></figure>

Significant differences between the illite and smectite hydrated interfaces are observed in the shape of the O<sub>w</sub> density profile of the water molecules (Figure S3). At the illite surface, a small feature appeared before the main maximum, indicating the presence of a layer of stronger adsorbed H<sub>2</sub>O molecules. The center of the distribution coincides for both clay surfaces; minor peaks after the main one are closer to the surface of illite. A more structured profile of O<sub>w</sub> at the illite surface is clearly the result of higher layer charge and its localization on the tetrahedral sites of the TOT layer closer to the interface. No differences in the interfacial water structure are observed between the two ClayFF versions neither for the illite nor for the smectite surfaces (Figure S3).
## 3.2. Radial Distribution Functions

The local environment of the adsorbed uranyl ions was also quantitatively characterized through the analysis of the calculated radial distribution functions (RDFs) and running coordination numbers (CNs) of H₂O and OH⁻ species around the U atom. RDFs and CNs of U-Ow (Ow is O atom of H₂O) and U-Ohw (Ohw is O atom of OH⁻ in solution) pairs for both the smectite and illite simulations are shown in Figures 4 and 5. Additional RDF and CN functions of U, Ca, and K could be found in Supplementary Materials (Figures S5–S11). The calculated RDFs for water virtually coincide with previously calculated RDFs for UO₂²⁺ at the hydrated quartz interface [26]. Uranyl cation at the surface of clay minerals has a uniform coordination of ~4.5 water molecules at a distance of 2.50 Å. This value is very close to the parameters for outer-sphere uranyl complexes on smectites at low pHs obtained with EXAFS (2.38–2.50 Å) [17–25].

<figure><img src="image_7.png" /><figcaption>**Figure 4.** (a) Radial distribution functions for U-Ow pairs; (b) running coordination numbers for U-Ow pairs.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** (a) Radial distribution functions for U-Ohw pairs; (b) running coordination numbers for U-Ohw pairs.</figcaption></figure>

To the best of our knowledge, there are no experimental EXAFS studies of uranyl sorbed on illite clays due to (i) difficulties in obtaining illite samples without smectite impurities and (ii) the relatively low sorption capacity of illite towards UO₂²⁺, which complicates obtaining a good enough EXAFS spectrum. From this point of view, our theoretical study sheds new light on the process of uranyl adsorption on the illite surface. It is interesting that no major differences in U-H₂O local coordination were observed between the illite and smectite interfaces, as well as between the two versions of the force field. The coordination of water molecules to uranyl cation is similar with the ClayFF-orig and ClayFF-MOH force fields, in accordance with the observed close similarity of the water density profiles for both force field versions (Figure S3). Comparing the Ow profiles for
the illite and smectite surfaces, we clearly observed the more structured distribution of water on illite. However, the interfacial water structuring does not affect the coordination of water molecules near the uranyl ions adsorbed as outer-sphere complexes.

Changes in the local U coordination are notable in the analysis of U–O<sub>hw</sub> RDF and CN (Figure 5). The distance between the U and O atoms of the hydroxyl ion is 2.33 Å and it is the same for both versions of the force field and for both surfaces. When comparing the CNs (Figure 5b), it is clear that the modification of the force field affects the local U coordination. With ClayFF-orig, the local environment of the UO<sub>2</sub><sup>2+</sup> complexes is similar on the illite and smectite surfaces, and the CN for the OH<sup>−</sup> group equals ~0.6. According to these data, the time-averaged complex speciation could be represented as [UO<sub>2</sub>(H<sub>2</sub>O)<sub>4.3</sub>(OH)<sub>0.7</sub>]<sup>1.3+</sup> for both surfaces. The non-integer CN means that there are two fractions of uranyl complexes: hydrolyzed and non-hydrolyzed, and the increase in CN for the hydroxyl group means that the fraction of hydrolyzed complexes increases. In simulations with the ClayFF-MOH model, it became possible to differentiate the composition of the complexes between the illite and smectite surfaces. On the smectite surface, ~0.9 OH<sup>−</sup> anion is coordinated to each uranyl cation, meaning that hydrolyzed complexes dominate and uranyl is almost completely present as [UO<sub>2</sub>OH(H<sub>2</sub>O)<sub>x</sub>]<sup>+</sup> (x = 4–5). The time-averaged speciation is [UO<sub>2</sub>(H<sub>2</sub>O)<sub>4.1</sub>(OH)<sub>0.9</sub>]<sup>1.1+</sup>, which means that around 90% of species are present as hydrolyzed complexes. In the illite simulations with the ClayFF-MOH model, the fraction of hydrolyzed complexes noticeably decreases although the CN still remains around 0.6. The rest of the uranyl cations are present as non-hydrolyzed aqua complexes [UO<sub>2</sub>(H<sub>2</sub>O)<sub>5</sub>]<sup>2+</sup>. The time-averaged speciation in this case results in [UO<sub>2</sub>(H<sub>2</sub>O)<sub>4.4</sub>(OH)<sub>0.6</sub>]<sup>1.4+</sup> composition. It should be noted that in each case, the sum of the CNs of the ligands gives precisely five, the most common CN around the U equatorial sphere in solution.

The speciation of U sorbed at the clay basal surfaces was thought to be particularly independent of the clay mineral structure, namely the magnitude and localization of the negative charge in the TOT layer. Similar local environments at the interfaces with two end-members of illite–smectite clay were observed with the ClayFF-orig version of the force field. However, a more accurate parametrization of the structural hydroxyls in the ClayFF-MOH model highlights the different nature of the adsorbed uranyl complexes on the illite and smectite surfaces. We might also notice that despite the observed differences, hydrolyzed uranyl complexes are still predominant on both surfaces. Therefore, the uranyl behavior in models with two end-members suggests that in the intermediate cases of mixed-layer illite–smectite clay minerals, the U speciation would also be similar.

The RDF of U–Si pairs at the illite and smectite surfaces reveals that the main maximum is located at distances higher than 4 Å from the surface (Figure 6), indicating only a weak interaction of uranyl complexes with the surface in both cases. However, some differences in the fine structure of the main peak are observed. For the illite interface, the RDF has a maximum for U–Si and U–Al pairs due to the ambiguity in the location of the uranyl complexes above the substituted and non-substituted ditrigonal cavities with both versions of the ClayFF force field. For the ClayFF-MOH model, the peaks are more resolved, again highlighting the effect of a more accurate parametrization. For the smectite interface, where only octahedral substitutions are present deeper in the clay structure, the U–Si RDF shows a single broader peak for both versions of ClayFF, meaning that the effects of uranyl complex location above a specific surface site are less important.

Experimental EXAFS studies of uranyl sorbed on smectites as an inner-sphere edge complex do indicate the presence of Si/Al coordination shells at 3.1, 3.3, or 3.6 Å [19,21,23–25], which is also supported with DFT calculations of the clay edge complexes [37–41]. In our models, only adsorption at the external basal surfaces of clays is considered and therefore the comparison with those experimental results would be inconsistent. The implementation of the newer ClayFF-MOH force field version [44] virtually did not change the positions of the RDF peaks in comparison to the original ClayFF parametrization [42]. However, the intensities of the RDF peaks changed quite a lot, leading to a higher uranyl adsorption predicted with the ClayFF-MOH model. At the illite surface, for U–Al pairs with ClayFF-MOH, the peak intensity
is significantly higher, indicating that uranyl cations would be more preferably adsorbed at the substituted ditrigonal cavities. The double maximum of this distribution indicates the existence of two distinct surface configurations slightly off the center of the ditrigonal cavities. With the ClayFF-orig model, a lower intensity of the second RDF peak for U–Al pairs is observed. Generally, the RDF intensity of U–Al pairs is higher than the RDF intensity for the U–Si pairs, which again points to a higher uranyl adsorption at the substituted ditrigonal cavities as modeled with ClayFF-MOH.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** (a) Radial distribution functions for U–Si pairs from smectite-solution simulations; (b) radial distribution functions for U–Si/Al pairs from illite-solution simulations.</figcaption></figure>

### 3.3. *Free Energy of Uranyl Adsorption at Clay Surfaces*

The free energy of adsorption, Δ*A*, was calculated from the equilibrium NVT-ensemble MD trajectories using the metadynamics approach, as described in Section 2.3. It is important to note, however, that under ambient conditions, one can assume that Δ*A* ≈ Δ*G*, i.e., Gibbs free energy of adsorption is equal to the free energy from our metaD calculations [72]. The adsorption-free energy profiles for uranyl ions at the illite and smectite surfaces calculated using the ClayFF-orig and ClayFF-MOH models are shown in Figure 7. To be able to compare the results for different force fields and surfaces, in these constrained MD simulations, only the interaction of the uranyl ions themselves with the clay surfaces was controlled and analyzed, while the nearest presence of coordinating water molecules and, occasionally, of hydroxyl ions, was not controlled and not taken into account in the analysis. Therefore, the results presented below relate to UO₂²⁺ species without consideration of its instantaneous local coordination. At both clay surfaces, the first minimum is located at 4.4–4.5 Å from the surface. Small differences from the atomic density peak positions can be explained by the choice of the center of mass of the ditrigonal cavity as the main adsorption site and the absence of Si/Al substitutions in constrained MD simulations. Thus, the adsorption of uranyl in the unconstrained MD simulations occurs somewhat farther from the center of mass of the ditrigonal cavity, which is also noticeable on time-averaged atomic density distribution maps (Figures S1 and S2). The first minimum verifies that the outer-sphere basal surface complexes of uranyl are strongly preferable. The free energy minimum at the illite surface is at –14.15 kJ/mol with ClayFF-orig, which is a relatively high absolute value, suggesting more stable adsorption of uranyl on the illite basal surface. However, the newer ClayFF-MOH model predicts a much shallower free energy minimum at only –6.05 kJ/mol, thus indicating much weaker uranyl adsorption at the illite basal surface. This result seems to contradict our conclusions derived from the analysis of atomic density profiles (Figures 1 and 2) and radial distribution functions (Figure 6), which consistently indicated that ClayFF-MOH predicts higher adsorption of uranyl. The explanation of this contradiction may be in the fact that the constrained simulations of the adsorption-free energy profiles were only performed for uranyl cation without controlling its nearest environment, while the RDF and CN calculations (Figure 4) indicate that both hydrated UO₂²⁺ and hydrolyzed [UO₂OH]⁺ surface species can be important. Another possible explanation of the discrepancies between the constrained and unconstrained simulations might be that
in the constrained simulations, only adsorption above the unsubstituted ditrigonal cavities was considered. Further calculations will clarify the effects of the tetrahedral substitutions on the adsorption-free energy, as it was suggested earlier [72].

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Adsorption-free energy profiles for uranyl cation as functions of distance from basal clay surfaces for illite (red curves) and smectite (blue curves).</figcaption></figure>

For the smectite surface, the first free energy minimum is −9.55 kJ/mol with the ClayFF-orig model, but only −2.79 kJ/mol with the more recent ClayFF-MOH model. Again, the inconsistency of the atomic density and energy profiles is due to the presence of a significant fraction of hydrolyzed uranyl ions. As expected from simple electrostatic considerations, both models predict lower adsorption at the lower-charge smectite surface than at the higher-charge illite surface. This is consistent with the observed differences in the atomic density profiles (Figure 1).

The second free energy minima in Figure 7 corresponds to the second outer-sphere uranyl complexes farther away from the basal clay surfaces. For smectite, the second minimum is predicted at ~7.4 Å from the surface with both versions of ClayFF. The values of the minima are close, −1.31 kJ/mol and −0.41 kJ/mol with ClayFF-orig and ClayFF-MOH, respectively. This result is also consistent with the close values of the intensities for the uranium density profiles (Figure 1). The position of the second free energy minimum from the illite simulations is ~0.2 Å closer to the basal clay surface, located at 7.2 Å and 7.1 Å with the ClayFF-orig and ClayFF-MOH models, respectively. It is interesting to note that the energy value of −5.18 kJ/mol with ClayFF-MOH is close to the value of the first free energy minimum for this system. The earlier ClayFF-orig model predicts the second free energy minimum at −8.08 kJ/mol, which is significantly smaller by absolute value than the global minimum for this simulation (−14.15 kJ/mol). However, no differences in the intensities of the atomic density profiles around ~7 Å are observed, again pointing to a high importance of [UO₂OH]⁺ contribution.

It is worth emphasizing that for both the smectite and illite clay surfaces, ClayFF-MOH predicts the adsorption-free energy minima much less pronounced compared to the ClayFF-orig model. This can be explained by (i) the impact of the hydrolyzed complexes contributing to the calculated atomic density profiles but not taken into account in free energy calculations and (ii) the effect of stronger orientational constraints on the libration of structural hydroxyl ions within the ditrigonal cavities of the TOT clay structure with the ClayFF-MOH force field, which affects the general electrostatic strength of uranyl attraction to the surface.

4. Conclusions

The research on the quantitative macroscopic characterization of cation adsorption and retention in clay materials is quite substantial, but studies focused on the molecular mechanisms of these processes remain rather limited, despite growing interest from both fundamental geochemistry and from numerous engineering and environmental applications, including the technology of geological disposal and storage of radioactive waste [76].
In this work, we aimed to study, on the fundamental molecular scale, the effects of the magnitude and localization of the negative structural charge of clay TOT layers on the adsorption capacity for uranyl ions at hydrated basal clay surfaces and on the local structural environment of adsorbed uranyl ions. For this purpose, two new models of end-member clay structures in the illite-smectite mixed-layer series were constructed and investigated using the classical molecular dynamics computer simulation technique. The effect of the layer charge turned out to be not as significant as it was initially expected. Near both clay surfaces, uranyl ions are weakly adsorbed at similar distances in the form of outer-sphere non-hydrolyzed aqua complexes $[UO_2(H_2O)_5]^{2+}$ and partly hydrolyzed $[UO_2(H_2O)_{4-5}OH]^{+}$ with U-O<sub>w</sub> and U-O<sub>hw</sub> distances around 2.50 and 2.33 Å, respectively. Local uranyl coordinations evaluated from the calculated radial distribution functions revealed that the fraction of hydrolyzed species is predominant for both surfaces. This result suggests that in mixed-layer clay minerals, such as ISCz-1, local uranium coordination would also remain similar. The observed U-O<sub>w</sub> distances are very close to the results of EXAFS studies of outer-sphere uranyl complexes.

The other side of our work concerned the comparison of predicted uranyl adsorption made from the simulations using the original [42] and the more recently modified [44] ClayFF classical force field models. It is shown that the modified version, ClayFF-MOH [44], provides a more detailed description of the structure of uranyl surface complexes, while keeping many earlier results principally unchanged. The use of the modified force field allowed us to observe more clearly the changes in the adsorbed uranyl complex composition, which seemed impossible with ClayFF-orig. With ClayFF-MOH, we can also observe a decrease in the adsorption-free energy for $UO_2^{2+}$ cation over ditrigonal cavities without substitutions. As a result, the modified ClayFF-MOH model predicts a higher $UO_2^{2+}$ surface mobility. However, a detailed quantitative analysis of the adsorption-free energy of various uranyl complexes over various substituted sites at basal clay surfaces and, consequently, the particular contributions of all these factors to the overall $UO_2^{2+}$ retention and surface mobility, requires further systematic simulation studies, and will be reported separately.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/min14010109/s1, Figure S1: Time-averaged atomic density distributions of uranium at the basal surface of smectite; Figure S2: Time-averaged atomic density distributions of uranium on the basal surface of illite; Figure S3: Atomic density profiles of Ow atoms of H₂O molecules for both versions of ClayFF force field as functions of their distance from the basal clay surface; Figure S4: Atomic density profiles of HO atoms as functions of their distance from the basal clay surface modeled with two versions of the ClayFF force field; Figure S5: (a) Radial distribution functions for U-O<sub>b</sub> pairs; (b) Running coordination numbers for U-O<sub>b</sub> pairs (smectite model); Figure S6: (a) Radial distribution functions for Ca-O<sub>w</sub> pairs; (b) Running coordination numbers for Ca-O<sub>w</sub> pairs (smectite model); Figure S7: (a) Radial distribution functions for Ca-O<sub>hw</sub> pairs; (b) Running coordination numbers for Ca-O<sub>hw</sub> pairs (smectite model); Figure S8: (a) Radial distribution functions for U-O<sub>b</sub> pairs; (b) Running coordination numbers for U-O<sub>b</sub> pairs (illite model); Figure S9: (a) Radial distribution functions for U-O<sub>bts</sub> pairs; (b) Running coordination numbers for U-O<sub>bts</sub> pairs (illite model); Figure S10: (a) Radial distribution functions for K-O<sub>w</sub> pairs; (b) Running coordination numbers for K-O<sub>w</sub> pairs (illite model); Figure S11: (a) Radial distribution functions for K-O<sub>hw</sub> pairs; (b) Running coordination numbers for K-O<sub>hw</sub> pairs (illite model).

**Author Contributions:** Conceptualization, A.G.K. and I.E.V.; methodology, A.G.K.; investigation, E.V.T. and A.D.K.; data analysis A.D.K. and E.V.T.; writing—original draft preparation, A.D.K., E.V.T.; writing—review and editing, A.G.K. and I.E.V.; visualization, E.V.T.; supervision, A.G.K.; project administration, I.E.V.; funding acquisition, I.E.V. and A.G.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was financially supported by Russian Science Foundation, grant number 23-73-30006. A.G.K. also acknowledges the financial support of the industrial chair "Storage and Disposal of Radioactive Waste" at the IMT Atlantique, funded by ANDRA, Orano, and EDF. The APC was funded by publication vouchers from previous MDPI reviewing activity.
**Data Availability Statement:** All data are available from the authors upon request.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Dohrmann, R.; Kaufhold, S.; Lundqvist, B. The role of clays for safe storage of nuclear waste. In *Developments in Clay Science*, 2nd ed.; Elsevier Ltd.: Amsterdam, The Netherlands, 2013; Volume 5. [CrossRef]

2. Sellin, P.; Leupin, O.X. The use of clay as an engineered barrier in radioactive-waste management—A review. *Clays Clay Miner.* **2014**, *61*, 477–498. [CrossRef]

3. Grambow, B. Geological disposal of radioactive waste in clay. *Elements* **2016**, *12*, 239–245. [CrossRef]

4. Norris, S. Radioactive waste confinement: Clays in natural and engineered barriers—Introduction. *Geol. Soc. Spec. Publ.* **2017**, *443*, 1–8. [CrossRef]

5. Belousov, P.E.; Krupskaya, V.V. Bentonite clays of Russia and neighboring countries. *Georesources* **2019**, *21*, 79–90. [CrossRef]

6. Eisenhour, D.D.; Brown, R.K. Bentonite and Its Impact on Modern Life. *Elements* **2009**, *5*, 83–88. [CrossRef]

7. Krupskaya, V.V.; Biryukov, D.; Belousov, P.; Lekhov, V.; Romanchuk, A.Y.; Kalmykov, S. The use of natural clay materials to increase the nuclear and radiation safety level of nuclear legacy facilities. *Radioact. Waste* **2018**, *2*, 24–34.

8. Laverov, N.P.; Yudintsev, S.; Kochkin, B.; Malkovsky, V. The russian strategy of using crystalline rock as a repository for nuclear waste. *Elements* **2016**, *12*, 253–256. [CrossRef]

9. Wang, C.; Mushkin, V.F.; Khan, V.A.; Panamareva, A.N. A review of the migration of radioactive elements in clay minerals in the context of nuclear waste storage. *J. Radioanal. Nucl. Chem.* **2022**, *331*, 3401–3426. [CrossRef]

10. Ma, Z.; Pathegama Gamage, R.; Rathnaweera, T.; Kong, L. Review of application of molecular dynamic simulations in geological high-level radioactive waste disposal. *Appl. Clay Sci.* **2019**, *168*, 436–449. [CrossRef]

11. IAEA. *Status and Trends in Spent Fuel and Radioactive Waste Management*; IAEA Nuclear Energy Series; International Atomic Energy Agency: Vienna, Austria, 2022; 88p. Available online: https://www.iaea.org/publications/14739/status-and-trends-in-spent-fuel-and-radioactive-waste-management (accessed on 10 January 2024).

12. Missana, T.; Alonso, U.; Garcia-Gutierrez, M.; Albarran, N.; Lopez, T. Experimental Study and Modeling of Uranium (VI) Sorption onto a Spanish Smectite. *MRS Proc.* **2008**, *1124*, 705. [CrossRef]

13. Bachmaf, S.; Merkel, B.J. Sorption of uranium(VI) at the clay mineral-water interface. *Environ. Earth Sci.* **2010**, *63*, 925–934. [CrossRef]

14. Villa-Alfageme, M.; Hurtado, S.; El Mrabet, S.; Pazos, M.C.; Castro, M.A.; Alba, M.D. Uranium immobilization by FEBEX bentonite and steel barriers in hydrothermal conditions. *Chem. Eng. J.* **2015**, *269*, 279–287. [CrossRef]

15. Fan, Q.; Li, P.; Pan, D. Radionuclides sorption on typical clay minerals: Modeling and spectroscopies. In *Interface Science and Technology*, 1st ed.; Elsevier Ltd.: Amsterdam, The Netherlands, 2019; Volume 29. [CrossRef]

16. Philipp, T.; Shams Aldin Azzam, S.; Rossberg, A.; Huittinen, N.; Schmeide, K.; Stumpf, T. U(VI)sorption on Ca-bentonite at (hyper)alkaline conditions—Spectroscopic investigations of retention mechanisms. *Sci. Total Environ.* **2019**, *676*, 469–481. [CrossRef] [PubMed]

17. Thompson, H.A.; Brown, G.E.; Parks, G.A. XAFS spectroscopic study of uranyl coordination in solids and aqueous solution. *Am. Mineral.* **1997**, *82*, 483–496. [CrossRef]

18. Sylwester, E.R.; Allen, P.G.; Zhao, P.; Viani, B.E. Interactions of uranium and neptunium with cementitious materials studied by XAFS. *Mater. Res. Soc. Symp. Proc.* **2000**, *608*, 307–312. [CrossRef]

19. Hennig, C.; Reich, T.; Dähn, R.; Scheidegger, A.M. Structure of uranium sorption complexes at montmorillonite edge sites. *Radiochim. Acta* **2002**, *90*, 653–657. [CrossRef]

20. Catalano, J.G.; Brown, G.E. Uranyl adsorption onto montmorillonite: Evaluation of binding sites and carbonate complexation. *Geochim. Cosmochim. Acta* **2005**, *69*, 2995–3005. [CrossRef]

21. Walter, M.; Arnold, T.; Geipel, G.; Scheinost, A.; Bernhard, G. An EXAFS and TRLFS investigation on uranium(VI) sorption to pristine and leached albite surfaces. *J. Colloid Interface Sci.* **2005**, *282*, 293–305. [CrossRef]

22. Křepelová, A.; Reich, T.; Sachs, S.; Drebert, J.; Bernhard, G. Structural characterization of U(VI) surface complexes on kaolinite in the presence of humic acid using EXAFS spectroscopy. *J. Colloid Interface Sci.* **2008**, *319*, 40–47. [CrossRef]

23. Schlegel, M.L.; Descostes, M. Uranium uptake by hectorite and montmorillonite: A solution chemistry and polarized EXAFS study. *Environ. Sci. Technol.* **2009**, *43*, 8593–8598. [CrossRef]

24. Fernandes, M.M.; Baeyens, B.; Dähn, R.; Scheinost, A.; Bradbury, M. U(VI) sorption on montmorillonite in the absence and presence of carbonate: A macroscopic and microscopic study. *Geochim. Cosmochim. Acta* **2012**, *93*, 262–277. [CrossRef]

25. Troyer, L.D.; Maillot, F.; Wang, Z.; Wang, Z.; Mehta, V.S.; Giammar, D.E.; Catalano, J.G. Effect of phosphate on U(VI) sorption to montmorillonite: Ternary complexation and precipitation barriers. *Geochim. Cosmochim. Acta* **2016**, *175*, 86–99. [CrossRef]

26. Greathouse, J.A.; O'Brien, R.J.; Bemis, G.; Pabalan, R.T. Molecular dynamics study of aqueous uranyl interactions with quartz (010). *J. Phys. Chem. B* **2002**, *106*, 1646–1655. [CrossRef]

27. Greathouse, J.A.; Cygan, R.T. Molecular dynamics simulation of uranyl(VI) adsorption equilibria onto an external montmorillonite surface. *Phys. Chem. Chem. Phys.* **2005**, *7*, 3580–3586. [CrossRef] [PubMed]
28. Greathouse, J.A.; Cygan, R.T. Water structure and aqueous uranyl(VI) adsorption equilibria onto external surfaces of beidellite, montmorillonite, and pyrophyllite: Results from molecular simulations. *Environ. Sci. Technol.* **2006**, *40*, 3865–3871. [CrossRef]

29. Kerisit, S.; Liu, C. Molecular simulation of the diffusion of uranyl carbonate species in aqueous solution. *Geochim. Cosmochim. Acta* **2010**, *74*, 4937–4952. [CrossRef]

30. Kerisit, S.; Liu, C. Diffusion and adsorption of uranyl carbonate species in nanosized mineral fractures. *Environ. Sci. Technol.* **2012**, *46*, 1632–1640. [CrossRef]

31. Martorell, B.; Kremleva, A.; Krüger, S.; Rösch, N. Density functional model study of uranyl adsorption on the solvated (001) surface of kaolinite. *J. Phys. Chem. C* **2010**, *114*, 13287–13294. [CrossRef]

32. Teich-McGoldrick, S.L.; Greathouse, J.A.; Cygan, R.T. Molecular dynamics simulations of uranyl adsorption and structure on the basal surface of muscovite. *Mol. Simul.* **2014**, *40*, 610–617. [CrossRef]

33. Kalinichev, A.G.; Loganathan, N.; Ngouana-Wakou, B.F.; Chen, Z. Interaction of Ions with Hydrated Clay Surfaces: Computational Molecular Modeling for Nuclear Waste Disposal Applications. *Procedia Earth Planet. Sci.* **2017**, *17*, 566–569. [CrossRef]

34. Zhang, J.; Mallants, D.; Brady, P.V. Molecular dynamics study of uranyl adsorption from aqueous solution to smectite. *Appl. Clay Sci.* **2022**, *218*, 106361. [CrossRef]

35. Pérez-Conesa, S.; Martínez, J.M.; Marcos, E.S. Hydration and Diffusion Mechanism of Uranyl in Montmorillonite Clay: Molecular Dynamics Using an Ab Initio Potential. *J. Phys. Chem. C* **2017**, *121*, 27437–27444. [CrossRef]

36. Li, N.; Wang, J.; Wang, J.; Wang, Y.; Fu, Y.; Zhao, J. Atom-scale understanding the adsorption mechanism of uranyl in the interlayer of montmorillonite: Insight from DFT + U calculation. *Appl. Surf. Sci.* **2023**, *612*, 155910. [CrossRef]

37. Kremleva, A.; Krüger, S.; Rösch, N. Uranyl adsorption at (010) edge surfaces of kaolinite: A density functional study. *Geochim. Cosmochim. Acta* **2011**, *75*, 706–718. [CrossRef]

38. Kremleva, A.; Martorell, B.; Krüger, S.; Rösch, N. Uranyl adsorption on solvated edge surfaces of pyrophyllite: A DFT model study. *Phys. Chem. Chem. Phys.* **2012**, *14*, 5815–5823. [CrossRef]

39. Kremleva, A.; Krüger, S.; Rösch, N. Uranyl adsorption at solvated edge surfaces of 2:1 smectites. A density functional study. *Phys. Chem. Chem. Phys.* **2015**, *17*, 13757–13768. [CrossRef]

40. Kremleva, A.; Krüger, S.; Rösch, N. Toward a reliable energetics of adsorption at solvated mineral surfaces: A computational study of uranyl(VI) on 2:1 clay minerals. *J. Phys. Chem. C* **2016**, *120*, 324–335. [CrossRef]

41. Zhang, C.; Liu, X.; Tinnacher, R.M.; Tournassat, C. Mechanistic understanding of uranyl ion complexation on montmorillonite edges: A combined first-principles molecular dynamics-surface complexation modeling approach. *Environ. Sci. Technol.* **2018**, *52*, 8501–8509. [CrossRef]

42. Cygan, R.T.; Liang, J.J.; Kalinichev, A.G. Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field. *J. Phys. Chem. B* **2004**, *108*, 1255–1266. [CrossRef]

43. Cygan, R.T.; Greathouse, J.A.; Heinz, H.; Kalinichev, A.G. Molecular models and simulations of layered materials. *J. Mater. Chem.* **2009**, *19*, 2470–2481. [CrossRef]

44. Cygan, R.T.; Greathouse, J.A.; Kalinichev, A.G. Advances in Clayff molecular simulation of layered and nanoporous materials and their aqueous interfaces. *J. Phys. Chem. C* **2021**, *125*, 17573–17589. [CrossRef]

45. Tararushkin, E.V.; Pisarev, V.V.; Kalinichev, A.G. Equation of state, compressibility, and vibrational properties of brucite over wide pressure and temperature ranges: Atomistic computer simulations with the modified ClayFF classical force field. *Minerals* **2023**, *13*, 408. [CrossRef]

46. Tararushkin, E.V.; Smirnov, G.S.; Kalinichev, A.G. Structure and properties of water in a new model of the 10-Å Phase: Classical and ab initio atomistic computational modeling. *Minerals* **2023**, *13*, 1018. [CrossRef]

47. Kalinichev, A.G. Atomistic modeling of clays and related nanoporous materials with ClayFF force field. In *Computational Modeling in Clay Mineralogy*; Sainz-Díaz, C.I., Ed.; Association Internationale pour l'Etude des Argiles (AIPEA): Bari, Italy, 2021; Volume 3, pp. 17–52. [CrossRef]

48. Martin, R.T.; Bailey, S.W.; Eberl, D.D.; Fanning, D.S.; Guggenheim, S.; Kodama, H.; Pevear, D.R.; Środoń, J.; Wicks, F.J. Report of the clay minerals society nomenclature committee: Revised classification of clay materials. *Clays Clay Miner.* **1991**, *39*, 333–335. [CrossRef]

49. Drits, V.A.; Kossowskaya, A.G. Clay minerals: Smectites, mixed-layer silicates. *Acad. Sci. USSR Trans.* **1990**, *446*, 214.

50. Altaner, S.P.; Ylagan, R.F. Comparison of structural models of mixed-layer illite/smectite and reaction mechanisms of smectite illitization. *Clays Clay Miner.* **1997**, *45*, 517–533. [CrossRef]

51. Velde, B.; Meunier, A. *The Origin of Clay Minerals in Soils and Weathered Rocks*; Springer Science and Business Media LLC: Dordrecht, GX, The Netherlands, 2008; ISBN 9783540756330.

52. Drits, V.A.; Sakharov, B.A. *X-ray Structural Analysis of Mixed-Layer Minerals*; Nauka Publishing: Moscow, Russia, 1976; 256p. (In Russian)

53. Nadeau, P.H.; Wilson, M.J.; McHardy, W.J.; Tait, J.M. Interstratified clays as fundamental particles. *Science* **1984**, *225*, 923–925. [CrossRef]

54. Srodon, J.; Zeelmaekers, E.; Derkowski, A. The charge of component layers of illite-smectite in bentonites and the nature of end-member illite. *Clays Clay Miner.* **2009**, *57*, 649–671. [CrossRef]

55. Drits, V.A.; Zviagina, B.B.; McCarty, D.K.; Salyn, A.L. Factors responsible for crystal-chemical variations in the solid solutions from illite to aluminoceladonite and from glauconite to celadonite. *Am. Mineral.* **2010**, *95*, 348–361. [CrossRef]
56. Viani, A.; Gualtieri, A.; Artioli, G. The nature of disorder in montmorillonite by simulation of X-ray powder patterns. *Am. Mineral.* **2002**, *87*, 966–975. [CrossRef]

57. Loewenstein, W. The distribution of aluminum in the tetrahedra of silicates and aluminates. *Am. Mineral.* **1954**, *39*, 92–96.

58. Chipera, S.J.; Bish, D.L. Baseline studies of the clay minerals society source clays: Introduction. *Clays Clay Miner.* **2001**, *49*, 398–409. [CrossRef]

59. Greathouse, J.A.; Durkin, J.S.; Larentzos, J.P.; Cygan, R.T. Implementation of a Morse potential to model hydroxyl behaviour in phyllosilicates. *J. Chem. Phys.* **2009**, *130*, 134713. [CrossRef]

60. Androniuk, I.; Kalinichev, A.G. Molecular dynamics simulation of the interaction of uranium (VI) with the C–S–H phase of cement in the presence of gluconate. *Appl. Geochem.* **2020**, *113*, 104496. [CrossRef]

61. Berendsen, H.; Grigera, J.; Straatsma, T. The missing term in effective pair potentials. *J. Phys. Chem.* **1987**, *91*, 6269–6271. [CrossRef]

62. Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. *J. Comput. Phys.* **1995**, *117*, 1–19. [CrossRef]

63. Thompson, A.P.; Aktulga, H.M.; Berger, R.; Bolintineanu, D.S.; Brown, W.M.; Crozier, P.S.; in't Veld, P.J.; Kohlmeyer, A.; Moore, S.G.; Nguyen, T.D.; et al. LAMMPS—A flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. *Comp. Phys. Comm.* **2022**, *271*, 108171. [CrossRef]

64. Allen, M.P.; Tildesley, D.J. *Computer Simulation of Liquids*, 2nd ed.; Oxford University Press: New York, NY, USA, 2017; 626p. [CrossRef]

65. Frenkel, D.; Smit, B. *Understanding Molecular Simulation: From Algorithms to Applications*, 2nd ed.; Academic Press: San Diego, CA, USA, 2002; 638p.

66. Brown, W.M.; Kohlmeyer, A.; Plimpton, S.J.; Tharrington, A.N. Implementing molecular dynamics on hybrid high performance computers—Particle-particle particle-mesh. *Comp. Phys. Comm.* **2012**, *183*, 449–459. [CrossRef]

67. Brown, W.M.; Wang, P.; Plimpton, S.J.; Tharrington, A.N. Implementing molecular dynamics on hybrid high performance computers—Short range forces. *Comp. Phys. Comm.* **2011**, *182*, 898–911. [CrossRef]

68. Swope, W.C.; Andersen, H.C.; Berens, P.H.; Wilson, K.R. A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters. *J. Chem. Phys.* **1982**, *76*, 637–649. [CrossRef]

69. Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. *J. Chem. Phys.* **1984**, *81*, 511–519. [CrossRef]

70. Hoover, W.G. Canonical dynamics: Equilibrium phase-space distributions. *Phys. Rev. A* **1985**, *31*, 1695–1697. [CrossRef] [PubMed]

71. Tribello, G.A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. *Comput. Phys. Commun.* **2014**, *185*, 604–613. [CrossRef]

72. Loganathan, N.; Kalinichev, A.G. Quantifying the mechanisms of site-specific ion exchange at an inhomogeneously charged surface: Case of $Cs^+/K^+$ on hydrated muscovite mica. *J. Phys. Chem. C* **2017**, *121*, 7829–7836. [CrossRef]

73. Liu, X.D.; Lu, X.C.; Wang, R.C.; Zhou, H.Q. Effects of layer-charge distribution on the thermodynamic and microscopic properties of Cs-smectite. *Geochim. Cosmochim. Acta* **2008**, *72*, 1837–1847. [CrossRef]

74. Meleshyn, A. Adsorption of $Sr^{2+}$ and $Ba^{2+}$ at the cleaved mica-water interface: Free energy profiles and interfacial structure. *Geochim. Cosmochim. Acta* **2010**, *74*, 1485–1497. [CrossRef]

75. Krot, A.; Vlasova, I.; Trigub, A.; Averin, A.; Yapaskurt, V.; Kalmykov, S. From EXAFS of reference compounds to U(VI) speciation in contaminated environments. *J. Synchrotron Radiat.* **2022**, *29*, 303–314. [CrossRef]

76. Liu, X.; Tournassat, C.; Grangeon, S.; Kalinichev, A.G.; Takahashi, Y.; Marques Fernandes, M. Molecular-level understanding of metal ion retention in clay-rich materials. *Nat. Rev. Earth Environ.* **2022**, *3*, 461–476. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.