Probing atomic-scale processes at the ferrihydrite-water interface with reactive molecular dynamics

Abstract

Interfacial processes involving metal (oxyhydr)oxide phases are important for the mobility and bioavailability of nutrients and contaminants in soils, sediments, and water. Consequently, these processes influence ecosystem health and functioning, and have shaped the biological and environmental co-evolution of Earth over geologic time. Here we employ reactive molecular dynamics simulations, supported by synchrotron X-ray spectroscopy to study the molecular-scale interfacial processes that influence surface complexation in ferrihydrite-water systems containing aqueous ${\text {MoO}_4}^{2-}$. We validate the utility of this approach by calculating surface complexation models directly from simulations. The reactive force-field captures the realistic dynamics of surface restructuring, surface charge equilibration, and the evolution of the interfacial water hydrogen bond network in response to adsorption and proton transfer. We find that upon hydration and adsorption, ferrihydrite restructures into a more disordered phase through surface charge equilibration, as revealed by simulations and high-resolution X-ray diffraction. We observed how this restructuring leads to a different interfacial hydrogen bond network compared to bulk water by monitoring water dynamics. Using umbrella sampling, we constructed the free energy landscape of aqueous ${\text {MoO}_4}^{2-}$ adsorption at various concentrations and the deprotonation of the ferrihydrite surface. The results demonstrate excellent agreement with the values reported by experimental surface complexation models. These findings are important as reactive molecular dynamics opens new avenues to study mineral-water interfaces, enriching and refining surface complexation models beyond their foundational assumptions.

Graphic Abstract

Introduction

The mobility and bioavailability of nutrient and contaminant elements in soils, sediments and water are strongly influenced by atomic-scale processes at mineral-water interfaces. Surface adsorption or precipitation reactions can promote element attenuation, while desorption and dissolution reactions may drive element release [1,2,3]. Ferrihydrite (Fh) and other iron (oxyhydr)oxides are often important controls on aqueous concentrations of oxyanion-forming elements, which can have important implications for ecosystem health and functioning [4,5,6]. For example, molybdate (${\text {MoO}_4}^{2-}$) adsorption onto Fh surfaces has implications for biosynthesis of Mo-enzymes central to global biogeochemical cycles that have shaped Earth’s biological and environmental coevolution.

In aqueous systems, Fh typically bears a non-zero surface charge originating from the ionization of the surface hydroxyl (–OH) groups and the adsorption of dissolved ions [7, 8]. Oxyanions adsorption strongly influences physical, chemical, and structural properties of the Fh surface [9, 10], primarily due to the induced re-organizing of water molecules altering phase stability and its transformation pathways [11]. Furthermore, Fh surface charge modification through adsorption not only leads to a reorientation of interfacial species, affecting both their adsorption kinetics and diffusion across the interface but also alters the hydrogen bond hydrogen bond network at the interface. This alteration in the HB network additionally influences the acidity of protic species, which include any species capable of donating a hydrogen ion ($\text {H}^+$) [8, 12, 13].

Past studies of charged interfaces in water experiencing adsorption events have commonly relied on empirical studies, supported by various spectroscopic techniques and surface complexation models (SCM), while some recent studies have used more advanced methodologies [7, 14,15,16,17,18,19]. SCM account for surface charge along with solute-surface adsorption complex equilibrium constants to fit a model of surface complexes to batch adsorption experiments data [20, 21]. The models rely upon parameters to describe the equilibrium between aqueous species and mineral surface sites [22, 23] using parameters such as surface site density, surface potential, intrinsic acid–base constants, and apparent complexation constants, determined with the assistance of non-linear parameter estimation codes [23,24,25].

While SCM have proven valuable, they require fitting a large amount of data under the presumption that surfaces are static and homogeneous [25,26,27,28], and simplifying the effects of water chemistry [15, 29,30,31]. Reactive molecular dynamics (MD) simulations offer a complementary perspective by providing a realistic depiction of fluid-solid interactions. Reactive force fields such as ReaxFF allow one to simulate molecular dynamics while incorporating bond-breaking and formation events. They are more computationally efficient than ab initio MD [19, 32, 33] and density functional theory (DFT) methods for the Mo(VI)–$\text {FeO}_x$ interface [19, 34, 35]. ReaxFF force fields are trained using large data sets obtained from DFT calculations. The accuracy of ReaxFF force fields has been evaluated numerous time in the past by examining the capability of the force field to reproduce experimental data and DFT test sets for different types of systems including bulk solids, liquids and interfacial systems [36,37,38]. Reactive force-fields have proven successful in modeling solid-water interfaces [39], accurately describing the restructuring of clay surfaces and interlayers in water consistent with neutron scattering and nuclear magnetic resonance spectroscopy (NMR) data [40]. Furthermore, reactive force fields have been utilized to examine the reactions of water with solid surfaces. This includes describing reactivity trends of various crystal surfaces towards water dissociation in TiO₂ [41], and studying water dissociation at the surfaces of calcium silicate hydrate micropores and silica surfaces, successfully reproducing experimental properties, such as dipole moment, diffusion coefficients, and proton-transfer pathways similar to the Grotthuss mechanism [42, 43].

In this work we have used reactive MD paired with synchrotron X-ray spectroscopy to establish a model that can describe Fh-water interface with adsorbed aqueous ${\text {MoO}_4}^{2-}$. Our approach effectively integrates experimental findings with theoretical models [44]; we aim to (I) explore the effects of hydration and adsorption on the restructuring of the Fh by examining interfacial properties such as surface charge and interfacial HB network, (II) evaluate the accuracy of the model using synchrotron X-ray spectroscopy to compare the structure of solid in simulations with experimental data, and (III) directly calculating SCM constants from simulations using concentration-dependent free energy profiles for ${\text {MoO}_4}^{2-}$ adsorption and surface deprotonation (acidity). It is important to note that our work supplements the numerous previous studies on molybdate adsorption on Al and Fe oxides by including a quantitative description of outer-sphere complexes.

The choice of molybdenum in its aqueous molybdate form (${\text{MoO}_4}^{2-}$) as a model oxyanion is motivated by its dual role as an essential micronutrient and a pollutant [45, 46]. Bioavailable molybdenum is crucial in the function of oxidoreductase enzymes within the biogeochemical cycles of N, C, and S. Furthermore, mobile aqueous molybdenum species negatively impact water quality, underscoring its significant environmental and biological implications.

Materials and methods

Sample preparation

Synthesis of 2-line Fh followed the method of Schwertmann and Cornell [47], where 500 mL of a 0.2 M Fe(NO₃)₃⋅9H₂O solution was titrated to approximately pH 7.5 using 1 M KOH. The base addition was conducted over approximately 30 min and the resulting solution was stirred for approximately 1 h. Aliquots of the suspension were transferred into 50 mL conical tubes, which were centrifuged at 3000 rpm for 5 min. The supernatant was decanted and 40 mL ultrapure water was added and mixed. The centrifugation process and ultrapure water addition was repeated for a total of three times to remove residual salts. Solids were flash frozen and freeze dried, and powder X-ray diffraction confirmed phase identity.

Modeling setup

Fh slabs and surfaces were generated based on an initial structure reported by Michel et al. [48]. This structure served as the starting point for reactive MD simulations, which permits unrestrained restructuring in response to evolving chemical environment. The Fh supercell was terminated by oxygen atoms of ($1\bar{1}0$) plane, which accounts for 67–85% of the total surface area of Fh [49], to emulate geochemical conditions, using Pymatgen [50]. Symmetric slab with zero or no dipole moment was then selected for the preparation of $3.6 \times 3.5 \times 2.9\,\text {nm}^{3}$ Fh supercell, with the composition of $\text {Fe}_{960}$H₁₆₄O₁₇₂₈. Dangling surface oxygens were protonated by hydrogen atoms based on their acidity [51].

Reactive MD simulations

All simulations were carried out with LAMMPS package [52] using the ReaxFF pair_style as implemented in LAMMPS [53] ReaxFF parameters for Fh were obtained from the study by Aryanpour et al. [36], that accounts for Fe d-orbital electrons correlations for the training of the force field. The parameters for Mo [54], Na [55], and Na–Fe cross-terms [56] were adopted from reported studies.

The data files for simulations were prepared using the VMD TopoTools plugin and Moltemplate [57, 58]. Except in the microcanonical ensemble (NVE—constant number of atoms, volume, and energy) where a Berendsen thermostat was utilized, in all other cases Nosé–Hoover thermostat and barostat with damping constants of 100 fs, and 1000 fs were used to control the temperature and pressure, respectively. In all simulations, a time step of 0.25 fs was used to integrate Newton’s equation of motion by the velocity-Verlet algorithm. The charge polarization effects of all atoms were included in the simulation by calculating geometry-dependent atomic charges at every time step using the charge equilibration (QEq) approach with a cutoff radii 10.0 Å for all interactions [59].

The initial configuration of the slab was first minimized using the steepest descent algorithm while stopping tolerances for energy and force were set to $10^{-4}$ and $10^{-6}$ (kcal mol⁻¹ Å⁻¹), respectively. To prevent melting, the temperature of the system was gradually increased in a series of iterations in the canonical ensemble (NVT—constant number of atoms, volume, and temperature), starting with 10,000 steps of 0.1 Å maximum distance-limited NVT at 1 K (damping constant 1 fs), followed by NVT stages consisting of 10,000 steps at 1 K (damping constant 1 fs) and a 10,000 steps temperature ramp from 1 to 300 K (damping constant 100 fs). Finally, the slab was relaxed for 1 ns in the anisotropic NPT ensemble (constant number of atoms, pressure, and temperature) at 298 K and 1 atm in x and y directions in a simulation cell with the dimensions of approximately $4 \times 4 \times 6.7\,\text {nm}^{3}$ including a 4 nm gap in the z direction, in which the volume was subsequently filled with equilibrated aqueous solutions.

The initial configuration of solutions was prepared by randomly placing aqueous ions within the 3.0 Å region of the middle of the bulk aqueous solution to yield simulations of Fh and solutions containing different numbers of Na₂MoO₄ (Mo–n, n = 0, 1, 2, 3, 6, 8), corresponding to [Na₂MoO₄] = 0, 0.025, 0.05, 0.075, 0.15, and 0.2 (mol $\text {L}^{-1}$) with a water density of 33 (particle nm⁻³) using Packmol [60]. In the equilibration stage, solutions were first simulated while turning off ReaxFF mechanism. This was achieved by treating water O–H and molybdate Mo–O bonds as classical harmonic springs with a constant of 50 (${\text{kcal mol}}^{-1}$ Å⁻²). The spring constant was selected considering the effectiveness of pressure/temperature control and the stability of simulations with the used time step [61]. The equilibration of aqueous systems was performed for 0.2 ns while constraining O–H bonds in water and Mo–O bonds in molybdate, in stages comprising of 20 ps of 0.1 Å maximum distance-limited NVT at 298 K (damping constant 100 fs), 20 ps in the NVT 298 K (damping constant 50 fs), 60 ps in the NVT at 298 K (damping constant 50 fs), and 100 ps in the NPT at 298 K and 1 atm in z direction. The relaxed Fh slab and equilibrated aqueous systems were merged for the simulations of Fh-water-molybdate systems. Three independent simulations of each system were modeled to capture statistical variations for 1 ns in the NVT ensemble (with constrained water O–H and molybdate Mo–O), continued by 2 ns of reactive simulation in the NVT ensemble. The reactive simulations for each system were combined and used for further analysis using in house python codes based on MDAnalysis package [62, 63].

Umbrella sampling

Umbrella sampling is a technique that enables sampling of high-energy states, also known as rare events, that are otherwise not feasible for direct sampling [64] during simulations. It acts by enforcing an external harmonic restraint along a reaction coordinate, in a series of windows, represented by a collective variable, which, in this case, is the distance of Mo from the surface. This free energy along the collective variable is also referred to as the potential-of-mean-force (PMF).

We used Colvars module interfaced with LAMMPS to perform umbrella sampling [65]. A series of 20 initial configurations (windows) were generated from 0 to 20 Å from the surface in the z direction (width = 1 Å). A constraining force constant of 5 Kcal mol⁻¹ Å⁻² was used, while allowing the selected Mo atom to move in x and y directions. The selected force constant resulted in sufficient overlap of the position histograms (shown in Figure S20), allowing a smooth reconstruction of the PMF. Each window was run in NVT ensemble for at least 3 ns. The PMFs were obtained using an algorithm known as Weighted Histogram Analysis Method (WHAM) which is used to remove the biasing potential and reconstruct the PMF with error analysis performed using the Monte Carlo bootstrap method [66, 67].

Synchrotron X-ray measurements

All X-ray measurements were conducted at room temperature at Brockhouse Diffraction Sector (BXDS) at Canadian Light Source (CLS). The high-resolution X-ray diffraction (XRD) data were collected at Low Energy Wiggler Beamline (WLE) at the CLS. The measurements were conducted by using a Kapton capillary filled with unaged solid particles or particle dispersions in de-ionized water with no background electrolytes, with ${\text {[Mo]}_{\text{aq}}}$ concentration of 0.12 M. The capillary was mounted and spun during the measurements to improve uniformity. A fixed photon energy $\lambda$ = 0.819 Å with a Si (111) monochromator and Dectris Mythen2 X series 1 K detector system were used.

Pair distribution function (PDF) measurements were performed at the High Energy Wiggler (HEW) beamline (Energy range 20–94 keV; beam size 2–10 mm $\times$ $\sim$ 100 mm (W $\times$ H); energy resolution 1 $\times$ $10^{-3}$ $\Delta d/d$) at 54 keV (0.2265 Å) with a crystal Si (111) monochromator. The samples were loaded in Kapton capillary tubes and were spun during the measurements to improve uniformity. A fixed photon energy $\lambda$ = 0.819 Å with a Si (111) monochromator and Dectris Mythen2 X series 1K detector system were used. The data processing was calibrated using a Ni standard. The data was integrated from 2D to 1D using GSAS-II [68], and the PDFs were generated by the Fourier transform of total scattering structure factor S(Q) with a $Q_{max}$ of 22 Å⁻¹. The PDF of containers were measure for background signal removal.

Calculations of X-ray scattering

For XRD calculation, the thermal disorder was taken into account by averaging snapshots from MD trajectories separated by a 1 ps time step. To further enhance the interfacial signal, and emulate the experimental setup, the snapshots were limited to the interfacial region found by identification of truly interfacial molecules (ITIM) algorithm. XRD patterns were simulated at $\lambda$ = 0.819 Å with parameters for Fe(III) and Mo(VI) [69] using XrDebye module of Atomic Simulation Environment (ASE) package [70]. The PDF plots of selected MD snapshots were calculated from atomic coordinates using the Debye scattering equation implemented in Diffpy-CMI as DebyeCalculator class [71]. Differential pair distribution functions (d-PDFs) were calculated by subtracting a scaled reference PDF, Fh-water mixture from the PDF of Fh-${\text {[Mo]}_{\text{aq}}}$.

Results and discussion

Synchrotron X-ray spectroscopy and reactive molecular dynamics

The high-resolution XRD experiments and ensemble-averaged MD simulated XRD patterns (Fig. 1) for synthesized Fh (a3) and MD-relaxed Fh slab (a2) display poorly resolved, broad features suggesting the presence of 2-line nano-crystalline Fh with a small domain size (<10 nm), consistent with PDF data [72, 73]. Based on simulated XRD patterns and radial distribution function (RDF) (Figure S3), MD relaxation of the initial structure of Fh slab for MD simulations based on the Michel model [73] (a1), results in the loss of long-range order to yield a Fh phase (a2) that somewhat resemble the experimental Fh sample (a3) for intense peaks while lacking sharp diffraction peaks.

Fig. 1

Comparative XRD analysis of Fh: structural insights from simulated and synthesized environments with ${\text {[Mo]}_{\text{aq}}}$. XRD pattern for a1 Michel’s crystallographic structure of Fh, a2 MD-relaxed Fh slab, a3 synthesized Fh, b1 Fh slab and water (Mo-0), b2 synthesized Fh and water, c1 Fh slab with ${\text {[Mo]}_{\text{aq}}}$ (Mo-8), and c2 synthesized Fh with ${\text {[Mo]}_{\text{aq}}}$. The plots are offset in y-direction for clarity.

Fh restructures in response to contact with water, which is evident in the XRD patterns of both synthesized (b2) and simulated Fh (b1). These patterns exhibit increased broadening and disorder compared to their respective dry structures. Upon the introduction of ${\text {[Mo]}_{\text{aq}}}$, further adsorption-induced structural changes are observed, characterized by an increased structural disorder in the lower angle region, particularly around $\sim 10^{\circ }\,2\theta$ as seen in (c1). This corresponds to the dominant surface (1$\bar{1}$0) plane of Fh (a2) at $10.50^{\circ }\,2\theta$. The observed changes do not indicate a phase transformation. The increased disorder observed between the initial 0 K model/dry solid and the wet condition, in both the simulated and experimental XRD patterns, is likely due to interactions between dry Fh and water. In simulations, this occurs through long-range interactions, surface charge effects, and atomic exchanges at the surface.

The size of coherent scattering domains [72] determined by PDF signal attenuation of the Fh/water/$\text{MoO}_{4(\text{aq})}$ samples (Figure S5) suggests average crystallite size of 1.5–2 nm. The addition of water and molybdate did not produce a noticeable change in the PDF signal. The observed distances at 1.98 Å, and 3.06 Å align closely with the reported values of 2 Å for the Fe–O bond in tetrahedrally coordinated iron and 3.03 Å for the Fe–Fe bond in FeO₆ polyhedra [72].

The d-PDF (c-b), comprising of peaks around 2 and 3–3.5 Å (Fig. 2A), enhances interfacial signals by isolating the interfacial region from the bulk system. This focus on changes in the PDF of Fh/water introduced by ${\text {[Mo]}_{\text{aq}}}$ allows for a more detailed analysis of interfacial interactions. The simulation of the PDF from the Mo-9 trajectory, containing adsorbed Mo ($[\text {Mo}]_{\text{ads}}$), provides a bottom-up approach for interpreting the measured PDF. This simulated total PDF aligns with the d-PDF at peak positions, exhibiting peaks at around 2 Å and a broad peak with a shoulder at 3–3.5 Å region. Detailed analysis of individual scattering pairs reveals that the simulated PDF comprises peaks at around 2 Å region with contributions from Fe–O and Mo–O pairs with Fe–Fe mainly constituting the feature at around 3–3.5 Å.

Fig. 2

Pair-wise distances at the Fh-water interface: insights from d-PDF and RDF analyses with adsorbed ${\text {[Mo]}_{\text{aq}}}$. A d-PDF plot resulted from the subtraction of Fh/water (b) PDF from Fh/water/${\text {[Mo]}_{\text{aq}}}$ (c) PDF (c–b) along with simulated PDF plots from Mo-9 representing the total calculated PDF (total), and contributions from different scattering pairs Fe–Fe (scaled by 0.5), Fe–O (scaled by 0.5), and Mo–O (scaled by 30). B Calculated RDF plots of Fh in the presence of $[\text {Mo}]_{\text{ads}}$ for FevFe and Fe–O pairs. C Calculated RDF plots of ${\text {[Mo]}_{\text{aq}}}$ with Fh Fe (Fe–Mo) and O (Mo–$\text {O}_{\text{Fh}}$) atoms

The reactive MD approach, successfully reproduces the Fh structure in contact with water and at varying ${\text {[Mo]}_{\text{aq}}}$ concentrations. This is further substantiated by both the measured PDF and RDF, which confirm that the Fh bulk maintains its local structure throughout the simulations aligning with reported values of 1.98 and 3.04 Å for Fe–Fe and Fe–O distances, respectively [72]. Analysis of the RDF from MD trajectories (Fig. 2B), effectively determines pairwise distances during simulations. While actual surfaces have defects that are known to affect Fh local structure and adsorption [7, 74], in case of Fh with varying compositions, crystallinity, and defects, applying ReaxFF to simulate defects would require using larger systems, and a broader range of solids resulting in high the computational costs, making it out of reach for such large system for the moment.

RDF analysis additionally enables the monitoring of prevalent ${\text {[Mo]}_{\text{aq}}}$ distances to the Fh surface (Fig. 2C). Notably, the peak positions remaining unaffected across different ${\text {[Mo]}_{\text{aq}}}$ concentrations suggests the energetically favorable nature of ${\text {[Mo]}_{\text{aq}}}$ interactions with the surface. The first shell peak at 1.73 Å, indicative of inner-sphere, octahedrally-coordinated Mo–O, aligns with values from extended X-ray absorption fine structure (EXAFS) (1.77 Å) [75] and X-ray scattering (1.78 Å) [76]. This distance is comparable with the Mo–O bond length in hematite (1.73 Å) [19], and to As–O (1.68 Å) on Fh [77]. The model architecture allows for exchange processes that extend beyond proton transfer reactions and allows the formation of directly-bound ${\text {[Mo]}_{\text{aq}}}$ complexes which are a result of the exchange of oxygen atoms. RDF analysis also facilitates the observation of adsorbed, fully water-solvated ${\text {[Mo]}_{\text{aq}}}$ species. The second shell peak at 2–3 Å corresponds to outer-sphere adsorbed ${\text {[Mo]}_{\text{aq}}}$ complexes, and water solvation shell (Figure S15). The broad Mo–Fe feature (3.2–4.5 Å) encompasses both inner- and outer-sphere adsorbed ${\text {[Mo]}_{\text{aq}}}$ complexes, aligning with the EXAFS reported range of 3.28–3.53 Å for Fh Mo–Fe [78]. Notably, the proximity of outer-sphere complexes to inner-sphere complexes at these distances suggests the potential for exchange between these states depending on the concentration of ${\text {[Mo]}_{\text{aq}}}$ (vide infra).

Surface interactions and structural dynamics

At equilibrium, the calculated surface charge of the Fh slab reveals a net positively charged top and bottom layer. This positive charge primarily results from Fe atoms and surface protons at $\text {pH} = 7$. This observation aligns with the reported point of zero charge (PZC) for Fh (PZC = 7–9) [79]. The ReaxFF force field used in our simulations implicitly equilibrates charges within the system by minimizing the electrostatic energy with respect to atomic charges [59]. As a result, it cannot model explicit electron-transfer events. To assess the impact of ${\text {[Mo]}_{\text{aq}}}$ adsorption on Fh surface charge, a charge distribution plot was constructed by subtracting the charge of solid Fh from Fh/water/${\text {[Mo]}_{\text{aq}}}$ (Mo-8) system (Fig. 3A). Complementary to Figure S6, the inner-sphere adsorbed complexes of negatively charged ${\text {[Mo]}_{\text{aq}}}$ impacts the Fh surface charge in short-range, localized regions that upon charge equilibration results in the total surface charge of Fh to become more negative.

Fig. 3

Surface charge analysis of Fh: impacts of ${\text {[Mo]}_{\text{aq}}}$ adsorption on surface charge distribution. A Charge probability distribution of Fh in Mo-0 and Mo-8. B Calculated surface charge as a function of the number of ${\text {[Mo]}_{\text{aq}}}$. C Calculated equilibrium atomic charge difference between Mo-0 and Mo-8. Circled areas show the spots with the highest charge perturbations

The charge probability distribution for Fh (Fig. 3B) displays distinct characteristics, with positively-charged Fe and H atoms, alongside negatively-charged oxygen atoms. Notably, the O atom charge distribution exhibits a doublet-like feature. Upon detailed analysis, the more negative peak is attributed to surface oxygens, while the relatively positive one corresponds to bulk Fh oxygen atoms, in water. The ITIM algorithm analysis, both in the presence (Mo-8) and absence (Mo-0) of ${\text {[Mo]}_{\text{aq}}}$, reveals that the more negatively-charged oxygens are located at the interfacial regions on the top or bottom of the Fh slab (Figures S7 and S8). In the case of Mo-8, the adsorption of negatively-charged ${\text {[Mo]}_{\text{aq}}}$ shifts the charge distribution of Fh oxygen atoms towards a more pronounced interfacial character, further indicating a charge transfer to the surface.

Calculating the surface charge as a function of ${\text {[Mo]}_{\text{aq}}}$ concentration (Fig. 3C) indicates that water alone minimally alters the relative surface charge of Fh. As the concentration of ${\text {[Mo]}_{\text{aq}}}$ increases, the surface charge of Fh becomes more negative. However, at higher ${\text {[Mo]}_{\text{aq}}}$ concentrations, the relative surface charge does not continue to decrease but instead reaches a plateau, as observed in the Mo-6 and Mo-8 regimes. This plateau phase is characterized by larger error margins in surface charge, particularly in Mo-8, indicative of non-classical overcharging behavior. This behavior is influenced by pronounced ion-ion association, which affect the structure of the interfacial water hydrogen bond network, competition for adsorption sites on Fh, and the polymerization equilibria of ${\text {[Mo]}_{\text{aq}}}$ species [80].

The structural dynamics of Fh during simulations are further elucidated by analyzing the Fe–O–Fe bond angle probability distribution (Figure S9). The initial nano-crystalline Michel Fh structure [73] comprises both tetrahedral ($109.5^{\underline{o}}$) and octahedral ($90^{\circ }$, $180^{\circ }$) Fe atoms. During MD relaxation and upon contact with water, the sharp features of angle probability distribution broaden, yet the tetrahedral-octahedral coordination environment of Fe is retained. This observation aligns with calculations suggesting that at 298 K, for a simulated Fh with a diameter of $\sim$ 3 nm, the Michel model remains slightly more stable [81]. In our simulations, the initial Michel’s Fh structure does not undergo a phase transformation into a fully octahedral Fe structure, also referred to as the Manceau Fh model. This highlights the stability and persistence of the mixed coordination environment under simulated conditions.

The presence of ${\text {[Mo]}_{\text{aq}}}$, at both low and high concentrations simulated (Mo-3 and Mo-8), does not appear to induce any noticeable changes in the Fe–O–Fe angle distribution (Figure S10). This observation highlights the interfacial nature of ${\text {[Mo]}_{\text{aq}}}$ interactions with Fh, suggesting that the ${\text {[Mo]}_{\text{aq}}}$ primarily affects the surface rather than altering the bulk structure of Fh.

The chemistry of water at ferrihydrite interface and solution

We also explored non-covalent interactions, especially the HB network of water at the Fh-water interface. These interactions influence the structure and dynamics at the solid-water interface [82], impacting geochemical processes such as crystallization [83] and isotopic fractionation [84]. In contrast to classical MD, where hydrogen bonds are not predefined and are instead determined based on energy and geometric criteria [85], the ReaxFF force field declares hydrogen bonding a priori in our reactive simulations. This approach enables more realistic simulations, particularly effective in predicting the short-range cohesive forces within water [86].

The average number of HBs ($\langle N_H \rangle$) (Fig. 4A) gradually increases with distance from the Fh surface. This trend results from the restricted mobility and fewer accessible configurations of water molecules at the interface, leading to a reduced number of HBs compared to bulk water [87]. For the Fh/water system (Mo-0), $\langle N_H \rangle$ approaches that of bulk water at approximately 6 Å from the surface, with a transition region extending from about $\sim$ 2.5 to 6 Å. In the presence of ${\text {[Mo]}_{\text{aq}}}$ (Mo-8), this interfacial-to-bulk HB transition region extends slightly further from the surface, before HB dynamics returns to its unperturbed bulk water state. This effect is most pronounced in the 2–5 Å region from the surface, correlating with the coexistence of inner- and outer-sphere ${\text {[Mo]}_{\text{aq}}}$ complexes, which can act as both HB donors and acceptors.

Fig. 4

Effect of ${\text {[Mo]}_{\text{aq}}}$ on the HB network at the interface of Fh-water: Average number of hydrogen bonds as a function of distance from the surface, A for Fh for Mo-0, and Mo-8. B for surface oxygen $\text {O}_f$ donors and water oxygen $\text {O}_w$ acceptors. Inset) shows surface, water, and ${{\text{MoO}}_{4}}^{2-}$ can act as HB donors and acceptors. For the latter, the simulations show the formation of protonated ${{\text{MoO}}_{4{\text{aq}}}}^{2-}$ species that can act as HB donors (not shown in this image)

${\text {[Mo]}_{\text{aq}}}$ influences the hydrogen HB donor abilities of both surface and water oxygens ($\text {O}_f$ and $\text {O}_w$), with its impact being more pronounced on water oxygens ($\text {O}_w$) in forming HBs with Fh surface oxygens (Figure S11). In the presence of ${\text {[Mo]}_{\text{aq}}}$ (Mo-8), the $\text {O}_w$ donor population shifts away from the surface, compared to the Fh/water system (Mo-0). This shift leads to a decrease in the average number of HBs involving $\text {O}_w$ donors ($\text {O}_w$...$\text {O}_f$), a consequence of the competitive adsorptive interactions between ${\text {[Mo]}_{\text{aq}}}$ and Fh, limiting water’s access to the surface. Regarding Fh surface donors ($\text {O}_f$) (Fig. 4B), the presence of ${\text {[Mo]}_{\text{aq}}}$ in Mo-8 extends the range within which the surface can act as a donor. ${\text {[Mo]}_{\text{aq}}}$ effectively out-competes water as an HB acceptor up to 5 Å from the surface, shaping a new HB network with water.

For the HB acceptor abilities of ${{\text{MoO}}_{4{\text{(aq)}}}}^{2-}$ oxygen atoms ($\text {O}_m$) (Figures S12 and S13), at distances greater than 5 Å, ${\text {[Mo]}_{\text{aq}}}$ is fully solvated and primarily accepts HBs from water. As the concentration of ${\text {[Mo]}_{\text{aq}}}$ increases from Mo-3 to Mo-8, a broader range of HB interactions between $\text {O}_f$...$\text {O}_m$ emerges, indicating the onset of ${\text {[Mo]}_{\text{aq}}}$ adsorption. ${\text {[Mo]}_{\text{aq}}}$ HB donor abilities (Figure S14), suggest that the reactive simulation predicts the formation of protonated ${\text {[Mo]}_{\text{aq}}}$ species with lifetimes long enough to form HBs. Adsorbed ${\text {[Mo]}_{\text{aq}}}$ out-competes water in HB donor abilities up to 5 Å from the surface.

The HB acceptor and donor abilities of $\text {O}_m$ emphasize the crucial role of outer-sphere complexes in mediating surface interactions through HBs [88]. The disruption of HB networks is energetically costly, and the surface HB abilities determine its affinity to water. Additionally, interfacial HB between solutes and water likely reduce its adsorption energy [89]. Consequently, the model predicts various pathways that affect the adsorption free energy.

The reactive simulations provide insights into the behavior of species in bulk solution, particularly highlighting the exchange of oxygen atoms between water and molybdate groups (Figure S15). This exchange is evident when monitoring the distance of Mo with water oxygens within a 4 Å shell. Over time, and as hydration of ${\text {[Mo]}_{\text{aq}}}$ progresses, water oxygens segregate into two distinct groups around ${\text {[Mo]}_{\text{aq}}}$: directly bonded oxygens at $\sim$ 1.78 Å, forming Mo = $\text {O}_{w}$ bonds, consistent with the reported experimental value [76], and the first hydration shell observed at 2–2.5 Å. This is also reflected in RDF of Mo and water oxygen (Mo–Ow) (Figure S16). The results suggest that at equilibrium in solution, Mo–$\text {O}_{m}$ bonds are susceptible to hydrolysis, indicating the dynamic and reactive nature of the hydration process.

Furthermore, the model accurately describes the behavior of solvents and counter-ions. The RDF of water ($\text {O}_w$–$\text {O}_w$) (Figure S17) gives a good reproduction of both experimental and ReaxFF values at room temperature [37]. Notably, there is a slight increase in the height of the first peak of the water $\text {O}_w$–$\text {O}_w$ RDF. This difference may be attributed to stronger interactions between water molecules, potentially heightened by nano-confinement effects between hydrophilic Fh layers [37]. In addition, the RDF of Na–$\text {O}_w$ (Figure S18), suggests a hydrated Na+ cation with a water solvation shell at 2.50 Å. This observation aligns with the reported values for ReaxFF electrolyte systems [55], further validating the model accuracy in describing the solvation structure of Na⁺ in aqueous environments.

Bridging reactive molecular dynamics with surface complexation models

Energetics of adsorption

The electron density profiles of ${\text {[Mo]}_{\text{aq}}}$ as a function of distance from the Fh surface (Fig. 5A) offer a kinetic perspective on ${\text {[Mo]}_{\text{aq}}}$ distribution. This surface is defined as the top interfacial layer of Fh, identified using the ITIM algorithm in each simulation snapshot. The dynamic and rough nature of the Fh surface necessitates this approach.

Fig. 5

Free energy landscape of ${\text {[Mo]}_{\text{aq}}}$ adsorption at the ferrihydrite-water interface. A Electron density profile of ${\text {[Mo]}_{\text{aq}}}$ from the surface of Fh. The surface is defined by the uppermost layer of Fh depicted with the blue region at zero. B Free energy of adsorption of ${\text {[Mo]}_{\text{aq}}}$ as a function of distance from the surface of Fh. Green and red-shaded areas show the regions where inner- and outer-sphere adsorbed complexes exist. The arrows denote the exchange between the two complexes

At infinite dilution (Mo-1), the non-zero electron density at a range of distances of 2.5–10 Å from the surface and farther in the bulk solution, indicates that ${\text {[Mo]}_{\text{aq}}}$ moves between the adsorbed and diffused states. As the concentration of ${\text {[Mo]}_{\text{aq}}}$ increases from Mo-2 to Mo-8, the inner-sphere complexes move closer to the surface at a distance of 2–4 Å, and outer-sphere complexes are positioned in 4–7 Å from the surface. At higher concentrations (Mo-6, Mo-8), more ${\text {[Mo]}_{\text{aq}}}$ are also present at farther distances to up to 10–15 Å and rather than in the bulk solution signifying the strength of adsorptive interactions between ${\text {[Mo]}_{\text{aq}}}$ and Fh. At lower concentrations (Mo-2, Mo-3), the electron density peaks corresponding to inner- and outer-sphere complexes of ${\text {[Mo]}_{\text{aq}}}$ show minimal or no overlap, suggesting less exchange of ${\text {[Mo]}_{\text{aq}}}$ between these two adsorbed states. However, at higher concentrations (Mo-6, Mo-8), there is a notable overlap between the peaks for inner- and outer-sphere complexes. This overlap indicates a facilitated exchange between the two states.

As previously discussed in the Methods section, to comprehensively access the full range of the free energy landscape and the probability distribution of ${\text {[Mo]}_{\text{aq}}}$ as a function of distance from the surface (Fig. 5B) (Figure S21), umbrella sampling is employed. This technique is capable of sampling states that are high in energy, with energy barriers exceeding $k_BT$ (2.479 $\text {kJ mol}^{-1}$ at 298 K), which are often referred to as rare events.

At infinite dilution (Mo-1) the free energy shows a flat profile with small energy barriers between different adsorbed states and bulk states, consistent with the electron density profile. At relatively low concentrations of ${\text {[Mo]}_{\text{aq}}}$ (Mo-2-3), two distinct adsorbed states of—inner- and outer-sphere complexes—are predominantly observed. Inner-sphere complexes, forming at distances of up to 3 Å from the surface, are more energetically favorable by about 15–20 $\text {kJ mol}^{-1}$ compared to outer-sphere complexes at distances > 3–7.5 Å. This energy difference results in a lower desorption rate ($k_{IS \rightarrow OS}$) from the surface consistent with electron density profiles for Mo-3 mostly consisting of inner-sphere complexes and Mo-2 having two distinct absences of electron density at the onset of $\sim$ 2.5 Å. Meanwhile, outer-sphere complexes are located at 4–5 Å from the surface, characterized by a smaller or negligible energy barrier and a higher rate of transition to inner-sphere complexes ($k_{OS \rightarrow IS}$). The observed energy barrier between the inner- and outer-sphere states in both directions is likely attributable to the solvation/desolvation of ${\text {[Mo]}_{\text{aq}}}$, the displacement and reorientation of water in the solvation sphere of ${\text {[Mo]}_{\text{aq}}}$ and at the Fh surface, and the disruption of the hydrogen bond network at the interface [90].

At higher ${\text {[Mo]}_{\text{aq}}}$ concentrations (Mo-6 to Mo-8), the adsorption free energy ranges between − 15 and − 10 $\text {kJ mol}^{-1}$. In this regime, the energy barrier between inner- and outer-sphere complexes disappears entirely. This phenomenon is likely due to the effects of surface charge screening by the adsorbed ${\text {[Mo]}_{\text{aq}}}$ species [91, 92].

Acidity of surface protons

We have selected Fh surface protons based on their equilibrium charge (high, low, medium) (Figure S22) after reactive MD in water (Mo-0) to calculate the $\text {pK}_{\text {a}}$ of the corresponding surface oxygen sites using umbrella sampling. Surface acidity, determined by the $\text {pK}_{\text {a}}$ of surface hydroxyl groups is foundational in predictive SCM. These $\text {pK}_{\text {a}}$’s which are generally found by titration measurements are not capable of assigning individual $\text {pK}_{\text {a}}$’s for different surface protons and are challenging due to mineral dissolution during the titration. The surface reactivity strongly depends on the surface acidity and is critical in understanding the adsorption of solutes at water-mineral interfaces. Computational approaches to calculating unique surface species $\text {pK}_{\text {a}}$ do not rely on fitting experimental parameters and can yield valuable information on surface chemical behavior and its interface with water. The idea behind these computational methods is to calculate the free energy difference between the protonated/deprotonated surface sites mostly through ab initio MD, and PMF techniques [93, 94]. The studied deprotonation sites are selected based on their chemical environment (coordination number-connected to 1 or 2 Fe or species Al–OH vs Si–OH).

The summarized results (Table 1) show surface deprotonation is an endergonic process in water at $\text {pH} = 7$. The results also show that there is no direct relationship between the initial charge of the proton and its $\text {pK}_a$ value giving an emphasis to other factors, such as the local chemical environment of the proton, its interface, (de)hydration effects, and hydrogen bonding affecting the final deprotonated state of the hydrogen atoms. The calculated $\text {pK}_a$ values agree well with reported charge distribution multi-site ion complexation (MUSIC) model calculations for ($1\bar{1}0$) surface of Fh with $\text {pK}_a =7.17$ slightly underestimated while being within 1 $\text {pK}_a$ unit.

Table 1 Fh surface acidity calculations: umbrella sampling results versus MUSIC models

Adsorption constants

The adsorption constants for different ${\text {[Mo]}_{\text{aq}}}$ complexes with the surface of Fh, calculated by SCM, linear free energy relationships (LFER) and the current study are summarized in Table 2. Although the SCM calculations for adsorbed inner-sphere complexes are performed at pH > 5 and with taking into account the acidity of (HO)MoO3− ($\text {pK}_{\text {a}}$ = 4.2), partially or fully protonated Mo species needed to be introduced to improve the model. The calculated adsorption constant for ${\text {MoO}_4}^{2-}$ adsorbed as ${\text {FeOHMoO}_4}^{-2}$, from reactive simulations free energy profile, was found to be 3.14, slightly higher than the value by linear free energy relationships (LFER) calculation. Using the adsorption PMF results (Fig. 5B), the adsorption constants for ${\text {[Mo]}_{\text{aq}}}$ in Mo-1-8 were calculated. The weight-averaged $\text {LogK}_\text {a}$ over all trajectories was then calculated using the Equations S3-S4. The calculated $\text {LogK}_\text {a}$ of 3.39 is attributable to the greater influence (lower error) of data points in regimes of higher concentration, potentially due to a better representation of experimental conditions. The observed slight difference with SCM fitting methodologies can be linked to the use of different parameters in deriving the values. However, the general agreement with the reported values suggests a holistic approach that can be further expanded by varying simulation conditions.

Table 2 Concentration-dependent adsorption constants of ${\text {[Mo]}_{\text{aq}}}$ onto Fh: insights from reactive MD and comparisons with SCM models

Conclusions

Global biogeochemical N, C, and S cycles that have shaped the environmental and biological co-evolution of Earth are dependent upon the mobility and the bioavailability of aqueous nutrient and contaminant elements, such as oxidoreductase enzymes cofactors requiring cellular uptake of aqueous molybdate (${\text {MoO}_4}^{2-}$) ions [45, 46].

Geochemical reactions involving Fe- and S-bearing minerals have regulated ${\text {MoO}_4}^{2-}$ mobility and availability over geologic time [1, 2]. Interfacial reactions with Fe(III) (oxyhydr)oxide ($\text {FeO}_x$) phases often limit aqueous ${\text {MoO}_4}^{2-}$ concentrations in soils, sediments, and waters [4, 5]. Unraveling these complex reaction mechanisms and pathways is important for understanding ${\text {MoO}_4}^{2-}$ availability in modern environments and paleoenvironments. SCM commonly used to describe the adsorption of nutrient and pollutant elements rely on assumptions that can be furthered refined to provide new insights into these interfacial interactions. Here, we employ reactive MD and synchrotron X-ray techniques to show how hydration and adsorption of aqueous ${\text {MoO}_4}^{2-}$ affects the surface charge of Fh which initiates restructuring towards a more disordered phase response observable in simulations and experiments. The Fh surface charge change due to interfacial interactions is concentration-dependent and affects a wide range of chemistry at the interface such as hydrogen bonding and even the dynamics and energetics of exchange between inner- and outer-sphere adsorbed complexes with important implications for the bioavailability and transport of the aqueous species, as well as, implications in isotopic fractionation mechanisms [96]. The accuracy of our model was further tested by deriving SCM parameters for Fh surface acidity and the adsorption constants of ${\text {[Mo]}_{\text{aq}}}$ showing excellent agreement with the literature. This new approach enables informing and improving SCM where experimental work is not feasible or for expanding the chemical conditions explored for accurately describing adsorption via SCM or our current model. Our work is an initial step towards using reactive MD to study adsorption at water-solid interfaces. While it offers valuable insights into adsorption dynamics, it still has limitations that need further investigation [97]. We will continue to rely on established SCM and experimental techniques to improve and validate our modeling approach.