Early Oxidation Processes on the Greigite Fe 3 S 4 (001) Surface by Water: A Density Functional Theory Study

: Greigite (Fe 3 S 4 ), the sul ﬁ de counterpart of the spinel-structured oxide material magnetite (Fe 3 O 4 ), is a mineral widely identi ﬁ ed in anoxic aquatic environments and certain soils, which can be oxidized, thereby producing extremely acid solutions of sulfur-rich wastewaters, so-called acid mine drainage (AMD) or acid rock drainage (ARD). Here we report a computational study of the partial replacement of sulfur (forming H 2 S) by oxygen (from H 2 O) in the Fe 3 S 4 (001) surface, derived from density functional theory calculations with on-site Coulomb approach and long-range dispersion corrections (DFT+ U − D2). We have proposed three pathways for the oxidation of the surface as a function of H 2 O coverage and pH. Di ﬀ erent pathways give di ﬀ erent intermediates, some of which are followed by a solid-state di ﬀ usion of the O atom. Low levels of H 2 O coverage, and especially basic conditions, seem to be essential, leading to the most favorable energetic landscape for the oxidation of the Fe 3 S 4 (001) surface. We have derived the thermodynamic and kinetic pro ﬁ le for each mechanism and plotted the concentration of H 2 S and protons in aqueous solution and thermodynamic equilibrium with the stoichiometric and partially oxidized Fe 3 S 4 (001) surface as a function of the temperature. Changes in the calculated vibrational frequencies of the adsorbed intermediates are used as a means to characterize their transformation. We have taken into account statistical entropies for H 2 S and H 2 O and other experimental parameters, showing that this mineral may well be among those responsible for the generation of AMD.


INTRODUCTION
Extremely acidic sulfur-rich wastewaters are a current worldwide problem. So-called acid mine drainage (AMD) or acid rock drainage (ARD) is associated with natural weathering of rock formations 1 but aggravated, in particular, by existing and historic human activities such as the mining industry. 2−4 Once mining or processing operations expose metal sulfide compounds, in particular pyrite (cubic FeS 2 ), to weathering elements such as O 2 and H 2 O 5 as well as certain microorganisms, 6−8 the minerals steadily oxidize. This process results in reduction of the water pH 9 and potentially high concentrations of toxic metallic and metalloid elements in solution, 10 depending on the initial composition of the exposed minerals.
A number of chemical reactions leading to oxidation and dissolution of metal sulfides have been studied, and several mechanisms have been suggested, depending on the minerals and the oxidizing agent present. One of the most remarkable mechanisms, called the polysulfide mechanism, is characteristic of acid-soluble sulfides, such as sphalerite (ZnS), where protons attack the mineral and produce H 2 S. Sulfur-oxidizing bacteria further oxidize the H 2 S to SO 4 2− regenerating the protons. 6 Thiosulfate, the other important mechanism characteristic of nonacid-soluble persulfides, such as FeS 2 , 11 proceeds initially by an attack of aqueous Fe 3+ ions to the metal sulfide mineral, thereby generating protons, thiosulfate (S 2 O 3 2− ) and Fe 2+ . This Fe 2+ is then reoxidized by iron-oxidizing bacteria, while S 2 O 3 2− decomposes into elemental sulfur and SO 4 2− ions. 6−8 Iron sulfides are the predominant sulfides found in anoxic marine sediments 12 and therefore one of the main causes of AMD. 5,13,14 FeS 2 is a stable mineral in these environments 12,15,16 and considerable research effort has been devoted to study its oxidation, 7,8,17−23 but other usually coexisting iron sulfides have been largely overlooked. Among other minerals, the spinel-structured greigite (Fe 3 S 4 ) is a metastable intermediate in the formation pathway of FeS 2 from mackinawite (tetragonal FeS), 24−28 which has a long environmental persistence 29 and can also be found in aquatic environments 30−35 and soils 36,37 as well as in magnetotactic bacteria 38 and gastropods. 39 Fe 3 S 4 solubility has been rationalized mostly as a digestion, ignoring other important oxidation processes that can take place in aqueous medium. It has been suggested that, during digestion, the mineral forms 0.75 parts of H 2 S and disproportionates to Fe 2+ and orthorhombic sulfur. 40 However, elemental sulfur was not formally identified even after exposing Fe 3 S 4 to a mixture of mineral acids and reducing agents and was therefore assumed to be the sulfur content not recovered (25% of S in Fe 3 S 4 ). In another study of the dissolution of Fe 3 S 4 in distilled water under an H 2 S atmosphere, orthorhombic sulfur was again considered as one of the products even though it was not detected. 41 This assumption was based on the fact that the oxidation potential of this process was close to that of the halfcell involving H 2 S and elemental sulfur.
Purity of the samples is another issue affecting the accurate measurement of Fe 3 S 4 solubility, which is usually difficult to control due to the unavoidable copresence of FeS and FeS 2 . For example, less sulfur can be recovered from Fe 3 S 4 aged for 1 week than from freshly prepared samples as it slowly oxidizes to pyrite, requiring harsher conditions to dissolve. 40 So far, due to the metastable nature of this mineral, any attempt to quantify Fe 3 S 4 solubility has been difficult to reproduce, with highly uncertain and low sulfur recoveries. 42 The analysis of Fe 3 S 4 in a mixture of iron sulfides has also been found to be inaccurate as the dissolution stoichiometry is difficult to control. 40,43−48 Moreover, the apparent value measured for Fe 3 S 4 by Berner 41 could be artificially increased, since each iron sulfide phase present has a different solubility. 49 Given the importance of Fe 3 S 4 as a contributor in the formation of AMD and the experimental inaccuracies in the determination of its solubility, in this paper we present a fresh description of the early oxidation stages of this mineral. Taking into account the structural differences between FeS 2 and Fe 3 S 4 and the decisive role of the persulfide group in dictating the mechanism of the oxidation reactions of metal sulfides, 6 we have used density functional theory (DFT) calculations to investigate the early oxidation processes of Fe 3 S 4 via a polysulfide mechanism, in order to explain this mineral's lability in a disturbed aqueous medium. We propose three mechanisms to account for the partial replacement of S by O atom in the top layer of the Fe 3 S 4 (001) surface, which is the most prominent plane in solvothermal synthesized nanoparticles. 50 We have provided the simulated infrared (IR) spectra of the intermediates along the oxidation mechanisms. We have also applied thermodynamic arguments to examine the pH conditions and the H 2 S concentration, in aqueous solution and temperatures at which these species are in equilibrium with the stoichiometric and oxidized Fe 3 S 4 (001) surface, which is relevant to the geochemical formation of AMD.

Calculation Details.
We have performed spinpolarized calculations with the Vienna ab initio simulation package (VASP). 51−54 All simulations were carried out within the periodic plane-wave DFT framework. The projector augmented wave (PAW) method was used to describe the core electrons and their interaction with the valence ones. 55,56 The frozen core of the Fe, S, and O elements was defined up to and including the 3p, 2p, and 1s electrons, respectively. At the level of the generalized gradient approximation (GGA), the exchange-correlation in the form of the Perdew−Wang 91 (PW91) 57,58 functional was used together with the spin interpolation of Vosko et al. 59 The long-range dispersion interactions were added via the semiempirical method of Grimme (D2), 60 using the global scaling factor parameter optimized for the Perdew−Burke−Ernzerhof (PBE) 61,62 functional, s 6 = 0.75, which has shown to provide good results in the modeling of a number of sulfides. 63−68 Brillouin zone integrations of the surface slabs were performed using a Γcentered Monkhorst−Pack grid 69 of 4 × 4 × 1 k-points. The isolated H 2 O and H 2 S molecules were modeled in a cell with broken symmetry which was big enough to avoid intermolecular interactions, sampling only the Γ point of the Brillouin zone. In order to increase the integration efficiency in the reciprocal space, the partial occupancies for all calculations were determined using the tetrahedron method with Blochl corrections. 70 Kohn−Sham (KS) valence states were expanded in a plane-wave basis set with the kinetic energy's cutoff fixed at 600 eV. The DFT + U 71 version of Dudarev et al. 72 was used for the description of the localized and strongly correlated d Fe electrons. On the basis of previous work, we have chosen a U eff of 1.0 eV. 50,73−75 Electronic density optimization was stopped when the total energy difference between two consecutive selfconsistent loop steps was below 10 −5 eV. Atomic positions were relaxed to their ground state using the conjugate-gradient method and were considered converged when the Hellmann− Feynman forces on all atoms were smaller than 0.02 eV·Å −1 . The dimer method was used to search the transition states (TS), 76,77 which were characterized by a single imaginary frequency along the reaction coordinate. Higher cutoff values and k-point grids, as well as a lower self-consistent energy threshold, were tested to ensure energies were converged within 1 meV per atom. Fe 3 S 4 has a spinel crystal structure characterized by the space group Fd3̅ m. 35 The face-centered cubic unit cell is composed of 32 sulfur anions, which are almost regularly close packed along the [111] direction. This S arrangement generates 8 tetrahedral (A) and 16 octahedral (B) holes per unit cell occupied by Fe cations, see Figure 1, part a. Fe 3 S 4 is a 2−3 type spinel, which classification is based on the cation oxidation state. This spinel has an inverse cation distribution  The Journal of Physical Chemistry C Article were set antiparallel, with the Fe ions in high-spin states, in line with previous studies. 73,74 2.2. Surface Models. Surface calculations were carried out using a symmetric and nondipolar stoichiometric slab. 83 This well-defined surface, which is only an approximate model of a real surface, is essential for calculating accurate and realistic surface properties. 84−88 The geometry optimized topmost layer of the (001) surface terminates in a bulk-like structure containing single rows in the [110] direction of 5-coordinated Fe B ions alternating every two single rows of S ions with cubic packing, see Figure 1, parts b and c. Beneath this layer, there are 1.5 monolayers (ML) of 4-coordinated Fe A also forming rows parallel to the S ones. From this layer, 0.5 ML of Fe A with a (√2 × √2)R 45°symmetry originates from above the surface, which moved inward during relaxation. The surface was modeled with a slab composed of 8 formula units of Fe 3 S 4 and exposing a surface area of 93.53 Å 2 . The slab model was formed by eight atomic layers, in which the five bottom-most layers were kept frozen at their relaxed bulk positions, while the remaining top layers were allowed to optimize. A vacuum of 12 Å was added between the periodic slabs in the direction perpendicular to the surface. Convergence of the slab energy within 1 meV per atom was further tested with different slab and vacuum thickness as well as a number of relaxed layers. In order to enhance the electronic convergence, dipole corrections perpendicular to the surface were included in our simulations to account for any dipole created by the chemical species added in the relaxed surface of the slab. 89,90 A Bader charge analysis was used to partition the electron density grid into atomic charges. 91−93 Charge density flow diagrams were constructed by subtracting from the electronic charge density of the total adsorbate−surface system, the sum of the charge densities of the isolated adsorbate and clean surface in the same geometry.
For further information, the optimized lattice (a) and sulfur parameters (u), the total magnetization of saturation (M S ), the atomic spin moments (m s ) and the atomic charges (q) of bulk Fe 3 S 4 are shown in Table SI 2.3. Simulation of the Infrared Spectra. Vibrational frequencies were determined for the isolated H 2 O and H 2 S molecules and for the reaction intermediates on the Fe 3 S 4 (001) surface for each state of each mechanism. We have used central a The presented vibrational modes are the asymmetric stretching (ν asym ), symmetric stretching (ν sym ) and bending (δ) modes. The Journal of Physical Chemistry C Article finite differences, by allowing each ion to move in the direction of each Cartesian coordinate by a small negative and positive displacement, to calculate the Hessian matrix. The prediction of the vibrational frequencies from the second derivative of the potential energy with respect to the atomic positions takes into account only the curvature at the minimum of the well. In this region of the energy surface, the potential has a harmonic behavior, meaning that the vibrational energy levels are equally spaced. In contrast, the fundamental modes observed experimentally usually correspond to the transition between the ν = 0 and the ν = 1 vibrational energy levels of an anharmonic potential energy surface. To compensate for the approximate nature of the calculated vibrational frequencies, they are usually corrected with an empirical scale factor. 94−98 We have reported in Table 1 the calculated values for the fundamental vibrational modes (unscaled and scaled) of the isolated H 2 O and H 2 S molecules. The quality of the calculated quantities can be assessed through comparison with the experimental values for the gas phase molecules, 99 also shown in Table 1. The scaling factor is calculated as c = ∑(ω exp · ω calc )/∑ω calc 2 , while its uncertainty is obtained from u = (∑(ω calc 2 (c − ω exp /ω calc ) 2 ))/∑ω calc 2 , where ω calc and ω exp are the calculated and experimentally observed IR wavenumbers, respectively. 94 The value of c = 0.98 ± 0.01 for our set of two molecules indicates a good predictive behavior of the PW91 functional for vibrational frequencies. The largest difference between the unscaled (scaled) calculated wavenumbers and experiments is 130.2 (−43.8) cm −1 and the smallest is 4.2 (19.3) cm −1 , while the uncertainty of the largest unscaled wavenumber is 38.9 cm −1 . Taking into account that the PW91 functional provides sufficiently accurate IR frequencies, we are confident of our predicted vibrational wavenumbers for the adsorbed intermediates presented in this study. 94−98 In subsequent sections, only the unscaled wavenumbers are presented.
We have only allowed the displacement of atoms with dangling bonds and adsorbed species for the simulation of the IR spectra. The coupling between the adsorbates and surface phonons, all of which appear below 400 cm −1 , was neglected. The intensities of the fundamental IR bands, the square root of the dynamic dipole moments, were determined by projecting the derivative of the dipole moment onto the basis of normal modes. Each vibrational band was smeared with a Gaussian function.
2.4. Thermodynamic and Kinetic Profiles. We have modeled three mechanisms for the initial stages of the Fe 3 S 4 (001) surface oxidation at various H 2 O coverage and pH regimes. Note that in our approximated model we have neglected the water environment and the zero point energy corrections. The relative energy of the system (ΔE j ) for every state j along the proposed thermodynamic and kinetic profile with respect to the energies of the pristine relaxed Fe 3 S 4 (001) surface slab (E slab ) and an isolated H 2 O molecule in vacuum (E H 2 O ) was calculated according to the following equation, where n represents the number of H 2 O molecules. The adsorption energy (E ads i ) of the species i is determined from the difference between ΔE j of the adsorbed state j and the previous one j − 1. The dissociation energy (E diss i ) is obtained from the difference between ΔE j of the dissociated state j and the state j − 2. The activation energy (E Aj ) is calculated from the difference between ΔE j of the transition state j and the intermediate of the state j − 1. The energy required to substitute one S atom or SH group in the mineral surface by an The Journal of Physical Chemistry C Article adsorbed OH group (E exch O→S ) is derived from the difference between ΔE j of the states after and before substitution.

RESULTS AND DISCUSSION
3.1. Oxidation of the Fe 3 S 4 (001) Surface. We have considered three main mechanisms for studying the oxidation of the Fe 3 S 4 (001) surface through the partial exchange of S from the top layer by O from the water producing H 2 S, see Figure 2, part a, Figure 3, part a, and Figure 4, part a. For pathway 1 (2) Figure 5. Although all Fe B are equivalent in this surface, the existence of the Fe A ions and the newly absorbed H 2 O molecule generates four nonequivalent types of S ions around the hydroxylated Fe B ion. The dissociated H lies between the S next to the hydroxylated Fe B and at the tetrahedral cavity of the Fe A row in the [110] direction. The charge analysis of state 2 indicates that, upon dissociation, an overall charge of 0.50 e − has been transferred, mainly from the Fe B to the OH group and from the protonated S to form a covalent bond with the H. This can be represented graphically in a charge density difference plot, illustrating electron density changes on the dissociative adsorption of a H 2 O molecule. Figure 6 shows that charge is mostly localized on the O−H group and along the S− H bond. At state 2, the stretching mode (ν) of the adsorbed OH and incrusted SH groups are red-shifted with respect to the asymmetric (ν asym ) and symmetric (ν sym ) stretching modes of the respective isolated molecules. This indicates the weakening of the O−H bond, which together with the bending mode (δ) of the H−S−Fe and H−O−Fe groups support the H 2 O dissociation on the Fe 3 S 4 (001) surface, see Figure 7. We have found that configuration 2 is 0.34 eV above the reference energy (the pristine Fe 3 S 4 (001) plus an isolated H 2 O molecule). We have also tested unsuccessfully the possibility of complete H 2 O splitting owing to the higher basicity of O 2− than S 2− . This behavior is similar to that found on the (001) surface of the isomorphic FeNi 2 S 4 . 66 The reaction profile for the partial exchange of S by an O atom on the Fe 3 S 4 (001) surface, according to mechanism 1, is shown in Figure 2. In this mechanism, only one H 2 O molecule per unit cell is involved, equivalent to a coverage of 0.25 ML. Note that we have defined a full ML as the amount of H 2 O The Journal of Physical Chemistry C Article required to hydrate all four Fe B ions in the top surface layer per unit cell. In a concerted step, the SH group migrates to the top of the nearby Fe B ion in the same layer, while O replaces S in the structure, thereby increasing its coordination number by three Fe ions. This is an exothermic process with an E exch O→S of −0.24 eV and an activation energy (E A3 ) of 0.80 eV. Upon exchange of SH by OH, the O atom lies at 0.63 Å below the average position of the surface S atoms, in part due to its smaller ionic radius, see state 4 in Figure 5. Note that 0.44 Å is the difference in Shannon's radii between the hexacoordinated ions, S 2− and O 2− . 101 The O−H distance remains at 0.98 Å, with the H lying at the same level of the top S layer on the cavity of the Fe A row along the [110] direction. The charge analysis shows that during substitution of the labile SH by the more nucleophilic OH group, there is a transfer of 0.10 e − from the adsorbed S to the structural O. Figure 6 shows the charge transfer mechanism at the state 4 with respect to the oxidized surface and dissociated H 2 S. In agreement with the atomic electronegativities, the electronic density located along the O− H bond comes primarily from the H, whereas the adsorbed SH relocates its electrons to the Fe B −S and S−H bonds. At state 4, the vibrational mode of the OH group is red-shifted by 56.7 cm −1 and the SH group is blue-shifted by 118.8 cm −1 with respect to the IR bands of state 2, see Figure 7, parts a and b.  The next step in this mechanism is the generation of H 2 S following the migration of the adsorbed H, see state 6 in Figure  2. In the resulting structure, the O remains slightly below the top S layer, but moves horizontally by 0. 81 102 According to the charge analysis, the formation of the H 2 S molecule is accompanied by a reduction of the charge of the migrating H by 0.51 e − , which is provided by the SH and O−Fe B −S groups (where S is next to O along the [1̅ 1̅ 0] direction). Note that the fact that the charges of the two H atoms become very similar is an indication of the weak adsorption of the H 2 S molecule, in agreement with the geometric similarities to the gas phase molecule. The step from state 4 to 6 is endothermic by 0.45 eV and its activation energy (E A5 ) is 0.29 eV smaller than E A3 indicating that the OH−SH exchange is the determinant step in this pathway. At state 6, the O−H stretching mode disappears, in agreement with the dissociation of this bond, see Figure 7, part a. We have been able to detect the H 2 S bending along with the asymmetric and symmetric stretching modes, corroborating the formation of this molecule, see Figure 7, parts b and c. The stretching bands are red-shifted with respect to the values calculated for the isolated molecule, suggesting that H 2 S is interacting with the surface. H 2 S desorption is an unfavorable process, as state 7 is 1.06 eV above the reference and 0.51 eV above state 6, which affects negligibly the surface structure and the charge of the atoms of the surface and H 2 S. In the reverse process, the dissociation of a H 2 S molecule adsorbed on a partially oxidized Fe 3 S 4 (001) has an activation energy of 0.06 eV, just above the thermal energy (2k B T) at 298 K and therefore is thermodynamically and kinetically more favorable.
According to this pathway, the position for the S substitution was dictated by the most favorable H adsorption site. Nevertheless, the O moves to coordinate the Fe A next to the cavity along the [110] direction in a solid state diffusion process that stabilizes the slab by an additional 0.14 eV after overcoming a barrier of 2.14 eV. This high activation energy for O diffusion suggests that this step is highly unlikely to take place under static conditions. In the structure of this final configuration, the O atom moved toward the vacuum by 0.21 Å but still remained below the S top layer. The charge analysis reveals that, after O diffusion, there has been 0.29 e − charge transfer from the Fe B ions (from the surface and subsurface layers) formerly coordinating the O to the surface Fe B and Fe A ions currently coordinating the O atom, consistent with an oxidation state of 3+.
3.1.2. Pathway 2. A sensible alternative is to consider the H migration away from the Fe B −OH center. This situation leads us to pathway 2, shown in Figure 3, which considers the partial oxidation of the surface by OH groups (coverage of 0.25 ML). Note that in agreement with the reference, the energies of states 2 to 6 have been calibrated by considering the energy of one proton located onto the most stable position on the Fe 3 S 4 (001) surface. We present the structures, charge density flow and IR spectra of the intermediates for pathways 2 and 3 in the Supporting Information, due to their similarities to pathway 1.
The hydroxylated surface, see stage 2 in Figure 3, is 0.09 eV lower in energy than the equivalent stage of pathway 1, corroborating the thermodynamic feasibility of the H migration. The Fe B −O bond length is 0.02 Å shorter, see Figure SI-2, suggesting a slightly stronger interaction between these atoms in the absence of the second H atom, in agreement with the charge analysis shown in Figure SI In the absence of a protonated sulfur, the hydroxyl O replaces the S, leading to the most thermodynamically favorable product. The energy of this state is 0.15 eV lower than the reference system, in contrast to the 0.10 eV higher in pathway The Journal of Physical Chemistry C Article 1. The activation energy required to take our system from state 2 to 4 is 0.17 eV lower than in pathway 1; see Figure 3. The replaced S atom forms a hydrogen bond with the hydroxyl H; see state 4 in Figure SI-2. The O−H bond distance increases by 0.02 Å, while the O is displaced outward from the surface by 0.07 Å, as well as along the [110] direction, compared to the same stage in pathway 1. The Fe B holding the adsorbed S atom is also moved along the [11̅ 0] direction. The charge analysis shows a rearrangement involving mainly the adsorbed S, the incorporated O and the Fe A ion coordinating it; see Figure SI The next step along the oxidation process is the formation of the SH group. The protonated S moves to rest almost straight atop the Fe B ion, which relaxes to its original position; see state 6 in Figure SI This suggests that its dipole moment is different than in pathway 1, which agrees with the SH group interacting more strongly with the surface.
In the last step of pathway 2, the SH group interacts with a coadsorbed H atom, leading to the formation and desorption of the H 2 S molecule, which leaves a partially oxidized surface slab in its most favorable geometry. We have also considered a variation of pathway 2, where the exchange of S by O takes place after the OH dissociation; see Figure 3. Nevertheless, the energies of the intermediates are higher than in the original pathway and therefore it is not representative as it is less likely to lead to Fe 3 S 4 (001) surface oxidation.
3.1.3. Pathway 3. Next, we have investigated the presence of a second H 2 O molecule (0.50 ML) on the oxidation mechanism, see Figure 4. Previous work on the analogue spinel FeNi 2 S 4 has suggested that the adsorption of a second H 2 O molecule after dissociation of the first is most likely to take place on the Fe B ion next to the protonated S, see state 2 in Figure SI-5. 66 The bond angle of the newly added H 2 O is increased by 3.4°with respect to the calculated and experimental value for the isolated molecule. 102 The O−H bonds of the second H 2 O molecule are 0.97 and 1.04 Å long. While the shortest length indicates no bond alteration, 102 the elongated OH bond is forming a hydrogen-bond with the neighboring hydroxyl group. The addition of a second H 2 O molecule stabilized significantly the system, making this configuration the one with the lowest energy in this study.
In order to generate a second SH group within close proximity to a subsurface Fe A ion, we have stretched the second H 2 O until dissociation; see state 4 in Figure SI-5. The two SH groups are bridging both hydroxylated Fe B ions, a configuration similar to the one found on FeNi 2 S 4 . 66 Hence, the new S−H bond is nearly perpendicular to the surface. The OH···OH hydrogen bond is stretched, leading to a reduction of the intrabond length. The configuration with two adsorbed OH groups and two surface SH is less thermodynamically stable than state 2. The charge difference for the migrating H is similar to the value calculated in previous pathways, see Figure  SI The protonated S binding the subsurface Fe A is the most energetically favorable to be exchanged by an OH. After the exchange, the Fe A ion holding the adsorbed SH group has migrated outward from the surface by 1.23 Å occupying an octahedral position; see state 6 in Figure SI-5. The structural OH group lies below the topmost atomic layer as in the two previously proposed pathways. An OH···OH distance of 1.38 Å indicates a hydrogen-bond between these OH groups, in agreement with the charge density flux plot of state 6 in Figure SI-6. The energy released during the SH to OH exchange step and the activation energy are the largest among the three pathways studied.
In the next step toward the Fe 3 S 4 (001) oxidation, the formation of a coadsorbed H 2 O molecule takes place as a result of the deprotonation of the incorporated OH; see state 8 in Figure 4. This process is endothermic and requires the lowest activation energy barrier, making it the less demanding step among the three pathways. After the H migration, the surface O atom moves inward by 0.09 Å, while the H 2 O moves away from the surface, almost presenting itself as an isolated molecule, 102 see Figure  In the final step of the H 2 S formation, one of the H atoms from the H 2 O may migrate to the SH group, see state 10 in Figure 4. As this step is energetically and kinetically unfavorable, the SH groups are likely to remain in the surface, leading to their accumulation, similar to pathway 2. At this stage, the S of the H 2 S molecule moves 0.18 Å outward and the Fe B −O distance is reduced by 0.29 Å; see Figure SI The simulated IR spectra for pathway 3 are more complex than previous ones due to a larger number of intermediates, see Desorption of the H 2 S molecule is an endothermic process, see state 11 in Figure 4. The release of the H 2 S molecule is accompanied by the formation of a hydrogen-bond between the OH group and the surface O atom, see Figure SI-5. However, apart from this new interaction, the impact of the H 2 S desorption on the surface structure is minimal.
The O diffusion process further stabilizes the system by 0.40 eV, i.e. nearly three times more exothermic than the solid diffusion process of pathway 1. However, the barrier to the transition state in pathway 3 is 6.6 times higher than in pathway 1, making this solid state transformation step less likely than in pathway 1. In the final configuration, the O occupies a position comparable to the one in the previously discussed pathways. There is, however, a lasting impact of the oxidation through the final mechanism on the position of the labile Fe A ion, which remains in the octahedral cavity.
The Journal of Physical Chemistry C Article

Thermodynamics of H 2 S and its Ionization
Products in Aqueous Solution. Because of the interest in preventing H 2 S release from iron−sulfur compounds, which leads to the acid mine drainage, we have investigated the equilibrium concentration of these species in aqueous solution. We have assumed that our system behaves ideally, which allows us to replace thermodynamic equilibrium constants by concentration quotients and activities by molal concentrations. We have calculated the concentration of aqueous H 2 S from the oxidation of the Fe 3 S 4 (001) surface and the pH of this solution as a function of temperature, in the range from 293 to 373 K. In order to calculate the concentrations of the species in aqueous solution in equilibrium with the pristine and partially oxidized Fe 3 S 4 (001) surfaces, we have considered the process of partial oxidation of this surface by a H 2 O molecule, according to eq 2, where the equilibrium constant (K o1 ) is equal to the [H 2 S] in equilibrium with the solid phases.  Table 2. We have not considered the second dissociation of H 2 S as its constant K a2 103 is at least 8 orders of magnitude smaller than K a1 104 in the temperature range considered here. We have also ignored the autodissociation of H 2 O as its ionic product (K W ) is also at least 6 orders of magnitude smaller than K a1 in the range of temperatures of interest. 105 Consequently, the [H + ] is controlled by the first dissociation of H 2 S, whereas the contribution of H + from the other processes is negligible.
In order to link the gas phase state of the isolated molecules in our DFT simulations with the states in eq 2, we have  Table 2.
The combination of eq 2, 4 and 5, leads to the reaction of the partial oxidation of the Fe 3 S 4 (001) surface, eq 6 where H 2 O and H 2 S are in the reference state of our DFT calculations: K o is then calculated as shown in eq 7, from the Gibbs free energy (G) of the partial oxidation of the Fe 3 S 4 (001) surface, where R and T are the ideal gas constant and the temperature of interest, respectively.
The change in the standard Gibbs free energy (ΔG ⊖ ) for the partial oxidation of the Fe 3 S 4 (001) surface, according to eq 6, was calculated from ΔG ⊖ = ΔH ⊖ − TΔS ⊖ , where ΔH ⊖ is the enthalpy of this process, per H 2 S molecule formed, ΔS ⊖ is the change in entropy and T is the temperature. The enthalpies of the three pathways investigated in this study are obtained directly from our calculations, where we assume that the enthalpy values do not depend on the temperature. We have additionally assumed that the entropies of the solid phases remain largely unchanged through the oxidation process and     The partial oxidation of the Fe 3 S 4 (001) surface in a wet environment, represented by the concentration of H 2 S as a function of the temperature, is shown in Figure 8, part a. This graph only shows the curve associated with pathway 2 (ΔH ⊖ = 0.92 eV), which is the most kinetically and thermodynamically favorable, since the curves related to the other the pathways are almost coincident with this one. As expected from an endothermic process, the increment of temperature brings an exponential increase of the concentration of products: [H 2 S] is equal to 2.77·10 −3 mol·kg −1 at T = 293 K and 3.65·10 −2 mol· kg −1 at T = 393 K. Figure 8, part b, shows how the pH decreases with temperature, as is also expected from a solution with an increasing concentration of a weak acid. In Figure 8, parts a and b, we have also added two auxiliary lines showing the behavior of hypothetical pathways whose ΔH ⊖ are approximately ±10 times bigger than pathway 2. These lines illustrate how noticeably the enthalpy of the process affects the dependence of [H 2 S] and pH with T. The trend of pH with respect to temperature shown in Figure 8, part b, is therefore also valid for any sulfide that dissolves in H 2 O through an oxidation mechanism with an enthalpy within the range −10.0 to +10.0 eV. As the experimental pH = 4.03 of a solution of Fe 3 S 4 at 25°C 41 apparently increases by the effect of the more soluble coexisting FeS phase, 49 our results agree with this expectation. Figure 8, part b, shows that the initial oxidation stages of pure Fe 3 S 4 by H 2 O will produce a solution with pH = 4.7 at 25°C. Since the pH is always below 4.81 in the range of temperatures from 293 to 373 K and large Fe 3 S 4 deposits are present in aquatic environments, from our calculations it would appear that these mineral deposits could be significant contributors to acid mine drainage.

CONCLUSIONS
In this paper, we have proposed and modeled three different pathways for the early steps of the Fe 3 S 4 (001) surface oxidation promoted by H 2 O. We have used DFT methods with a Hubbard Hamiltonian and empirical long-range dispersion corrections to calculate the reactants, products, intermediates and transition states in the thermodynamic and kinetic energy profiles. In each of these mechanisms, a surface S atom is replaced by an O atom from H 2 O, producing a H 2 S molecule. We have found that the step where the OH group replaces the S atom always take place before the OH deprotonation. However, for pathway 2, we could also model an alternative reaction route, but there all intermediates and transition states were between 1.56 and 0.37 eV higher in energy than in the former mechanism. We have found that pathways 1 and 2 are the most effective routes to initiate the Fe 3 S 4 (001) surface oxidation. Although in pathway 1 the total ΔH ⊖ (per H 2 S molecule) is 0.14 eV higher than in pathway 2, the ratedetermining step in pathway 1 is 0.41 eV smaller than in pathway 2. The presence of the dissociated H + from the H 2 O in the vicinity of the reactive site in pathway 1 directs the oxidation reaction toward a kinetic product where the solid state exchange of S by O becomes the step with the highest activation energy. On the other hand, the absence of the dissociated H + in the vicinity of the reactive site in pathway 2 leads to the most thermodynamically favorable oxidized state and the O−H dissociation becomes the rate-limiting step. When two H 2 O molecules are adsorbed near the reactive center, pathway 3, they direct the oxidation reaction to a kinetic product 0.42 eV higher in energy than the one obtained in pathway 1. The highest activation energy calculated in this study is associated with the dissociation of the second H 2 O molecule in pathway 3, making this an unlikely process. For any of the proposed mechanisms, H migration (i.e., from H 2 O or OH) and H 2 S desorption are endothermic processes, while the introduction of the OH group on the surface is an exothermic step. The energy balance of these processes, especially for pathway 2 operating under basic conditions, makes the partially oxidized Fe 3 S 4 (001) surface likely to remain hydroxylated as this is the most stable intermediate. Changes in the values of the wave numbers of the calculated normal vibrational modes of the adsorbed intermediates indicate their transformation. The calculated [H 2 S] in aqueous solution, and therefore pH, in chemical equilibrium with the solid phases at a range of temperatures, shows that Fe 3 S 4 may be among those minerals responsible for acid mine or rock drainage.
Future work will include the consideration of the Fe 3 S 4 oxidation by O 2 and a mixture of H 2 O and O 2 . In addition, we The Journal of Physical Chemistry C Article will also carry out a microkinetic analysis of all the steps taking place during the initial oxidation stages of Fe 3 S 4 , which will allow us to include temperature and concentration effects and will mimic the dissolution of this mineral within a realistic processes framework.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.6b00216.
Table showing the experimental and optimized lattice and sulfur parameters, total magnetization of saturation, and the atomic spin moments and atomic charges of bulk