Adsorption Equilibrium and Mechanism and of Water Molecule on the Surfaces of Molybdenite (MoS2) Based on Kinetic Monte-Carlo Method

Abstract

The oxidation/weathering of molybdenite (MoS₂) is too slow to be monitored, even under pure oxygen and high temperatures, while it proceeds rapidly through humid air. The adsorption of water molecules on molybdenite is necessary for the wet oxidation/weathering of molybdenite. Therefore, we employ kinetic Monte Carlo modeling to clarify the adsorption isotherm, site preferences and kinetics of water on different surfaces of molybdenite. Our results indicate that (1) the adsorption capacity and adsorption rate coefficient of H₂O on the (110) surface are significantly larger than those on the (001) surface at a temperature of 0~100 °C and a relative humidity of 0~100%, suggesting that the (110) surface is the predominant surface controlling the reactivity and solubility of molybdenite in its interaction with water; (2) the kinetic Monte Carlo modeling considering the adsorption/desorption rate of H₂O, dissociation/formation rate of H₂O and adsorption/desorption of dissociated H indicates that the adsorption and dissociation of H₂O on the (110) surface can be completed in one microsecond (ms) at 298 K and in wet conditions; (3) the adsorption and dissociation of H₂O on molybdenite are not the rate-limiting steps in the wet oxidation/weathering of molybdenite; and (4) kinetic Monte Carlo modeling explains the experimental SIMS observation that H₂O and OH (rather than H⁺/H⁻ or H₂O) occupy the surface of MoS₂ in a short time. This study provides new molecular-scale insights to aid in our understanding of the oxidation/weathering mechanism of molybdenite as the predominant mineral containing molybdenum in the Earth’s crust.

1. Introduction

Molybdenum (4d⁵5s¹ at the outermost shell) is a typical transition metal and exists mainly as Mo(VI) and Mo(IV) in nature. Therefore, molybdenum acts as an important tracer to construct the paleo-ocean condition due to its redox sensitivity to the aqueous environment [1,2]. Molybdenum is an essential trace element for organisms because it can form a series of important enzymes [3,4]. Its abundance in the Earth’s crust ranges from 0.6 to 1.4 mg/kg [5]. Molybdenite (MoS₂) has been considered as the only mineral category with commercial significance, and it usually coexists with other sulfide minerals (e.g., pyrite, chalcopyrite) in the crust [6,7]. Most of the molybdenum minerals, such as molybdenite (MoS₂), ferrimolybdite (Fe₂(MoO₄)₃·nH2O), molybdite (CaMoO₄), kamilokite (Fe₂Mo₃O₈), wulfenite (PbMoO₄) and drysdallite (Mo(S,Se)₂), are nearly insoluble. Therefore, the weathering/oxidation of molybdenite is the fundamental geochemical process controlling the release of molybdenum from the lithosphere to the hydrosphere.

The oxidation of molybdenite (MoS₂) cannot be monitored by Raman spectroscopy, even when exposed to pure oxygen and temperatures up to 340 °C [8]. However, the wet oxidation/weathering of molybdenite’s surface can be monitored by X-ray photoelectron spectroscopy (XPS) and atomic force microscopy (AFM) after it has been immersed in water for one hour [9] and exposed to humid air [10,11]. All these factors suggest that the adsorption of water molecules on the surface of molybdenite is a precondition for the weathering/oxidation of molybdenite [9]. However, the adsorption isotherm, site preferences and kinetics of water molecules on the different surfaces of molybdenite at the size and time scale of geochemical interest are still poorly understood.

A series of studies applying density functional theory (DFT) were performed to determine the adsorption site stability for H₂O adsorption on the surface of the monolayer of molybdenite [12,13,14], because the adsorption and dissociation of H₂O on the molybdenite surface is difficult to control and discriminate in experiments. According to the adsorption energy of water molecules at different sites, the adsorption site preferences can be compared [12,13,14]. The adsorption of H₂O on the basal plane of the (001) monolayer of MoS₂ shows a positive adsorption enthalpy change, indicating a repulsive interaction between the free H₂O molecule and the perfect MoS₂ (001) surface [13,14]. By comparing the adsorption energy of H₂O above two types of S atoms (S atom with a single S-Mo bond and S atom with double S-Mo bonds) and a Mo atom, the adsorption of H₂O above the S atom was considered unstable due to the repulsive interaction between the S and O atoms, while H₂O adsorption above Mo was considered stable due to the attractive interaction between the O and Mo atoms [14]. Moreover, there are also other possible adsorption geometries to decrease the repulsive interaction between O and S atoms, e.g., the bond between H and S atoms instead of O and S atoms. Therefore, the adsorption site preference of H₂O on the surface of MoS₂ needs more detailed research. Theoretical modeling through first-principle calculation is a useful tool to explore the adsorption mechanism, but present studies are mainly focused on the mono-layer of MoS₂ (001 surface). The surface properties largely depend on the size scale of the mineral surface. In a microscopic view, the interaction among the atoms in the long range can also affect the adsorption site. Therefore, there is always a difference in interest and scale when directly applying physical calculation to solve problems with more geochemical interest. For example, the natural mineral MoS₂(2H) is a three-dimensional object and has three surfaces, i.e., (001), (110), (010). The lack of studies of the water adsorption mechanism (e.g., sorption isotherm, adsorption energy and kinetics) on these surfaces inhibits our fundamental understanding of the oxidation/weathering of molybdenite. To clarify these basic adsorption reactions regarding the interaction at the water–molybdenite surface, we employ annealing dynamics and kinetic Monte Carlo modeling to investigate these fundamental problems.

2. Method

2.1. Thermodynamic Calculation of Eh-pH Diagrams and Reaction Forward Modeling

To fully understand the effect of oxidation on the surface of molybdenite, an Eh-pH diagram is drawn based on the PHREEPLOT [15]. Moreover, forward reaction modeling is performed using the code PHREEQC [16].

2.2. Molybdenite Crystal and Geometry Optimization

Molybdenite mainly exists as 2H-MoS₂ (space group: P63/mmc) in nature. To date, the reported molybdenite crystal categories include (1) the octahedral coordination of Mo (1T) and (2) the trigonal prismatic coordination of Mo (2H, 3R and 4H), which depends on the stacking sequences in S-Mo-S at the z-axis [17]. However, molybdenite in nature exists only as 2H-MoS₂ (space group: P63/mmc) and 3R-MoS₂ (space group: C5-R3m) [18,19,20]. Moreover, 2H-MoS₂ is the predominant crystal category of molybdenite in the deposit [18,19,20]. Therefore, 2H-MoS₂ (Figure 1) was selected for this study. The calculations in this study were performed using the code of Material Studio (Accelrys). To obtain the primary cell of 2H-MoS₂, geometry optimization was performed using the Gaussian-Plane Wave (GPW) method. Different sets of functional and cutoff energy were comparatively used to calculate the lattice constant and band gap. The lattice constant, calculated using the functional PBE, successfully converged, and its result was closer to the experimental value than LDA and PW91 (Table 1). The functional PBE and energy cutoff of 398 eV were selected. The K-point was set as 5 × 5 × 1 in the calculations at the reciprocal space.

2.3. Sorption Isotherm

The dissociation of H₂O into OH and H has been considered to mainly occur at the adsorption site above the Mo atom [13,14,24]. If the adsorption of H₂O above Mo is the most stable adsorption site (with the lowest energy after adsorption), it is very difficult to explain the factors that drive the H₂O above Mo to break the strong O-H bond and dissociate into OH and H. The adsorption of water molecules on the mineral surface is affected by the surface layers in the supercell, while the initial multi-atom system can result in exponential growth in the calculation cost and time. The universal force field (UFF) was used in the molecular dynamics modeling of H₂O adsorbed on the (001) and (110) surfaces of molybdenite. UFF is constructed using the general rules only based on the element, its hybridization and its connectivity [25], and it has been widely used in calculating the interactions among atoms in molecular dynamics simulation. The 5 × 5 × 2 supercell of each surface, (001) and (110), was constructed in the calculation in order to more closely reflect the real size of molybdenite and consider the long-range interactions among atoms. At a constant temperature, average numbers of H₂O adsorbed for each cell of the (001) and (110) surfaces of molybdenite were sampled and calculated from the atmospheric pressure from 1 bar to 10 bar. The temperature was increased from 273 to 398 K (0~120 °C) with a step of 25 K. The summation method was based on the Ewald Group and atom-based van der Waals. At each simulation step, the metropolis method was used to sample the ensembles. The probability density of the canonical ensemble is shown in Equation (1).

(1)$\mathsf{\rho }\left(\mathsf{\Gamma }\right)=\frac{\exp \left(-\mathsf{\beta }\mathrm{E}\left(\mathsf{\Gamma }\right)\right)}{\int \mathrm{d}\mathsf{\Gamma }\prime \exp \left(-\mathsf{\beta }\mathrm{E}\right))}$

where E(Γ) is the potential energy (eV) of the system in state Γ and β = 1/(K_bT), in which K_b is the Boltzmann constant and T is the thermodynamic temperature (K).

2.4. Kinetic Monte Carlo and Rate Coefficient Calculation

The most likely site and geometry for H₂O adsorbed on the (001) and (110) surfaces of MoS₂ were explored through the annealing process, which was realized by increasing the temperature of the system in steps to 100,000 K, and then automatically decreasing it to 100 K. In the molecular dynamics process, we also used the UFF, but with the Monte Carlo method to sample the energy of each state. The process was repeated in 3 cycles with 15,000 steps in each cycle to explore the state with lower energy.

Based on the Monte Carlo kinetics, rate activation energy and rate coefficient calculated, we investigated the dissociative adsorption rate of water at the mineral surface. Kinetic Monte Carlo (KMC) modeling is based on the rate coefficient or activation energy of each step reaction to simulate the concentration and rate change as a function of time, being comparable to laboratory settings. Its visualization result can help us to understand the reaction process at the atomic scale. According to the experimental observation of OH and H₂O on the mineral surface [24] and the theoretical prediction of H desorption as H₂ from the surface, in this study, KMC is employed to combine the three-stage reactions, which requires the rate coefficient or activation energy of H₂O adsorption, H₂O dissociation and H₂ desorption. We did not consider the diffusion of these species at the surface, because the effect of diffusion on the overall reaction rate is not significant if the time scale concerned is in seconds or less. The adsorption rate coefficient of H₂O adsorption on the (001) surface was calculated based on Equation (2) [26]:

(2)${\mathrm{K}}_{{\mathrm{H}}_2\mathrm{O}-\mathrm{ads}}=\frac{{\mathrm{P}}_{{\mathrm{H}}_2\mathrm{O} }{\mathrm{A}}_{\mathrm{site}}б}{\sqrt{2{\mathsf{\pi }\mathrm{m}}_{{\mathrm{H}}_2\mathrm{O}}{\mathrm{K}}_{\mathrm{b}}\mathrm{T}}}$

where K_H2O-ads is the rate coefficient for H₂O adsorption at the surface of the mineral (s⁻¹). A_site is the area of a single adsorption site (m²), which can be estimated based on the geometry of the cleaved surface in the calculation. PH₂O is the pressure of water (atm). б is the sticking coefficient (s⁻¹). mH₂O is the molar mass of water molecules (g/mol). K_b is the Boltzmann constant (J/K) and T is the thermodynamic temperature (K).

The water vapor saturation pressure as a function of temperature is described with the Arden–Buck equation; see Equation (3) [27]. Its error compared to the experimental value [28] is less than 0.04% at the temperature range of 273.15 to 373.15 Kelvin.

(3)$\mathrm{P}=0.61121 \exp \left[\left(18.678-\frac{\left(T-273.15\right)}{234.5}\right)\ast \left(\frac{\left(T-273.15\right)}{257.14+\mathrm{T}}\right)\right]$

P is the water vapor saturation pressure (KPa) and T is the thermodynamic temperature (K). Therefore, we can obtain the water adsorption rate coefficient reaching water vapor saturation at a certain temperature as Equation (4):

(4)${\mathrm{K}}_{{\mathrm{H}}_2\mathrm{O}-\mathrm{ads}}=\frac{0.00603 \exp \left[\left(18.678-\frac{\left(T-273.15\right)}{234.5}\right)\ast \left(\frac{\left(T-273.15\right)}{257.14+\mathrm{T}}\right)\right]{\mathrm{A}}_{\mathrm{site}}б}{\sqrt{2{\mathsf{\pi }\mathrm{m}}_{{\mathrm{H}}_2\mathrm{O}}{\mathrm{K}}_{\mathrm{b}}\mathrm{T})}}$

where K_b is the Boltzmann constant (J/K) and T is the thermodynamic temperature (K).

The rate coefficient of H₂O dissociation at the (001) MoS₂ surface has been reported previously [14]. The rate coefficient of H₂ desorption from the (001) MoS₂ surface has not been reported and no direct equation is available. The H₂ desorption from the (001) MoS₂ surface can be expressed as in Reaction 2. When obtaining the equilibrium between H₂ adsorption and desorption on the (001) surface, the product of the H₂ adsorption rate coefficient and H₂ desorption rate coefficient will be equal to its equilibrium constant. Therefore, Equation (5) can be used to calculate the rate coefficient of H₂ desorption from the (001) MoS₂ surface (the detailed process is summarized in Supplementary Materials Section S1 as an elementary reaction.

2H-MoS₂ = MoS₂ + H₂ (Reaction 2)

(5)${\mathrm{K}}_{\mathrm{eq}}=\frac{{\mathrm{K}}_{\mathrm{forward}}}{{\mathrm{K}}_{\mathrm{reverse}}}=\frac{{\mathrm{K}}_{{\mathrm{H}}_2-\mathrm{desorption}}}{{\mathrm{K}}_{{\mathrm{H}}_2-\mathrm{adsorption}}}$

where K_eq is the equilibrium constant of H₂ adsorption and desorption at the (001) surface; K_forward and K_reverse are the rate coefficients (s⁻¹) of the forward and reverse reaction, respectively; K_{H2-desorption} and K_{H2-adsorption} are the rate coefficients of H₂ desorption and adsorption, respectively. The detailed calculation is presented in the Supplementary Materials Section S1. Supplementary Materials Section S1. The rate coefficient of H₂ adsorption was also calculated using the same equation, Equation (2). Based on these parameters, the overall and step reaction rate in H₂O’s dissociative adsorption as a function of time was predicted (Figure 2).

3. Results

3.1. Oxidation of Molybdenite Based on Thermodynamic Calculation

Based on the Eh-pH diagram of the Mo-S-H₂O system at 25 °C and reference atmospheric pressure, we can find that MoS₂ is stable under the acid and reductive conditions (Figure 2) and can be oxidized and transformed into Mo (VI) species in the surface environment. The oxidation energy difference varies from −0.9 to −2.4 eV based on the DFT calculation [29]. In acidic water (pH < 4.3) in the relatively oxidative condition, MoS₂ can also remain stable (Figure 3). The forward path modeling of MoS₂ oxidized by O₂ indicates that oxidation will improve its solubility unlimitedly if there is sufficient O₂/air and will decrease the pH (Figure 3). The increase in Mo (VI) and decrease in pH will form a series of aqueous species and complexes (Figure 3), such as MoO₄²⁻, HMoO₄⁻, H₂MoO₄, Mo₇O₂₄⁶⁻. All of these suggest that the oxidation of molybdenite is very important in understanding the Mo aqueous species during oxidation.

2MoS₂ + 6H₂O + 9O₂ = 2MoO₄²⁻ + 12H⁺ + 4SO₄²⁻ (Reaction 1)

3.2. Sorption Isotherm of H₂O on the (001) and (110) Surfaces

All of the movement of H₂O at the surface of MoS₂ is random in the simulation. The sorption result of H₂O on the (001) and (110) surfaces of MoS₂ from 1 bar to 10 bar at 298 K is presented (Figure 4). The sorption capacity comparison and the isotherms at different temperatures and pressures are presented (Figure 5 and Figure 6). According to the number of water molecules adsorbed on the supercells of the (001) and (110) surfaces from 1 to 10 bar (Figure 4), we can calculate the average H₂O molecules adsorbed for each cell (Figure 5). H₂O molecules adsorbed on the (001) surface are clearly distributed in different layers in the space, while those water molecules adsorbed on the (110) surface are relatively evenly distributed in the space (Figure 4). The average loadings of H₂O molecules for the cells of the surfaces (001) and (110) at 298 K both increased with the H₂O fugacity in the range of 1 to 10 bar (Figure 5). Moreover, the H₂O molecule sorption capacity of the (110) surface at 298 K was much larger than that of the d(001) surface, and the capacity difference between the two surfaces also increased significantly with the H₂O fugacity (Figure 5). Based on the average loading of H₂O molecules sorbed on the cell and the stoichiometry mass of H₂O and the cell, we can obtain the sorption isotherm from the (001) and (110) surfaces as a function of H₂O fugacity and temperature (Figure 6). There were significant and highly positive linear relationships between the H₂O capacity on the two surfaces and the H₂O fugacity (Figure 6). When the temperature increased from 273 to 398 K, the H₂O sorption capacity on both the (001) and (110) surfaces decreased significantly, as shown by the decrease in slope for each isotherm line at the two surfaces (Figure 6). The sorption capacity of H₂O on the (110) surface was stronger than that of the (001) surface at all temperatures and H₂O fugacity values investigated (Figure 6).

3.3. Adsorption Energy Derived from the Annealing Dynamics of H₂O on the (001) and (110) Surfaces

The annealing dynamics from 100,000 K to 100 K indicated that there were two possible adsorption sites on the (001) surface, 001-G1 and 001-G2, whose adsorption energy was −2.51 and −2.23 eV, respectively (Figure 7), and there were four possible adsorption sites on the (110) surface—110-G1, 110-G2, 110-G3 and 110-G4—whose adsorption energy was −2.56, −2.36, −2.15, −1.93 and −1.48 eV, respectively (Figure 7). The stability of different adsorption sites on the surfaces can be reflected by the potential energy. Correspondingly, the stability of the adsorption sites on molybdenite follows the order 110-G1 > 110-G2 > 001-G1 > 001-G2 > 110-G3 > 110-G4 > 110-G5, as shown by the increase in negative adsorption energy for H₂O in the order (Figure 7). Moreover, the four most stable geometries of H₂O adsorption sites are presented (Figure 8 and Figure 9). To clarify the geometry on the sites, the most stable adsorption sites on the (001) surface were (1) 001-G1, H₂O above the S atom at the edge; (2) 001-G2, H₂O above the S atom closing the edge. Both 001-G1 and 001-G2 were above the S atoms, with two H atoms oriented towards the S atom. The two most stable adsorption sites on the (110) surface were (1) 110-G1, H₂O above the S-Mo bond at the edge, with two H atoms oriented towards the S atom and an O atom oriented towards the Mo atom (Figure 8); (2) 110-G2, H₂O above the S-Mo bond close to the edge, with two H atoms oriented towards the S atom and an O atom oriented towards the Mo atom (Figure 9). The distance between two H atoms and an S atom for 001-G1, 001-G2, 110-G1 and 110-G2 was 3.311–3.313 3.128–3.198, 3.382–3.392, 3.326–3.355 Å, respectively, and only 001-G1 had two equal S-H distances (Figure 8 and Figure 9). The specific distances between O and H atoms from H₂O and Mo and H atoms from MoS₂ for the four most stable sites are presented (Figure 8 and Figure 9). The distance between the O atom and Mo atom of 001-G1 and 001-G2 could not be determined due to the uncertain direction and overly long distance (Figure 8).

3.4. Adsorption Rate Coefficient of H₂O on the (001) and (110) Surfaces as a Function of the Temperature and Humidity

Based on Equation (2), we can calculate the adsorption rate coefficient of H₂O on the mineral surface at certain partial pressure and temperature values; the A_site is the area of a single sorption site, which equals the reciprocal of the sorption site density. According to the area of the basal plane and the H₂O adsorption number per cell, we can calculate the area of a single sorption site at the (001) and (110) surfaces, respectively (Table 2). Combined with the bulk equation, we can calculate the saturation vapor pressure of H₂O at certain temperature ranges from 0 to 100 °C (Figure 10a) and calculate the corresponding adsorption rate coefficient at certain conditions. The adsorption rate coefficient of H₂O on the (110) surface is close to that on the (001) surface at a temperature from 0 to 25 °C, but it increases more significantly than the (001) surface with the temperature. When reaching 100 °C, the adsorption rate coefficient of H₂O on the (110) surface is approximately 16 times higher than that on the (001) surface. To explore the impact of humidity on the adsorption of H₂O on MoS₂, the adsorption rate coefficient of H₂O on the (001) and (110) surfaces of MoS₂ at 25 °C was presented as a function of relative humidity (Figure 10b). The adsorption rate coefficients of H₂O on the (110) and (001) surfaces both increased with the humidity.

3.5. Kinetics of H₂O Dissociative Adsorption on the (110) Surface

The adsorption mechanism and rate coefficient both indicate that H₂O adsorption on the (110) surface is always favored at all of the temperature and humidity levels investigated. Therefore, the surface reactions of H₂O on the (110) surface are selected in the KMC modeling, which includes the kinetics of (1) the adsorption of H₂O on the surface; (2) the dissociation of H₂O into OH and H on the surface; (3) the desorption of H as H₂ molecules from the surface. According to Equations (2) and (3), we can obtain the rate coefficient of the three-stage reaction. When reaching the condition with 20% relative humidity of H₂O in air at 298 K, the adsorption rate coefficient is 4.84 × 10⁷ S⁻¹. The rate coefficient of H₂O dissociation on the (110) surface of MoS₂ at 300 K is 2.8 × 10¹⁰ S⁻¹. The rate coefficient of H₂ desorption on MoS₂ is 2.38 × 10⁸ S⁻¹ (detailed calculation is presented in Supplementary Materials Section S1.). To more closely approach the real situation, the rate coefficient of H₂O desorption from MoS₂ was also calculated (detailed calculation is presented in Supplementary Materials Section S3). The metadynamic calculation provides the rate coefficient of H₂O dissociation on the MoS₂ surface [14]. Based on the reaction free energy of the dissociation, we can calculate the reverse reaction (the formation of H₂O from dissociated OH and H) rate coefficient as 2.9 × 10¹¹ S⁻¹ (detailed calculation is presented in Supplementary Materials Section S2). The average H₂O number adsorbed on the (110) surface is approximately 0.85 per cell at 298 K at the water vapor saturation pressure. Therefore, we constructed a surface composed of 32 × 32 adsorption sites, which represents a (110) surface composed of 37 × 37 cells, to predict the overall surface reactions as a function of time in seconds.

Reaction scheme in kinetic Monte Carlo and rate coefficient are presented in Table 3. Based on the reaction rate coefficient of these reactions, the kinetic model of the surface reaction on the (110) surface shows that H₂O is firstly adsorbed on the surface and 20% of the adsorption sites are occupied by H₂O, and only a few OH and H atoms appear at the surface at 10–20 ns (Figure 11 and Figure 12). The concentration of H₂O adsorbed on the (110) surface increases at a reaction time of less than 40 ns, while it gradually decreases afterwards. The concentration of OH on the surface increases rapidly with time and reaches approximately 90% at a reaction time of 1 ms (Figure 11 and Figure 12). The concentration of H adsorbed on the surface is very low and decreases to nearly zero in 1 ms (Figure 11 and Figure 12).

4. Discussion

4.1. Adsorption Capacity and Rate Comparison of H₂O on the Different Surfaces of Molybdenite and Its Implications for the Mineral Reactivity

The distribution of H₂O molecules in the space above the surface is different between surfaces (001) and (110), indicating that the sorption capacity and mechanism of H₂O on the two surfaces are significantly different (Figure 4). The H₂O adsorption on the (001) surface is multi-layer adsorption, while the H₂O adsorption on the (110) surface is single-layer adsorption (Figure 4). The periodic adsorption of H₂O above the (001) surface (Figure 4) suggests that it is probably the result of the periodic repulsive force between the O atom of H₂O and the S atom of MoS₂. The sorption capacity of H₂O on the two surfaces increases with the H₂O fugacity but decreases with the temperature (0–120 °C), indicating that the relatively wet and cold environment is beneficial for the H₂O molecules’ adsorption capacity on the mineral surface (Figure 7), while a very low temperature can decrease the adsorption rate of H₂O on the surface (Figure 10). The sorption isotherm of H₂O shows that the sorption capacity of the (110) surface is much stronger than that of the (001) surface (Figure 5 and Figure 6), indicating that the (110) surface can obtain more water adsorbed on its surface in the weathering of MoS₂ and its interaction with water.

Although a temperature decrease can improve the sorption capacity of H₂O on the MoS₂ surface, it can significantly decrease the adsorption rate of H₂O on the MoS₂ surface, suggesting that the effect of temperature on water sorption on MoS₂ is complicated. The adsorption rate coefficient of H₂O on the (110) surface was always larger than that on the (001) surface at all the temperatures (0~100 °C) and relative humidity (0~100%) investigated (Figure 10). Therefore, both the adsorption capacity and rate coefficient of H₂O on the (110) surface are stronger than those on the (001) surface, suggesting that the (110) surface is the predominant surface that controls the solubility and reactivity of MoS₂ in its interaction with water, because the (110) surface can attract many more water molecules, and the atoms of Mo and S and their bonds are also exposed to the H₂O molecules.

4.2. Adsorption Mechanism of H₂O on the Surfaces of Molybdenite

Compared to the possible adsorption geometry optimization using DFT techniques at 0 K, the annealing dynamics method explores many more possible adsorption geometries at a larger temperature range from 10 to 100,000 K. The most stable adsorption sites of H₂O on the (001) and (110) surfaces derived from annealing dynamics are much more stable than the sites derived from optimization with DFT. For example, the most stable adsorption site on the Mo and S atoms at the edge of the (001) surface (equivalent to the (110) surface) has the lowest adsorption energy of −0.55 eV [14], which is much higher than that of the adsorption sites with adsorption energy lower than −2.5 eV, e.g., 001-G1 and 110-G1 (Figure 7). The adsorption sites of H₂O on the (001) surface have been considered unstable (positive adsorption energy) due to the repulsive interaction between S and O atoms [14], but the annealing dynamics indicated that H2O molecules adsorbed on the (001) surface can also remain stable (e.g., 001-G1 and 001-G2 with negative adsorption energy lower than −2.2 eV) by means of changing the direction of H₂O molecules to direct the two H atoms towards the S atom at the surface in order to lower the total energy and repulsive interaction (Figure 8). For the (110) surface, the adsorption site 001-G1 has larger distances between Mo and O atoms and between S and H atoms than 001-G2, which causes the O atom (of H₂O) and S atom (of MoS₂) to receive more electrons from the Mo atom (of MoS₂) and H atom (of H₂O). However, the 001-G1 adsorption site is more stable than 001-G2, as shown by the lower adsorption energy (Figure 7), because 001-G1 has larger distances between O and S and between H and Mo atoms (Figure 9b,d), suggesting that the repulsive interactions between O and S and between H and Mo atoms are the major factors that control the energy level of the adsorption geometry. These results suggested that the annealing dynamics and larger supercells with more periodic clusters can provide more possible and stable sites than the common cluster optimization, which provides new, complete insights in understanding the adsorption mechanism at the level of geochemical interest.

4.3. Kinetic Modeling of Dissociative Adsorption of H₂O on the (110) Surface and Its Implication for the Weathering of Molybdenite in the Surface Environment

We selected the (110) surface to predict the progress of H₂O adsorption and dissociation as a function of time due to the great reactivity of the (110) surface of molybdenite. Although a previous study showed that the H₂O adsorbed on molybdenite can dissociate into OH and H at room temperature [15], the effects of the adsorption/desorption rate of water, the formation rate of OH and H into H₂O and the desorption and sorption rate of H on the overall reaction rate have not been considered. To more closely reflect the real situation, we calculated and considered the reaction rates of all of these step reactions in the KMC modeling, which was firstly applied to predict the dissociation adsorption of the H₂O reaction on the (110) surface as a function of time (Figure 11 and Figure 12). The modeling result indicated that the adsorption and dissociation of H₂O into OH and H can occur in the time scale of ns (10⁻⁹ s) at 25 °C and in a humid environment (Figure 11 and Figure 12). Moreover, the OH can occupy most adsorption sites of the (110) surface in one microsecond (Figure 11 and Figure 12), suggesting that OH covers the surface of MoS₂ in a very short time, and the dissociative adsorption of H₂O is not the rate-limiting step in the oxidation of MoS₂ at a normal temperature. This explains the experimental SIM observation that H₂O and OH (rather than H⁺/H⁻ or H₂) occupied the surface of MoS₂ in a short time [24]. In the overall reaction of MoS₂ oxidation (Reaction 1), the adsorption/capture of molecular oxygen on the OH shell around molybdenite may be an important reaction controlling the rate of MoS₂ oxidation/weathering, which needs further detailed investigation.

5. Conclusions

By employing molecular dynamics and kinetic Monte Carlo modeling, we investigated the adsorption isotherm of H₂O on the (001) and (110) surfaces of molybdenite at temperatures ranging from 0 to 120 °C and H₂O fugacity ranging from 1 to 10 bar, which indicated that an increase in fugacity and decrease in temperature can promote adsorption, but a very low temperature can lower the adsorption rate coefficient. The adsorption geometry of H₂O on the (001) and (110) surfaces at 298 K was explored using annealing dynamics, which provides a more stable geometry with lower adsorption energy, relative to a previous report using DFT geometry optimization. The adsorption rate coefficient and capacity of H₂O on the two surfaces indicated that the (110) and (010) surfaces are the predominant surfaces that control the solubility and reactivity of MoS₂ in its interaction with water. The reaction time of 1 ms (10⁻³ s) is sufficient for H₂O to adsorb and dissociate into OH at the (001) surface, suggesting that the adsorption and dissociation of H₂O on the MoS₂ surface is not the rate-limiting step in MoS₂ oxidation/weathering.

Author Contributions

Writing: R.W., Writing and visualization: X.W., Concept and method: Z.Z., Funding acquisition: S.N., J.D. and D.W. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All the data used is included in the article and Supplementary Materials.

Acknowledgments

Thanks should be given to the anonymous reviewers and academic editors.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

KMC	Kinetic Monte Carlo
XPS	X-ray photoelectron spectroscopy
AFM	Atomic force microscopy
DFT	Density functional theory
SIMS	Secondary ion mass spectroscopy
GPW	Gaussian-Plane Wave
PBE	Perdew–Burke–Eenzerhof
LDA	Local-density approximation
PW91	Perdew-Wang, 1991
UFF	Universal force field

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Supercell of MoS2(2H) crystal (5 × 5 × 2) (a): three-dimensional, (b): top view (001), (c): side view (110) (yellow ball is sulfur atom and blue ball is molybdenum atom).

Figure 2. Eh-pH diagram of Mo-S-H2O system at 25 °C and 1 bar.

Figure 3. Oxidation process of molybdenite and speciation change as a function of the number of reaction steps.

Figure 4. Sorption saturation status of H2O molecule on the supercell (5 × 5 × 2) of (001) surface at (a) 1 bar and (b) 10 bar, and (110) surface at (c) 1 bar and (d) 10 bar (one red point represents one water molecule; yellow ball represents sulfur atom; blue ball represents molybdenum atom).

Figure 5. Sorption capacity change of H2O on the (001) and (110) surfaces of MoS2 at 1 and 10 bar.

Figure 6. Sorption isotherm of H2O molecules on the (a): (001) and (b): (110) surfaces as a function of temperature.

Figure 7. Adsorption energy of H2O at different adsorption sites on the (001) and (110) surfaces.

Figure 8. Geometry of the two most stable adsorption sites for H2O on the (001) surface: (a,b): 001-G1; (c,d): 001-G2 (yellow ball represents sulfur atom; blue ball represents molybdenum atom).

Figure 9. Geometry of the two most stable adsorption sites for H2O at the (110) surface: (a,b): 110-G1; (c,d): 110-G2 (yellow ball represents sulfur atom; blue ball represents molybdenum atom).

Figure 10. Adsorption rate coefficient of water molecules on molybdenite as a function of temperature (a) and relative humidity (b).

Figure 11. KMC modeling of the dissociative adsorption of H2O on the molybdenite (001) surface in one microsecond at 298 K.

Figure 12. Relative concentration variance of H2O, OH, H on the molybdenite (110) surface in one microsecond at 298 K.

Table 1

Comparison of DFT calculation using different functional and experimental parameters in MoS2 cell geometry optimization.

Source	Lattice Constant—a (Å)	Lattice Constant—c (Å)	Mo-S Bond Length (Å)	S-S Bond Length (Å)	Band Gap (eV)
LDA	3.136	12.052	2.381	3.110	0.74
PBE	3.168	12.616	2.407	3.128	1.102
PW91	3.180	12.690	2.411	3.126	1.023
Dickinson and Pauling, 1923 [21]	3.15	12.3			1.29
Swanson et al., 1955	3.16	12.295			1.29
Bronsema et al., 1986 [22]	3.16	12.296
Zubavichus et al., 1998 [23]			2.39	3.16

Table 2

Adsorption rate coefficient of H₂O (K_H2O-ads) on the (001) and (110) surfaces of molybdenite at temperature (t) from 0 to 100 °C and corresponding water vapor saturation pressure (P_vp).

	T (°C)	Pvp (atm)	A-Site (Å2)	KH2O-ads (S−1)
(001)	0	0.006	3.23	2.97 × 107
	25	0.0313	4.63	2.13 × 108
	50	0.1218	4.57	7.84 × 108
	75	0.3806	6.27	3.24 × 109
	100	1	6.09	7.98 × 109
(110)	0	0.006	5.68	5.22 × 107
	25	0.0313	7.31	3.35 × 108
	50	0.1218	8.92	1.53 × 109
	75	0.3806	12.10	6.25 × 109
	100	1	11.61	1.52 × 1010

Table 3

Reaction scheme in kinetic Monte Carlo dynamics and rate coefficient.

Reaction	Rate Coefficient (S−1)
MoS2 + H2O = MoS2-H2O	4.84 ×107
MoS2-H2O = MoS2 + H2O	3.62 × 1010
MoS2-H2O = MoS2-OH + H	2.8 × 1010
MoS2-OH + H = MoS2-H2O	2.9 × 1011
2H-MoS2 = MoS2 + H2	2.38 × 108
MoS2 + H2 = 2H-MoS2	1.20 × 108

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/molecules27248710/s1. Section S1: Desorption rate of H₂ from MoS₂ unit cell based on geometry optimization; Section S2: Formation and dissociation rate coefficient of H₂O above MoS₂; Section S3: Adsorption and desorption rate coefficient of H₂O above MoS₂.

References

1. Gardner, C.B.; Carey, A.E.; Lyons, W.B.; Goldsmith, S.T.; McAdams, B.; Trierweiler, A.M. Molybdenum, vanadium, and uranium weathering in small mountainous rivers and rivers draining high-standing islands. Geochim. Cosmochim. Acta; 2017; 219, pp. 22-43. [DOI: https://dx.doi.org/10.1016/j.gca.2017.09.012]

2. Greaney, A.T.; Rudnick, R.L.; Gaschnig, R.M.; Whalen, J.B.; Luais, B.; Clemens, J.D. Geochemistry of molybdenum in the continental crust. Geochim. Cosmochim. Acta; 2018; 238, pp. 36-54. [DOI: https://dx.doi.org/10.1016/j.gca.2018.06.039]

3. Mendel, R.R. Biology of the molybdenum cofactor. J. Exp. Bot.; 2007; 58, pp. 2289-2296. [DOI: https://dx.doi.org/10.1093/jxb/erm024] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17351249]

4. Schwarz, G.; Mendel, R.R.; Ribbe, M.W. Molybdenum cofactors, enzymes and pathways. Nature; 2009; 460, 839. [DOI: https://dx.doi.org/10.1038/nature08302] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19675644]

5. Wedepohl, K.H. The composition of the continental crust. Geochim. Cosmochim. Acta; 1995; 59, pp. 1217-1232. [DOI: https://dx.doi.org/10.1016/0016-7037(95)00038-2]

6. Hess, F.L. Molybdenum Deposits, A Short Review; Government Publishing Office: Washington, DC, USA, 1924.

7. LeAnderson, P.J.; Schrader, E.L.; Brake, S.; Kaback, D.S. Behavior of molybdenum during weathering of the Ceresco Ridge porphyry molybdenite deposit, Climax, Colorado and a comparison with the Hollister deposit, North Carolina. Appl. Geochem.; 1987; 2, pp. 399-415. [DOI: https://dx.doi.org/10.1016/0883-2927(87)90024-2]

8. Yamamoto, M.; Einstein, T.L.; Fuhrer, M.S.; Cullen, W.G. Anisotropic etching of atomically thin MoS2. J. Phys. Chem. C; 2013; 117, pp. 25643-25649. [DOI: https://dx.doi.org/10.1021/jp410893e]

9. Zhang, X.; Jia, F.; Yang, B.; Song, S. Oxidation of molybdenum disulfide sheet in water under in situ atomic force microscopy observation. J. Phys. Chem. C; 2017; 121, pp. 9938-9943. [DOI: https://dx.doi.org/10.1021/acs.jpcc.7b01863]

10. Ross, S.; Sussman, A. Surface oxidation of molybdenum disulfide. J. Phys. Chem.; 1955; 59, pp. 889-892. [DOI: https://dx.doi.org/10.1021/j150531a020]

11. Buck, V. Preparation and properties of different types of sputtered MoS2 films. Wear; 1987; 114, pp. 263-274. [DOI: https://dx.doi.org/10.1016/0043-1648(87)90116-5]

12. Raybaud, P.; Hafner, J.; Kresse, G.; Kasztelan, S.; Toulhoat, H. Ab initio study of the H2–H2S/MoS2 gas–solid interface: The nature of the catalytically active sites. J. Catal.; 2000; 189, pp. 129-146. [DOI: https://dx.doi.org/10.1006/jcat.1999.2698]

13. Ataca, C.; Ciraci, S. Dissociation of H₂O at the vacancies of single-layer MoS₂. Phys. Rev. B; 2012; 85, 195410. [DOI: https://dx.doi.org/10.1103/PhysRevB.85.195410]

14. Ghuman, K.K.; Yadav, S.; Singh, C.V. Adsorption and dissociation of H2O on monolayered MoS2 edges: Energetics and mechanism from ab initio simulations. J. Phys. Chem. C; 2015; 119, pp. 6518-6529. [DOI: https://dx.doi.org/10.1021/jp510899m]

15. Kinniburgh, D.G.; Cooper, D.M. PhreePlot: Creating Graphical Output with PHREEQC; Centre for Ecology and Hydrology: Gwynedd, UK, 2011.

16. Parkhurst, D.L.; Appelo, C.A.J. Description of Input and Examples for PHREEQC Version 3—A Computer Program for Speciation, Batch-Reaction, One-Dimensional Transport, and Inverse Geochemical Calculations; U.S. Geological Survey Techniques and Methods: Denver, CO, USA, 2013.

17. Raybaud, P.; Hafner, J.; Kresse, G.; Toulhoat, H. Ab initio density functional studies of transition-metal sulphides: II. Electronic structure. J. Phys. Condens. Matter; 1997; 9, 11107. [DOI: https://dx.doi.org/10.1088/0953-8984/9/50/014]

18. Stotu, U.; Marnntall, T. Molybdenite polytypes in theory and occurrence. II. Some naturally-occurring polytypes of molybdenite. Am. Mineral.; 1970; 55, pp. 1857-1875.

19. Luck, J.M.; Allègre, C.J. The study of molybdenites through the ¹⁸⁷Re-¹⁸⁷Os chronometer. Earth Planet. Sci. Lett.; 1982; 61, pp. 291-296. [DOI: https://dx.doi.org/10.1016/0012-821X(82)90060-7]

20. Peng, J.; Zhou, M.F.; Hu, R.; Shen, N.; Yuan, S.; Bi, X.; Du, A.; Qu, W. Precise molybdenite Re–Os and mica Ar–Ar dating of the Mesozoic Yaogangxian tungsten deposit, central Nanling district, South China. Miner. Depos.; 2006; 41, pp. 661-669. [DOI: https://dx.doi.org/10.1007/s00126-006-0084-4]

21. Dickinson, R.G.; Pauling, L. The crystal structure of molybdenite. J. Am. Chem. Soc.; 1923; 45, pp. 1466-1471. [DOI: https://dx.doi.org/10.1021/ja01659a020]

22. Bronsema, K.D.; De Boer, J.L.; Jellinek, F. On the structure of molybdenum diselenide and disulfide. Z. Für Anorg. Und Allg. Chem.; 1986; 540, pp. 15-17. [DOI: https://dx.doi.org/10.1002/zaac.19865400904]

23. Zubavichus, Y.V.; Slovokhotov, Y.L.; Schilling, P.J.; Tittsworth, R.C.; Golub, A.S.; Protzenko, G.A.; Novikov, Y.N. X-ray absorption fine structure study of the atomic and electronic structure of molybdenum disulfide intercalation compounds with transition metals. Inorg. Chim. Acta; 1998; 280, pp. 211-218. [DOI: https://dx.doi.org/10.1016/S0020-1693(98)00209-6]

24. Karolewski, M.A.; Cavell, R.G. SIMS study of Cs/MoS₂ (0001): II. Chemisorption of O₂, H₂O, HCOOH, CO₂ and CS₂. Surf. Sci.; 1989; 219, pp. 261-276. [DOI: https://dx.doi.org/10.1016/0039-6028(89)90212-4]

25. Rappe, A.K.; Casewit, C.J.; Colwell, K.S.; Goddard, W.A.; Skiff, W.M. UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations. J. Am. Chem. Soc.; 1992; 114, pp. 10024-10035. [DOI: https://dx.doi.org/10.1021/ja00051a040]

26. Jansen, A.P.J. An Introduction to Kinetic Monte Carlo Simulations of Surface Reactions; Springer: Berlin, Germany, 2012.

27. Buck, A.L. New equations for computing vapor pressure and enhancement factor. J. Appl. Meteorol.; 1981; 20, pp. 1527-1532. [DOI: https://dx.doi.org/10.1175/1520-0450(1981)020<1527:NEFCVP>2.0.CO;2]

28. Haynes, W.M. CRC Handbook of Chemistry and Physics; CRC Press: Boca Raton, FL, USA, 2014.

29. Liang, T.; Sawyer, W.G.; Perry, S.S.; Sinnott, S.B.; Phillpot, S.R. Energetics of oxidation in MoS2 nanoparticles by density functional theory. J. Phys. Chem. C; 2011; 115, pp. 10606-10616. [DOI: https://dx.doi.org/10.1021/jp110562n]