# Direct prediction of isotopic properties from molecular dynamics trajectories: Application to sulfide, sulfate and sulfur radical ions in hydrothermal fluids

Marc Blanchard, Elsa Desmaele, Gleb S Pokrovski, Carlos Pinilla, Merlin

Méheut, Rodolphe Vuilleumier

► To cite this version:

Marc Blanchard, Elsa Desmaele, Gleb S Pokrovski, Carlos Pinilla, Merlin Méheut, et al. Direct prediction of isotopic properties from molecular dynamics trajectories: Application to sulfide, sulfate and sulfur radical ions in hydrothermal fluids. Chemical Geology, 2024, 661, pp.122202. 10.1016/j.chemgeo.2024.122202. hal-04744419

HAL Id: hal-04744419

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

Submitted on 18 Oct 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Direct prediction of isotopic properties from molecular dynamics trajectories: Application to sulfide, sulfate and sulfur radical ions in hydrothermal fluids

Marc Blanchard<sup>a, *</sup>, Elsa Desmaele<sup>b</sup>, Gleb S. Pokrovski<sup>a</sup>, Carlos Pinilla<sup>c</sup>, Merlin Méheut<sup>a</sup>, Rodolphe Vuilleumier<sup>b</sup>

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

Check for
updates

## ARTICLE INFO

Keywords:

Sulfur isotopes

Sulfur radical ions

Sulfate

Sulfide

Molecular dynamics

Density functional theory

Hydrothermal fluid

## ABSTRACT

In hydrothermal fluids, disulfur ($S_2^{2-}$) and trisulfur ($S_3^{2-}$) radical anions have been observed to coexist with the major hydrogen sulfide and sulfate species. These radical ions have potentially important effects on the solubility, transport and fractionation of metals and sulfur, with consequences for ore deposit formation and, more generally, for geochemical cycles of metals and volatiles. It is therefore essential to know the intrinsic isotopic properties of these important sulfur species in order to use sulfur isotopes for tracing different geological processes. Here, the theoretical equilibrium isotopic properties of the disulfur and trisulfur radical ions are computed and compared to the hydrogen sulfide ($H_2S$) and sulfate ion ($SO_4^{2-}$), using, for the first time, a first-principles molecular dynamics (MD) approach. The isotopic properties are calculated directly from molecular dynamics trajectories using the vibrational density of states and the atomic kinetic energy, and then compared to the more established method based on sampling of several snapshots. This comparison allowed us to validate the new modelling method and to assess its advantages and limitations. The predicted equilibrium isotope fractionation in terms of $^{34}$S/$^{32}$S between $S_2^{2-}$ and $S_3^{2-}$ is small, i.e. <1‰, with a slight enrichment in the heavier isotope for $S_3^{2-}$, over the temperature range 200–500 °C. Both radical ions are slightly depleted in the heavier isotope, by 1 to 2‰, relative to aqueous $H_2S$. Our results help tuning sulfur isotope fractionation models used for tracing the origin and evolution of hydrothermal fluids. Our method opens large perspectives for using the rapidly growing body of MD simulation data in geosciences on structure and stability of aqueous complexes to assess in parallel element isotope fractionations from MD-generated trajectories.

## 1. Introduction

Sulfur isotopes ($^{32}$S, $^{33}$S, $^{34}$S, and $^{36}$S) are among the ‘traditional’ stable isotopes most widely used in geochemistry whose study began in the 1940s (e.g., Thode et al., 1949). Since then, sulfur isotopes have received numerous applications for tracing a plethora of geological processes from the mantle to the Earth’s surface (e.g., Thode et al., 1961; Ohmoto and Rye, 1979; Canfield, 2001; Johnston, 2011). Most hydrothermal, magmatic and metamorphic systems involve sulfur-bearing fluids, which are at the origin of metal sulfide ore. The understanding of the processes and mechanisms controlling mineral resources on Earth requires robust knowledge of sulfur speciation in such fluids and the

associated isotope fractionations. The speciation of sulfur in geological fluids has long been believed to be controlled by sulfide ($H_2S$, $HS^-$ and $S^{2-}$), sulfate ($HSO_4^-$ and $SO_4^{2-}$), and sulfur dioxide ($SO_2$, mostly in vapor phase). The large differences in oxidation state among these major sulfur forms (from -2 to +6) are mainly responsible for the large sulfur isotope fractionations observed in nature. Other intermediate-valence sulfur species such as sulfite ions ($HSO_3^-$ and $SO_3^{2-}$), polysulfide dianions ($S_n^{2-}$) and their protonated counterparts ($HS_n^-$ and $H_2S_n$) have been invoked, in complement to $HS^-$ and $H_2S$ in fluids and vapors within limited temperature-pressure-compositional ranges, as potential transporting agents for sulfur and chalcophile metals (e.g., Seward, 1973; Berndt et al., 1994; Pokrovski et al., 2009, 2014; Walters et al., 2020; Hu et al.,

<footnote><sup>a</sup> Géosciences Environnement Toulouse (GET), Université de Toulouse, CNRS, IRD, UPS, CNES, 14, avenue Edouard Belin, 31400 Toulouse, France</footnote>

<footnote><sup>b</sup> PASTEUR Laboratory, Département de Chimie, Ecole Normale Supérieure, PSL University, Sorbonne Université, CNRS, 24 rue Lhomond, 75005 Paris, France</footnote>

<footnote><sup>c</sup> Departamento de Física y Geociencias, Universidad del Norte, Barranquilla, Colombia</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: marc.blanchard@get.omp.eu (M. Blanchard).</footnote>
2022; references therein), but their stability and abundance versus the major sulfate and sulfide forms remain to be quantified. Recently, additional intermediate-valence species, the disulfur (S₂⁻) and trisulfur (S₃⁻) radical anions, have emerged in aqueous fluids at elevated temperatures and pressures as demonstrated by in-situ Raman spectroscopy both in nature and experiment (Pokrovski and Dubrovinsky, 2011; Jacquemet et al., 2014; Truche et al., 2014; Pokrovski and Dubessy, 2015; Barré et al., 2017; Schmidt and Seward, 2017; Colin et al., 2020). In contrast to the 'traditional' intermediate-valence sulfur forms mentioned above, the radical ions appear to be stable across a wide range of pressure and temperature, from 150 to at least 800 °C and 30 kbar. They had gone unnoticed because they undergo very fast breakdown to sulfide, sulfate and native sulfur on fluid cooling to ambient conditions, at which their equilibrium concentrations in water become very small. Since this discovery, the physicochemical properties and role of these radical ions in geological processes have been explored. For example, it has been shown that the radical ions partition at least 10 times more than sulfate and hydrogen sulfide into the fluid phase during silicate melt degassing at depth (Colin et al., 2020). Sulfur radicals also enhance the solubility of metals like gold, platinum and palladium in hydrothermal fluids with consequences for the formation of ore deposits (Pokrovski et al., 2015, 2021a, 2022; Laskar et al., 2022). Furthermore, trisulfur ion complexes with iron may act as precursors during pyrite formation in hydrothermal systems (Pokrovski et al., 2021b). The investigation of their isotopic properties has also begun. Kokh et al. (2020) experimentally explored the role played by the sulfur radical ions on sulfur isotope signatures of molecular sulfur (S⁰) and pyrite in hydrothermal systems. These authors showed that S₃⁻, S⁰, and H₂S have similar isotope compositions, within 1‰ of δ³⁴S values, at temperatures between 300 and 450 °C, but did not find any significant mass independent fractionation (Δ³³S < 0.1‰). From a theoretical point of view, Tossell (2012) investigated the stability of S₃⁻ and its complexes with Cu⁺. Using quantum chemistry calculations based on molecular clusters, Tossell (2012) reported the first equilibrium isotopic fractionation factors between S₃⁻, S₂⁻ and other sulfur species. His theoretical work was followed by others such as Eldridge et al. (2021) who focused on the isotopic properties of a large number of sulfur species relevant to hydrothermal environments, including radical ions and polysulfide diions. In their study, the solvation of sulfur-bearing molecules was considered explicitly through encapsulation in a static cluster of water molecules.

Computational chemistry approaches have become efficient tools to investigate isotopic fractionations (e.g., Blanchard et al., 2017 for a review). While minerals are relatively simple to treat due to their long-range order, modelling of solvated chemical species represents a certain degree of difficulty. To go beyond static calculations, which are based on the construction of a molecular cluster and performed at the athermal limit, it is necessary to employ molecular dynamics techniques to accurately capture the effect of temperature and the fluctuation in time of the molecular geometry. However, once the molecular dynamics trajectory is acquired at a given temperature, the isotopic fractionation is usually determined by sampling snapshots for which the atomic positions are relaxed at zero temperature before calculating the vibrational frequencies – a necessary ingredient to obtain the isotopic fractionations. Although this method has been validated in a number of studies (e.g., Kowalski and Jahn, 2011; Dupuis et al., 2015; Pinilla et al., 2015; Ducher et al., 2018), it is not totally satisfactory. This is because it requires heavy additional calculations following the molecular dynamics, and, especially, it erases significant contributions from the configurational disorder as well as from the temperature effect.

To overcome these fundamental limitations, here we employed an alternative method to directly calculate the isotopic fractionation from molecular dynamics trajectories. Our method is based on the atomic kinetic energy, which is derived from the vibrational density of states of the element of interest. Such an approach has recently been applied to study iron isotopes in metallic alloys relevant to the Earth's core (Pinilla

et al., 2021) as well as in silicate melts (Rabin et al., 2023). Here, for the first time, we evaluated this method in the case of hydrothermal fluids. The structural, vibrational and mass-dependent isotopic properties of the radical ions S₂⁻ and S₃⁻ were determined and compared to the better-known hydrogen sulfide (H₂S) and sulfate ion (SO₄²⁻). The latter has been chosen as a proxy of the sulfate isotopic system, even though SO₄²⁻ and HSO₄⁻ species have nearly equal fractions in most near-neutral hydrothermal fluids at the calculation conditions. Furthermore, there is no direct experimental work measuring the isotopic fractionation between HSO₄⁻ and SO₄²⁻ separately, so the common assumption is that the fractionation between the two species is within error. This is what we see, for example, between the experimental work of Syverson et al. (2015) with 100–1000 times more HSO₄⁻ than SO₄²⁻ at pH = 2, and of Kokh et al. (2020) with a predominance of SO₄²⁻ and NaSO₄ species at pH = 4–8. These two studies reported the same total fractionation between sulfate and H₂S, within a couple of permils.

## 2. Methods

### 2.1. First-principles molecular dynamics simulations

Four aqueous sulfur species (S₃⁻, S₂⁻, SO₄²⁻, and H₂S) were modeled by first-principles molecular dynamics (FPMD). Cubic boxes containing these species together with 128 water molecules were built with a volume adjusted in order to reach densities typical to natural hydrothermal fluids within the shallow crust (<5 km). The resulting cubic box sizes are 18.51 Å and 18.56 Å at 673 K for the two sulfur radical species, S₂⁻ and S₃⁻ respectively, while the box sizes for H₂S and SO₄²⁻ are 17.61 Å and 17.75 Å at 573 K. The two latter species were simulated at a lower temperature at which their stability is generally increased with respect to the radical ions (Pokrovski and Dubessy, 2015). Periodic boundary conditions were employed. The simulations were carried out with the CP2K code package (VandeVondele et al., 2005; Kühne et al., 2020) in the framework of the density functional theory (DFT) as implemented in the QUICKSTEP module, following the same approach and parameters as in a previous study on gold speciation in the presence of sulfur radical species (Pokrovski et al., 2015). The BLYP exchange-correlation functional (Becke, 1988; Lee et al., 1988) was used in combination with a Van der Waals correction DFT-D3 (Grimme et al., 2010). A plane-wave cutoff of 600 Ry for the electronic density of S₂⁻ and S₃⁻ and 800 Ry for H₂S and SO₄²⁻ were chosen as well as a triple-zeta valence doubly polarizable (TZV2P) basis set (VandeVondele and Hutter, 2007). The interaction between the ionic cores and the valence electrons was treated with the Goedecker-Teter-Hutter (GTH) norm-conserving pseudo-potentials (Goedecker et al., 1996). In systems with an odd number of electrons, spin polarization was taken into account. A neutralizing background charge was implicitly added for all charged systems.

The Born-Oppenheimer molecular dynamics was carried out in the NVT ensemble with a time step of 0.5 fs, and the system was thermostated by means of a CSVR thermostat (Bussi et al., 2007) with a time constant of 100 fs for simulation of the solvated radical species S₂⁻ and S₃⁻, and with a time constant of 500 fs for the lower temperature simulations of solvated H₂S and SO₄²⁻. For hydrogen, the mass of the heavy isotope deuterium was chosen, which improves the convergence of the simulation (e.g. Pokrovski et al., 2015, 2022), except for the H₂S simulation where it is important to have H–S chemical bonds with the right masses in order to derive the accurate sulfur isotopic properties. The simulations were run for at least 18 ps.

### 2.2. Isotope fractionation

The enrichment in the heavier isotope of one individual phase or aqueous species relative to another phase or species coexisting at equilibrium is characterized by the isotope fractionation factor *a*. This fractionation factor can be determined from the measured isotope ratios
of each species, but can also be calculated using first-principles methods. In the latter case, the equilibrium fractionation factor α (expressed in ‰) is obtained from the β-factors of each species (A and B)

$$1000\ln\alpha_{A-B} = 1000\ln\beta_{A} - 1000\ln\beta_{B} \quad (1)$$

Rigorously, this β-factor must itself be obtained from the reduced partition function ratios of each isotopologue of the chemical species under study and their relative abundances (Richet et al., 1977). In general, the β-factor is approximated by the geometric mean of the reduced partition function ratios of all the singly-substituted isotopologues, but it has been shown that a single calculation with all the sites of the atom of interest substituted at the same time leads to the same result (Méheut et al., 2007). It is important to note, however, that this approach does not apply to hydrogen isotopes (Méheut et al., 2007). In this study, we computed the β-factor directly from the FPMD trajectory by using the following expression based on the atomic kinetic energy, K (Polyakov, 1991)

$$\ln\beta = \frac{m - m'}{m} \left( \frac{K}{RT} - \frac{3}{2} \right) \quad (2)$$

where m and m' are the masses of the heavy and light isotopes respectively. In the present study, the method, hereafter referred to as “VDOS” method, is applied to sulfur stable isotopes (four in number: $^{32}$S, $^{33}$S, $^{34}$S and $^{36}$S) and results are given for the mass ratio 34/32. Rigorously, in eq. (2), K is the kinetic energy of the isotope m' ($^{32}$S in this case). However, considering the average mass of the sulfur in the simulation does not affect the results. This approach therefore corresponds to a simultaneous substitution of all the S sites and allows the β-factor to be obtained directly. The kinetic energy is obtained by integrating the vibrational density of states (VDOS) of the element of interest. It is computed from the Fourier transform of the velocities of atom i (from the space of time (t) to the space of frequencies (ω)):

$$v_i(\omega) = \int_0^{\omega_{tot}} v_i^a(t) e^{i\omega t} dt \quad (3)$$

where $t_{tot}$ is the total simulation time. The VDOS $f_n(\omega)$ for element n is then proportional to the sum of $|v_i(\omega)|^2 = \bar{v}_i(\omega)v_i(\omega)$, where $\bar{v}_i$ is the complex conjugate of $v_i$ and for i spanning the atoms of element n:

$$f_n(\omega) = \frac{1}{N_n} \sum_i^{N_n} \sum_{\alpha=(x,y,z)} \frac{|v_i(\omega)|^2}{T} \quad (4)$$

Following this relation, the average kinetic energy for element n is finally expressed as:

$$K = \frac{m_n}{k_B T} \int_0^{\omega_{cut}} e(\hbar\omega) f_n(\omega) d\omega \quad (5)$$

where $\omega_{cut}$ is a frequency at which all vibrational states have been integrated (we arbitrarily chose 100 THz, which is beyond all vibrational contributions for the four investigated systems), $m_n$ the mass of element n, and $e(\hbar\omega)$ is the Einstein function for the vibrational energy of a single harmonic oscillator at frequency $\omega$:

$$e(\hbar\omega) = \hbar\omega \left( \frac{1}{2} + \frac{1}{\exp(\hbar\omega/k_B T) - 1} \right) \quad (6)$$

In order to improve convergence, we set the temperature in eq. (5) as the temperature measured from the kinetic energy of element n (not the target temperature):

$$k_B T = \frac{1}{3} m_n \int_0^{\omega_{cut}} f_n(\omega) d\omega \quad (7)$$

The standard deviation associated with the calculation of the kinetic energy was estimated using a block averaging method. For blocks of different length, $\tau_b$, we calculate $K(\nu_{cut})$. The standard deviation as a function of block length is then expressed as:

$$\epsilon(\tau_b) = \frac{\sqrt{\frac{1}{n_b} \sum_{b \text{ in blocks}} (K_b - \langle K_b \rangle)^2}}{\langle K_b \rangle} \quad (8)$$

Isotopic fractionation is controlled by the bonding properties of the element of interest with its molecular environment. Therefore, another quantity that is commonly investigated is the interatomic force constant, which is linearly related with the β-factor (e.g. Dacher et al., 2018). The average force constant, $\langle F \rangle$, can also be obtained by integrating the VDOS of the sulfur atom:

$$\langle F \rangle = m_n \int_0^{\omega_{cut}} \frac{\omega^2 f_n(\omega) d\omega}{\int_0^{\omega_{cut}} f_n(\omega) d\omega} \quad (9)$$

The average force constant derived from this method will be referred as $\langle F_{VDOS} \rangle$. The standard deviation associated to the calculation of the force constant was estimated using the same block averaging method as for the kinetic energy.

Additionally, we determined the average interatomic force constant from a more conventional approach based on instantaneous force constants determined for more than fifteen snapshots evenly sampled along the FPMD trajectory. This approach will be referred to as “Snapshot” method and the average force constant will be noted $\langle F_{SNAP} \rangle$. In this case, the finite displacement method was used. The force constant is determined from the second derivative of the total energy with respect to the displacement of the atom in three orthogonal directions. This is equivalent to the first derivative of the forces acting on this atom with respect to its displacements. The atom of interest was displaced by <0.01 Å starting from the atomic positions of the snapshot without preliminary geometry optimization.

# 3. Results

## 3.1. Structural properties

The average geometries (bond lengths and bond angles) of the four aqueous species are given in Fig. 1 and compared with previous experimental and theoretical data in Table 1. The corresponding radial distribution functions are displayed in Fig. S1 (Supplementary Material). Our results show similar average S–S bond lengths for S$_3^-$ and S$_3^+$ . With 2.06 ± 0.05 Å, the average S–S bond length of S$_3^-$ is slightly longer here than in the ab initio molecular dynamics calculations of Mei et al. (2013), i.e. 2.02 ± 0.01 Å at 573 K and 2 GPa with the PBE functional, as well as in the athermal static calculations performed in Eldridge et al. (2021) where the sulfur species are embedded in water clusters, i.e. 2.03 ± 0.01 Å with the B3LYP functional. It should be noted that the large uncertainty associated with our measurements compared to other studies reflects the dynamic nature of our simulations. For comparison, Tossell (2012) performed static athermal calculations at several different quantum mechanical levels in vacuum or with an implicit

<figure><img src="image_3.png" /><figcaption>Fig. 1. Structures of stable S species from FPMD simulations (S$_2^-$, S$_3^-$, H$_2$S and SO$_4^{2-}$ from left to right) whose S, O and H atoms are shown in yellow, red and grey colour, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<caption>**Table 1**</caption>

Geometries (average bond lengths and angles with the corresponding 1σ standard deviations) obtained in this study for the four aqueous sulfur species compared with a selection of previous theoretical (Calc.) and experimental (Exp.) values; NA = not applicable.

<table><thead><tr><th>Species</th><th>Method, reference</th><th>Bond length (Å)</th><th>Bond angles (degree)</th></tr></thead><tbody><tr><td rowspan="3">S<sub>2</sub><sup>•−</sup></td><td>Calc. This study</td><td>S-S</td><td>NA</td></tr><tr><td>Calc. Eldridge et al. (2021)</td><td>2.06 ± 0.05</td><td>NA</td></tr><tr><td></td><td>2.032 ± 0.001</td><td>NA</td></tr><tr><td rowspan="4">S<sub>3</sub><sup>•−</sup></td><td>Calc. This study</td><td>S-S</td><td>S-S-S</td></tr><tr><td>Calc. Eldridge et al. (2021)</td><td>2.06 ± 0.05</td><td>114.9 ± 7.9</td></tr><tr><td>Calc. Mei et al. (2013)</td><td>2.031 ± 0.011</td><td>114.3 ± 1.1</td></tr><tr><td>Calc. Tossell (2012)</td><td>2.017 ± 0.004</td><td>112.9 ± 6.0</td></tr><tr><td rowspan="5">H<sub>2</sub>S</td><td>Calc. This study</td><td>S-H</td><td>H-S-H</td></tr><tr><td>Calc. Eldridge et al. (2021)</td><td>1.36 ± 0.02</td><td>92.5 ± 6.6</td></tr><tr><td>Calc. Eldridge et al. (2016)</td><td>1.347 ± 0.002</td><td>92.9 ± 0.5</td></tr><tr><td>Calc. Otake et al. (2008)</td><td>1.349 ± 0.001</td><td>92.6</td></tr><tr><td>Exp. Herzberg (1966)</td><td>1.355</td><td>93.2</td></tr><tr><td rowspan="6">SO<sub>4</sub><sup>2−</sup></td><td>Calc. This study</td><td>S-O</td><td>O-S-O</td></tr><tr><td>Calc. Eldridge et al. (2021)</td><td>1.52 ± 0.04</td><td>109.3 ± 5.9</td></tr><tr><td>Calc. Eldridge et al. (2016)</td><td>1.52 ± 0.01</td><td>109.5 ± 1.2</td></tr><tr><td>Calc. Otake et al. (2008)</td><td>1.52 ± 0.01</td><td>109.5 ± 1.0</td></tr><tr><td>Exp. McGinnety (1972)</td><td>1.512</td><td>109.5</td></tr><tr><td>Exp. McGinnety (1972)</td><td>1.486</td><td>NA</td></tr></tbody></table>

solvation model, and found a S-S distance of 1.99 Å for the calculation assumed as the most accurate one (CCSD pvQZ calculation). Both in present and previous calculations, S-S-S angles are similar within uncertainty (Table 1).

For the H₂S and SO₄²⁻ species, we found average bond lengths of 1.36 Å and 1.52 Å, respectively. These values, as well as the average bond angles, are very close to the results obtained from the thermal static calculations done with localized basis sets and B3LYP functional (Otake et al., 2008; and Eldridge et al., 2016, 2021, treating the effects of solvation either with the Polarizable Continuum Model or with explicit water molecules). Experimental bond lengths determined on gaseous compounds (Herzberg, 1966; McGinnety, 1972) are ~0.03 Å shorter than any of these calculated data.

The radial distribution functions (Fig. S1) also indicate that the four sulfur-bearing species are surrounded by 6 to 12 water molecules forming a hydration shell with their oxygen atoms located on average at 3.5–3.8 Å from the central sulfur atom.

## 3.2. Vibrational properties

Vibrational densities of state (VDOS) of sulfur, obtained from the Fourier transform of sulfur atomic velocities and thus capturing the anharmonicity related to the dynamical behavior of the system, are displayed in Fig. 2 for the four aqueous species. The average frequencies of the main vibrational densities of state (encompassing infrared and Raman active modes) are compared with theoretical and experimental available data in Table 2.

Overall, our calculated frequencies are slightly underestimated compared to the experimental Raman frequencies, for both bending and stretching modes. This underestimation, related to the DFT functional used, reaches 12% of the absolute Raman wavenumber value for S₂⁻ and S₃⁻ while H₂S displays the best agreement with an underestimation of only 2%. Previous theoretical studies (Otake et al., 2008; Tossell, 2012) employing the B3LYP functional in thermal static calculations could obtain frequencies closer to the experimental values. It should be noted that these calculated frequencies are harmonic, whereas the experimental frequencies are fundamental and contain both harmonic and anharmonic contributions. These theoretical frequencies use sometimes scaling factors (× 0.89–0.97) of relative amplitudes similar to the underestimations above.

## 3.3. Kinetic energy and force constants

The frequency-dependence of the normalized kinetic energy of sulfur obtained from the integration of the VDOS of each sulfur species (using Eq. 5) is shown in Fig. 2. Values of the normalized kinetic energies at the $v_{cut}$ frequency of 100 THz are reported in Table 3, associated to the standard deviation (2σ) estimated using a block averaging method. Similarly, the sulfur force constant is derived from the integration of the VDOS using Eq. 9 (Fig. S2) and the corresponding values are reported as 〈FVDOS〉 in Table 3.

Following the alternative “Snapshot” method, sulfur force constants were also computed for a set of evenly sampled snapshots and then averaged (Fig. 3). The average of instantaneous force constants called 〈F<sub>SNAP</sub>〉 is reported in Table 3. The method based on the sampling of snapshots is commonly used for investigating the isotopic properties of liquid systems (e.g., Kowalski and Jahn, 2011; Dupuis et al., 2015; Pinilla et al., 2015; Dacher et al., 2018).

## 3.4. Isotopic properties

Fig. 4 shows the sulfur β-factors of the four aqueous species calculated from the sulfur kinetic energy using Eq. 2 (i.e. from the “VDOS” method). Values of β-factors are calculated over the 473–773 K temperature range surrounding the target temperatures of the MD simulations. The equations of the temperature-dependence polynomial fits drawn in Fig. 4 are given in Table 4. The temperature extrapolation of the β-factor using Eq. 2 from molecular dynamics simulations conducted at a given temperature has been studied in detail in the work of Pinilla et al. (2021) on iron isotope fractionation in solid and liquid Fe–S alloys. Over the 300 K range considered here, the error due to extrapolation can be considered negligible compared to the error bars intrinsic to the “VDOS” method, as estimated above.

The present results confirm the common knowledge that when isotopic equilibrium is reached, the SO₄²⁻ species will be strongly enriched in the heavier S isotopes compared to H₂S (~30–10‰ over the 473–773 K temperature range). This large fractionation is primarily due to the contrasting sulfur redox states (S²⁻ versus S⁶⁺) as well as the different nature of the interatomic bonds of S with H and O. These differences in interatomic bonding result in a significant contribution of the vibrational frequencies around 1000 cm⁻¹ to the isotope fractionation of SO₄²⁻ (Fig. 2) due to the much larger sulfur force constant for SO₄²⁻ than for H₂S (Table 3). On the other hand, β-factors of the S₂⁻ and S₃⁻ radical species are lower by 1–2‰ than that of H₂S over this T range, with the following order of increasing β-factor values S₂⁻ < S₃⁻ < H₂S, with the error bars partly overlapping.

# 4. Discussion

## 4.1. Comparison between the “VDOS” and “Snapshot” methods

Fig. 5 compares sulfur interatomic force constants obtained from the “VDOS” method using Eq. (9) and from the “Snapshot” method (corresponding data are in Table 3). The two methods give the same force constant value (within error bars) for S₂⁻, S₃⁻ and SO₄²⁻. However, for H₂S, the “Snapshot” method gives a force constant value of 44% larger than the “VDOS” method. The mass of deuterium was chosen for the hydrogen atoms of all systems, with the exception of H₂S. This choice has no impact on the observation made here. Fig. S3 compares the radial distribution functions of a preliminary trajectory of D₂S + 128 D₂O, with that of the system with the mass of hydrogen. The water molecules are distributed in the same way around H₂S in both cases. The average bond lengths are identical for H₂S and D₂S (i.e. 1.36 Å). As a result, the force constant of S obtained by the “snapshot” method (F<sub>SNAP</sub>) is identical in both cases: 280 ± 28 N/m and 283 ± 15 N/m for the systems with the mass of hydrogen and deuterium, respectively. It should be noted, however, that the “VDOS” method requires the use of hydrogen mass to
<figure><img src="image_4.png" /><figcaption>Fig. 2. Vibrational density of state (VDOS) of sulfur for the four aqueous species (S₂⁻, S₃⁻, H₂S and SO₄²⁻ from top left to bottom right) and the normalized kinetic energy of sulfur (KE/k_BT) as a function of frequency up to the νcut frequency of 100 THz.</figcaption></figure>

obtain the correct force constant (FVDOS) and S β-factor. This method is based on vibrational frequencies, which are much lower for D₂S than for H₂S. From this and other preliminary results, we learnt that the “VDOS” method is strongly affected by the quality of the local temperature equilibration. If the temperature of the atoms of the H₂S molecule and thus of the H—S stretching vibrations is lower than the average temperature of the system, then the force constant obtained from the vibrational density of states is decreased. These observations suggest that the “Snapshot” method is more robust than the “VDOS” one. The “Snapshot” method is only sensitive to the atomic configurations whereas the “VDOS” method requires a good description of the configurations but also of the equipartition of the energy. The temperature set from the kinetic energy of S atoms, Eq. (7), only corrects for the lack of equipartition in an isotropic way, while the H-S bond is very directional. The latter method therefore requires a long and perfectly equilibrated MD trajectory in the case of H₂S.

Despite the precautions to have properly equilibrated trajectories, the final MD trajectory of H₂S shows significantly different “VDOS” and “Snapshot” force constants. Note that if we consider the last 20 ps of the trajectory rather than the 30 ps, the “VDOS” force constant goes from 194 ± 45 to 224 ± 40 N/m. The latter value remains within the

uncertainty but gets closer to the value of 280 ± 28 N/m given by the “Snapshot” method. This slight discrepancy leads us to estimate the isotopic fractionations associated with the H₂S force constant obtained by the “Snapshot” method.

Here, the snapshots of the MD trajectory were not relaxed before calculating the instantaneous force constants by the finite displacement method. This procedure allows us both to keep the configurational disorder and to partly take into account the anharmonicity contribution. However, to obtain the corresponding β-factors, it would be necessary to calculate the vibrational frequencies for each snapshot that has been previously relaxed in order to avoid the appearance of imaginary frequencies linked to structural instabilities. These structural relaxations and frequency calculations represent a considerable amount of additional work and, at the same time, erases part of the disorder.

Another possibility is to use the approximate formula derived from Bigeleisen and Mayer (1947) in order to obtain the β-factor directly from the force constant:

$$ \beta \approx 1 + \frac{m - m'}{mm'} \frac{\hbar^2 F}{8(kT)^2} \qquad (10) $$
<caption>**Table 2**</caption>

Vibrational frequencies (in cm⁻¹) obtained in this study from the main VDOS bands of the four aqueous sulfur-bearing species (Fig. 2) compared with a selection of previous theoretical and experimental values.

<table><thead><tr><th>Species</th><th>Method, reference</th><th>Wavenumber (cm<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="3">S<sub>2</sub><sup>•</sup></td><td>Calc. VDOS (this study)</td><td>523</td></tr><tr><td>Calc. Raman (Tossell, 2012)</td><td>557 (ss)</td></tr><tr><td>Exp. Raman (Pokrovski and Dubessy, 2015)</td><td>585 (ss)</td></tr><tr><td rowspan="3">S<sub>3</sub><sup>•</sup></td><td>Calc. VDOS (this study)</td><td>215</td></tr><tr><td>Calc. Raman (Tossell, 2012)</td><td>220 (b)</td></tr><tr><td>Exp. Raman (Pokrovski and Dubrovinsky, 2011)</td><td>240 (b)</td></tr><tr><td rowspan="3">H<sub>2</sub>S</td><td>Calc. VDOS (this study)</td><td>1165</td></tr><tr><td>Calc. Raman (Otake et al., 2008)</td><td>1168 (b)</td></tr><tr><td>Exp. Raman (Pokrovski and Dubessy, 2015)</td><td>ND</td></tr><tr><td rowspan="3">SO<sub>4</sub><sup>2-</sup></td><td>Calc. VDOS (this study)</td><td>–</td></tr><tr><td>Calc. Raman (Otake et al., 2008)</td><td>406 (b)</td></tr><tr><td>Exp. Raman (Pokrovski and Dubessy, 2015)</td><td>ND</td></tr></tbody></table>

* Data from Tossell (2012) correspond to the B3LYP CBSB7 PCM model, while data from Otake et al. (2008) correspond to the scaled frequencies obtained with the B3LYP 6-311G(d,p) PCM model. The experimental data for S<sub>2</sub><sup>•</sup> and S<sub>3</sub><sup>•</sup> were acquired at elevated temperatures (> 300 °C, 0.01–0.1 GPa) while for H<sub>2</sub>S and SO<sub>4</sub><sup>2-</sup>, the spectra were collected in the temperature range 200–500 °C in liquid-like fluids. ND = not detected, b = bending, ss = symmetric stretching, as = asymmetric stretching.

<caption>**Table 3**</caption>

Length of the MD trajectory (t in ps) used for data production (i.e. the 3 first picoseconds of initial equilibration were neglected for S<sub>2</sub><sup>•</sup>, S<sub>3</sub><sup>•</sup> and SO<sub>4</sub><sup>2-</sup>, as well as the 10 first picoseconds for H<sub>2</sub>S), average temperature of the MD simulations (⟨T⟩ in K), normalized kinetic energy of sulfur (KE/k<sub>B</sub>T), and average interatomic force constants of sulfur (⟨F<sub>VDOS</sub>⟩ and ⟨F<sub>SNAP</sub>⟩ in N/m) obtained either from the “VDOS” or “Snapshot” method. KE/k<sub>B</sub>T and ⟨F<sub>VDOS</sub>⟩ results are given at a frequency of 100 THz, i.e. an arbitrary frequency chosen beyond all vibrational contributions for the four investigated systems.

<table><thead><tr><th>Species</th><th>t</th><th>⟨T⟩</th><th>KE/k<sub>B</sub>T</th><th>⟨F<sub>VDOS</sub>⟩</th><th>⟨F<sub>SNAP</sub>⟩</th></tr></thead><tbody><tr><td>S<sub>2</sub><sup>•</sup></td><td>30.0</td><td>674</td><td>1.523 ± 0.006</td><td>77 ± 21</td><td>99 ± 18</td></tr><tr><td>S<sub>3</sub><sup>•</sup></td><td>15.5</td><td>669</td><td>1.532 ± 0.004</td><td>107 ± 19</td><td>106 ± 14</td></tr><tr><td>H<sub>2</sub>S</td><td>30.0</td><td>573</td><td>1.556 ± 0.012</td><td>194 ± 45</td><td>280 ± 28</td></tr><tr><td>SO<sub>4</sub><sup>2-</sup></td><td>30.0</td><td>574</td><td>1.861 ± 0.015</td><td>948 ± 45</td><td>968 ± 88</td></tr></tbody></table>

where *ħ* = *ħ*/2π, *m* and *m'* are the atomic masses of the two isotopes. Fig. 6 compares the β-factors of the four sulfur species calculated from the vibrational density of state and the atomic kinetic energy using Eq. 2 with the β-factors obtained from the “VDOS” force constant using Eq. 9 and 10. Eq. 10 is a high-temperature approximation diverging from the VDOS approach at low temperature. The amplitude of this divergence depends on the nature of the elements forming bonds with sulfur. The lighter these elements are, the bigger is the divergence. Results indeed show that this high-temperature approximation is valid for S–S bonds in the S<sub>2</sub><sup>•</sup> and S<sub>3</sub><sup>•</sup> radicals, while the divergence goes beyond the error bars for the S–O bonds of the SO<sub>4</sub><sup>2-</sup> species and is particularly large for the S–H bonds of H<sub>2</sub>S. Therefore, this approximation cannot be used here to estimate the H<sub>2</sub>S β-factor corresponding to the “Snapshot” force constant.

The alternative solution to the above high-temperature approximation is to take advantage of the one-to-one relationship between β-factors and interatomic force constants, as demonstrated for instance by Ducher et al. (2016). Thus, this linear correlation has been established here between sulfur β-factors and force constants, both using the “VDOS” method, at three temperatures covering the investigated range (Fig. 7A). It was then possible to apply these linear correlations to the H<sub>2</sub>S force constant derived from the “Snapshot” method, to obtain the corresponding β-factor and its temperature dependence (Fig. 7B). This complementary result will be used in the next section to discuss the isotopic data.

## 4.2. Isotopic properties

Fig. 8A compares the present results for the equilibrium isotopic fractionation between sulfate and sulfide (i.e. 1000 ln β(SO<sub>4</sub><sup>2-</sup>) – 1000 ln β(H<sub>2</sub>S)) with a selection of previous theoretical and experimental data.

Ohmoto and Lasaga (1982) estimated the equilibrium fractionation factor from a large set of available hydrothermal experiments of partial exchange of sulfur isotopes between aqueous sulfate (represented by HSO₄⁻, SO₄²⁻ ions and their Na⁺ ion pairs) and hydrogen sulfide conducted before the 1980s. Syverson et al. (2015) applied the three-isotope technique to multiple S isotope data during experiments of pyrite precipitation under hydrothermal conditions at acidic conditions, where HSO₄⁻ and H₂S coexist (300–350 °C, 500 bars). Kokh et al. (2020) performed quadruple S isotope analyses in the system S<sub>3</sub><sup>•</sup>-S<sub>2</sub><sup>•</sup>-S<sup>0</sup>-sulfatesulfide-bearing aqueous fluid - pyrite in a wide pH range where SO₄²⁻, NaSO₄⁻, H₂S, HS⁻ were present (300–450 °C, 300–800 bar). With regard to the atomistic modelling work, different approaches were taken. The studies of Otake et al. (2008), Tossell (2012), Eldridge et al. (2021) all correspond to static athermal non-periodic calculations but differing in the choice of the functionals, basis sets, solvation models and scaling factors for vibrational frequencies. For example, Otake et al. (2008) found, as best choice, the HF functional and the use of scaling factors. Tossell (2012) selected the B3LYP functional and applied an implicit solvation model (Polarizable Continuum Model). Eldridge et al. (2021) also used the B3LYP functional but adopted an explicit solvation by building “manually” different water molecule clusters around the sulfur species and taking the average sulfur β-factor value from several clusters. Our approach differs from those studies as it is the only one based on first-principles molecular dynamics (BLYP functional with a Van der Waals correction). Therefore, it enables capturing the dynamical disorder of the aqueous system at a given temperature, and partially accounting for anharmonicity related to deviations of the potential energy from equilibrium. Despite a good description of bond lengths and angles (Table 1; our values are comparable to those of Otake et al., 2008 and Eldridge et al., 2021), as well as vibrational frequencies (Table 2), our isotope fractionation values are lower (at the extreme limit of error bars) than the previous experimental and theoretical data set. Note also that the use of the “Snapshot” force constant for H₂S results in a smaller fractionation, which is further away from the experimental data.

Fig. 8B compares the isotopic fractionations S<sub>3</sub><sup>•</sup> - H₂S, and S<sub>2</sub><sup>•</sup> - H₂S. No experimental studies having directly measured equilibrium fractionation are available, because the radical ions cannot be preserved on fluid cooling. The only existing so far experimental estimate is that from the study of Kokh et al. (2020) who measured isotope fractionation between H₂S and native sulfur - the latter being produced by S<sub>3</sub><sup>•</sup> breakdown on fluid cooling. Their data suggest S<sub>3</sub><sup>•</sup> and H₂S to have similar isotopic signatures within 1‰ under hydrothermal conditions (300–450 °C), which is in good agreement with our “VDOS” method estimations. We can further compare our results with the theoretical estimates of Tossell (2012) and Eldridge et al. (2021). For these sulfur species, our data, despite the significant intrinsic uncertainty associated
<figure><img src="image_5.png" /><figcaption>Fig. 3. Evolution of the sulfur interatomic force constant calculated from snapshots sampled evenly along the MD trajectory without preliminary geometry optimization (thin line connecting the symbols) and cumulative average (bold line) for the four aqueous species.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Temperature-dependence of the sulfur β-factors for the four aqueous species, calculated using Eq. 2.</figcaption></figure>

with dynamic variability of the system, do frame the available theoretical data. More specifically, the “VDOS” method predicts a slightly lower equilibrium fractionation factors than the previous data while the fractionation corresponding to the “Snapshot” method is a little more pronounced (i.e. the difference between the “VDOS” and “Snapshot” methods equals 2.8‰ at 300 °C, and 1.7‰ at 450 °C). Another important

<caption>Table 4</caption>

Polynomial fits of $ ^{34}\text{S}/^{32}\text{S} $ β-factors as a function of temperature (473–773 K range) for the four aqueous species: 1000 $\ln \beta = ax + bx^2$ with $x = 10^6/T^2$ (T in K), and values of the β-factors (in ‰) at 573 K and 673 K.

<table><thead><tr><th>Species</th><th>a</th><th>b</th><th>1000 ln β at 573 K</th><th>1000 ln β at 673 K</th></tr></thead><tbody><tr><td>S<sub>2</sub><sup>-</sup></td><td>0.62219</td><td>-5.40 × 10<sup>-3</sup></td><td>1.84 ± 0.48</td><td>1.34 ± 0.35</td></tr><tr><td>S<sub>3</sub><sup>-</sup></td><td>0.85490</td><td>-6.88 × 10<sup>-3</sup></td><td>2.54 ± 0.32</td><td>1.86 ± 0.24</td></tr><tr><td>H<sub>2</sub>S</td><td>1.33312</td><td>-8.16 × 10<sup>-2</sup></td><td>3.27 ± 0.71</td><td>2.54 ± 0.55</td></tr><tr><td>SO<sub>4</sub><sup>2-</sup></td><td>7.52465</td><td>-1.80 × 10<sup>-1</sup></td><td>21.22 ± 0.88</td><td>15.74 ± 0.65</td></tr></tbody></table>

result of our study is the equilibrium fractionation between S<sub>2</sub><sup>-</sup> and S<sub>3</sub><sup>-</sup>. According to the study by Eldridge et al. (2021) and the present work, S<sub>3</sub><sup>-</sup> would be slightly enriched in the heavier isotope compared with S<sub>2</sub><sup>-</sup>, with an equilibrium fractionation of 0.44‰ in our work and 0.33‰ in Eldridge et al. (2021) study at 723 K. Only the value of S<sub>2</sub><sup>-</sup> estimated by Tossell (2012) is outside this trend with a β-factor of 6.1‰ smaller than that of S<sub>3</sub><sup>-</sup> at the same temperature.

The S–S bond lengths of S<sub>2</sub><sup>-</sup> and S<sub>3</sub><sup>-</sup> are similar (Fig. 1), and their vibrational frequencies are roughly comparable (Fig. 2). The interatomic force constants are close (Table 3) and the isotopic fractionation between the two species is weak. The analysis of each individual S atom of S<sub>2</sub><sup>-</sup> and S<sub>3</sub><sup>-</sup> shows that the inner S of S<sub>3</sub><sup>-</sup> has a slightly larger force constant than the outer S (i.e. 148 ± 20 N/m for the inner S vs 83 ± 13 N/m and 88 ± 20 N/m for the two outer S), whereas both S atoms in S<sub>2</sub><sup>-</sup> have a force constant similar to that of the outer S in S<sub>3</sub><sup>-</sup> (i.e. 99 ± 19 N/m). These data are in perfect agreement with the analysis of the reduced partition function ratio of each S atom by Eldridge et al. (2021) (i.e. 30.9‰ and 18.4‰ at 25 °C for the inner and outer S of S<sub>3</sub><sup>-</sup> to be compared to 19.0‰ for S<sub>2</sub><sup>-</sup>). It is therefore logic to conclude that the isotopic fractionation between S<sub>2</sub><sup>-</sup> and S<sub>3</sub><sup>-</sup> is weak, with a slight enrichment of S<sub>3</sub><sup>-</sup> in the heavier isotope related to the inner S atom.
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Comparison of sulfur interatomic force constants (in N/m) obtained from the “VDOS” or “Snapshot” methods (Table 3). The dashed line represents the 1:1 relationship.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Temperature-dependence of the sulfur β-factors for the four aqueous species. Values of β-factor is calculated either from the VDOS and atomic kinetic energy using Eq. 2 (bold lines) or from the “VDOS” force constant using Eqs. 9 and 10 (dashed lines).</figcaption></figure>

### 4.3. Geochemical implications

Our results extend the knowledge of sulfur radical species fractionation factors that are particularly hard to measure experimentally because of the unquenchable nature of these species in hydrothermal fluids (e.g., Kokh et al., 2020). Interestingly, our calculations predict “lighter” isotope signatures for both ions compared to H₂S at equilibrium (by 1 to 2‰ at 300–400 °C), which does not follow the general intuitive relationship that the heavier isotope must enrich species of sulfur in the higher oxidation state (formal oxidation states of S atoms in

$H_2S$, $S_3^-$ and $S_2^+$ are –2, –0.33 and –0.5, respectively). This deviation from the formal oxidation state trend is due to combined factors from different degree of covalency in sulfur bonding with different atoms and variable sulfur coordination and solvation degree.

The small but significant differences in the β-factors between H₂S and both radical species may partly affect the S isotope signatures of pyrite that precipitates from hydrothermal fluids of different composition and sulfur speciation. If $S_3^-$ and its complexes with Fe act as a precursor of pyrite formation at hydrothermal conditions (e.g., Pokrovski et al., 2021b), $S_3^+$ might transmit its lighter isotope signature to pyrite precipitating from S-rich fluids containing both sulfide and sulfate, in which radical ions are more abundant. In contrast, pyrite precipitating from more reduced, acidic and H₂S-dominated fluids, in which the radicals are have much lower concentrations, would be expected to reflect the isotopically heavier H₂S signatures, at least at earlier, out-of-equilibrium, stages of precipitation. This idealized scenario may, however, be complicated by other factors potentially affecting aqueous sulfur-pyrite isotope fractionation such as precipitation/equilibration rates, solution supersaturation degree, fractionation during growth of different crystallographic phases (e.g., Kokh et al., 2020 and references therein). Notwithstanding those factors, our new data and future MD-based calculations of isotopic properties of sulfur species may help tuning sulfur isotope fractionation models, thereby contributing to better understanding the complex framework of fluid-mineral interactions in hydrothermal systems.

More generally, our novel molecular dynamics-based approach for predicting mass-dependent sulfur isotope fractionation in hydrothermal fluids opens large perspectives for the use of MD approaches for assessing isotopic properties of other elements in geological fluids. In fact, the rapidly growing body of MD simulations of different complexes of economic metals with water, chloride and sulfur ligands in hydrothermal fluids relevant to ore formation has been until present primarily focused on determining metal species structures and stabilities to help interpreting spectroscopic and solubility data (e.g., Liu et al., 2013; Mei et al., 2013; Pokrovski et al., 2015, 2022; Laskar et al., 2022), to name a few). The same MD-generated trajectories in such calculations may therefore be used in parallel to assess metal and ligand isotopic properties as well, thereby increasing the impact of the MD approaches in geosciences.

### 5. Conclusions

In line with previous studies that applied the atomic kinetic energy method (“VDOS” method) to iron isotope fractionations (Pinilla et al., 2021; Rabin et al., 2023), here we have tested this new modelling approach on a study of aqueous complexes, such as sulfur species in hydrothermal fluids. Our work assessed the advantages and limitations of this promising method compared to a more commonly used method based on snapshot sampling from the molecular dynamics trajectory (“Snapshot” method).

The major advantage of the “Snapshot” method is that it relies only on structural configurations. It is therefore robust and can be satisfied with relatively short MD trajectories. On the other hand, this method requires additional calculations of structural relaxation and vibrational frequencies on each of the snapshots. These calculations, in particular the calculation of phonons, can be very demanding in terms of computing resources, depending on the size of the system. Moreover, the structural relaxation of the snapshots results in the partial loss of the effects of temperature and configurational disorder characterizing the system.

The “VDOS” method has the great advantage of being able to predict the equilibrium isotope fractionation factor (β-factor) directly from the molecular dynamics trajectory. Thus, the information contained in the configurational disorder (i.e. structural variability and part of the anharmonicity) is taken into account. However, this method requires an accurate molecular dynamics trajectory with a perfect thermal
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** A. Sulfur β-factors of the four aqueous species calculated at three temperatures using Eq. 2 and reported as a function of sulfur interatomic force constant obtained from Eq. 9 ("VDOS" method). The linear fit for a given temperature enables to estimate the β-factor corresponding to any force constant. B. Applying the linear fits of panel A, sulfur β-factor corresponding to the 〈F<sub>SNAP</sub>〉 of 280 N/m for H<sub>2</sub>S is shown as a function of temperature, together with the corresponding polynomial fit.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Temperature dependence of equilibrium sulfur isotope fractionation factors between SO<sub>4</sub><sup>2-</sup> and H<sub>2</sub>S (Panel A), between S<sub>3</sub><sup>•</sup> and H<sub>2</sub>S, and between S<sub>2</sub><sup>•</sup> and H<sub>2</sub>S (Panel B, in blue and red respectively). Fractionation factors are calculated using the H<sub>2</sub>S β-factor obtained by the "VDOS" or "Snapshot" method (solid and dashed lines, respectively). The corresponding polynomial fits are given in Table 5. These present study results are compared with previous theoretical (Otake et al., 2008; Tossell, 2012; Eldridge et al., 2021) and experimental (Ohmoto and Lasaga, 1982; Syverson et al., 2015; Kokh et al., 2020) data. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

equilibration for all elements of the system. Compared to the "Snapshot" method, the "VDOS" method yields β-factor values with larger uncertainties, which decrease only slightly with increasing the duration of the molecular dynamics simulations.

The "VDOS" method was applied to four aqueous sulfur species

(SO<sub>4</sub><sup>2-</sup>, H<sub>2</sub>S, S<sub>3</sub><sup>•</sup>, S<sub>2</sub><sup>•</sup>) under hydrothermal conditions. For the sulfate-sulfide fractionation (SO<sub>4</sub><sup>2-</sup> - H<sub>2</sub>S), our results are slightly lower than the majority of available experimental and theoretical studies. For the radical species S<sub>2</sub><sup>2-</sup> and S<sub>3</sub><sup>•</sup>, the equilibrium isotopic fractionation between them is small, i.e. <1‰ over the temperature range 200-500 °C,
<caption>**Table 5**</caption>

Polynomial fits of equilibrium S isotope fractionation factors as a function of temperature (473–773 K range) calculated relative to the H₂S β-factor obtained by the “VDOS” or “Snapshot” method: 1000 ln α = 1000 ln β(A) – 1000 ln β(B) = αx + bx² with x = 10⁶/T² (T in K).

<table><thead><tr><th>A-B pair</th><th>α</th><th>b</th></tr></thead><tbody><tr><td>SO<sub>3</sub><sup>-</sup> - H<sub>2</sub>SVDOS</td><td>6.1915</td><td>-0.09834</td></tr><tr><td>SO<sub>3</sub><sup>-</sup> - H<sub>2</sub>SSNAP</td><td>5.3519</td><td>-0.12004</td></tr><tr><td>S<sub>3</sub><sup>-</sup> - H<sub>2</sub>SVDOS</td><td>-0.47822</td><td>0.07467</td></tr><tr><td>S<sub>3</sub><sup>-</sup> - H<sub>2</sub>SSNAP</td><td>-1.31790</td><td>0.05297</td></tr><tr><td>S<sub>2</sub><sup>-</sup> - H<sub>2</sub>SVDOS</td><td>-0.71093</td><td>0.07615</td></tr><tr><td>S<sub>2</sub><sup>-</sup> - H<sub>2</sub>SSNAP</td><td>-1.55060</td><td>0.05445</td></tr></tbody></table>

with a slight enrichment in heavier isotope for S₃⁻. Both radical species are isotopically lighter than H₂S by a similar value.

### CRediT authorship contribution statement

**Marc Blanchard:** Conceptualization, Formal analysis, Writing – original draft, Writing – review & editing. **Elsa Desmaele:** Formal analysis. **Gleb S. Pokrovski:** Funding acquisition, Writing – review & editing. **Carlos Pinilla:** Conceptualization, Writing – review & editing. **Merlin Méheut:** Conceptualization, Writing – review & editing. **Rodolphe Vuilleumier:** Conceptualization, Supervision, Writing – review & editing.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.

### Acknowledgements

This research was funded by the Agence Nationale de la Recherche (grant RadicalS - ANR-16-31CE-0017) and CNRS through the MITI interdisciplinary programs (grant ExtremeS, Conditions Extremes 2024). Calculations were performed using the HPC resources from CALMIP (grant 2022 - P1037) and GENCI (CINES/IDRIS, Grant No. 2022- A0130802309). We thank D. Eldridge and an anonymous reviewer for their comments and M. Fiorentini for editorial handling.

### Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.chemgeo.2024.122202.

### References

Barré, G., Truche, L., Bazarkina, E.F., Michels, R., Dubessy, J., 2017. First evidence of the trisulfur radical ion S₃⁻ and other sulfur polymers in natural fluid inclusions. Chem. Geol. 462, 1–14.

Becke, A., 1988. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 38, 3098–3100.

Berndt, M.E., Buttram, T., Earley III, D., Seyfried Jr., W.E., 1994. The stability of gold polysulfide complexes in aqueous sulfide solutions: 100 to 150°C and 100 bars. Geochim. Cosmochim. Acta 58, 587–594.

Bigeleisen, J., Mayer, M.G., 1947. Calculation of equilibrium constants for isotopic exchange reactions. J. Chem. Phys. 15, 261–267.

Blanchard, M., Balan, E., Schauble, E.A., 2017. Equilibrium fractionation of non-traditional isotopes: a molecular modeling perspective. Rev. Mineral. Geochem. 82, 27–63.

Bussi, G., Donadio, D., Parrinello, M., 2007. Canonical sampling through velocity rescaling. J. Chem. Phys. 126, 014101.

Canfield, D.E., 2001. Biogeochemistry of sulfur isotopes. Rev. Mineral. Geochem. 43 (1), 607–636.

Colin, A., Schmidt, C., Pokrovski, G.S., Wilke, M., Borisova, A.Y., Toplis, M.J., 2020. In situ determination of sulfur speciation and partitioning in aqueous fluid-silicate melt systems. Geochim. Persp. Lett. 14, 31–35.

Ducher, M., Blanchard, M., Balan, E., 2016. Equilibrium zinc isotope fractionation in Zn-bearing minerals from first-principles calculations. Chem. Geol. 443, 87–96.

Ducher, M., Blanchard, M., Balan, E., 2018. Equilibrium isotopic fractionation between aqueous Zn and minerals from first-principles calculations. Chem. Geol. 483, 342–350.

Dupuis, R., Benoit, M., Nardin, E., Méheut, M., 2015. Fractionation of silicon isotopes in liquids: the importance of configurational disorder. Chem. Geol. 396, 239–254.

Eldridge, D.L., Guo, W., Farquhar, J., 2016. Theoretical estimates of equilibrium sulfur isotope effects in aqueous sulfur systems: Highlighting the role of isomers in the sulfite and sulfoxylate systems. Geochim. Cosmochim. Acta 195, 171–200.

Eldridge, D.L., Kamышny Jr., A., Farquhar, J., 2021. Theoretical estimates of equilibrium sulfur isotope effects among aqueous polysulfur and associated compounds with applications to authigenic pyrite formation and hydrothermal disproportionation reactions. Geochim. Cosmochim. Acta 310, 281–319.

Goedecker, S., Teter, M., Hutter, J., 1996. Separable dual-space Gaussian pseudopotentials. Phys. Rev. B 54, 1703–1710.

Grimme, S., Antony, J., Ehrlich, S., Krieg, H., 2010. A consistent and accurate ab initio parameterization of density functional dispersion correction (DFTD) for the 94 elements H-Pu. J. Chem. Phys. 132, 154104/1–9.

Herzberg, G., 1966. Electronic Spectra and Electronic Structure of Polyatomic Molecules. Van Nostrand, New York, p. 745.

Hu, M., Chou, I.M., Wang, R., Shang, L., Chen, C., 2022. High solubility of gold in H₂S-H₂O ± NaCl fluids at 100-200 MPa and 600-800 °C: a synthetic fluid inclusion study. Geochim. Cosmochim. Acta 330, 116–130.

Jacquetem, N., Guillaume, D., Zwick, A., Pokrovski, G.S., 2014. In-situ Raman spectroscopy identification of the S₃⁻ ion in S-rich hydrothermal fluids from synthetic fluid inclusions. Am. Mineral. 99, 1109–1118.

Johnston, D.T., 2011. Multiple sulfur isotopes and the evolution of Earth's surface sulfur cycle. Earth Sci. Rev. 106, 161–183.

Kokh, M.A., Assayag, N., Mounic, S., Cartigny, P., Gurenko, A., Pokrovski, G.S., 2020. Multiple sulfur isotope fractionation in hydrothermal systems in the presence of radical ions and molecular sulfur. Geochim. Cosmochim. Acta 285, 100–128.

Kowalski, P.M., Jahn, S., 2011. Prediction of equilibrium Li isotope fractionation between minerals and aqueous solutions at high P and T: an efficient ab initio approach. Geochim. Cosmochim. Acta 75, 6112–6123.

Kühne, T.D., Iannuzzi, M., Del Ben, M., Rybkin, V.V., Seewald, P., Stein, F., Laino, T., Khaliullin, R.Z., Schütt, O., Schiffmann, F., Golze, D., Wilhelm, J., Chulkov, S., Bani-Hashemian, M.H., Weber, V., Borstnik, U., Taillefumier, M., Jakobovits, A.S., Lazzaro, A., Pabst, H., Müller, T., Schade, R., Guidon, M., Andermatt, S., Holmberg, N., Schenter, G.K., Hehn, A., Bussy, A., Belleflamme, F., Tabacchi, G., Glöß, A., Lass, M., Bethune, I., Mundy, C.J., Plessl, C., Watkins, M., VandeVondele, J., Krack, M., Hutter, J., 2020. CP2K: an electronic structure and molecular dynamics software package - Quickstep: Efficient and accurate electronic structure calculations. J. Chem. Phys. 152 (19), 194103.

Laskar, C., Bazarkina, E.F., Kokh, M.A., Hazemann, J.-L., Vuilleumier, R., Desmaele, E., Pokrovski, G.S., 2022. Stability and structure of platinum sulfide complexes in hydrothermal fluids. Geochim. Cosmochim. Acta 336, 407–422.

Lee, C., Yang, W., Parr, R., 1988. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron-density. Phys. Rev. B 37, 785–789.

Liu, X., Sprik, M., Cheng, J., 2013. Hydration, acidity and metal complexing of polysulfide species: a first principles molecular dynamics study. Chem. Phys. Lett. 563, 9–14.

McGinnety, J.A., 1972. Redetermination of the structure of potassium sulfate and potassium chromate: the effect of electrostatic crystal force upon observed bond lengths. Acta Crystallogr. Sect. B: Struct. Sci. 28, 2845–2852.

Méheut, M., Lazzeri, M., Balan, E., Mauri, F., 2007. Equilibrium isotopic fractionation in the kaolinite, quartz, water system: predictions from first-principles density functional theory. Geochim. Cosmochim. Acta 71, 3170–3181.

Mei, Y., Sherman, D.M., Liu, W., Brugger, J., 2013. Complexation of gold in S-3-rich hydrothermal fluids: evidence from ab-initio molecular dynamics simulations. Chem. Geol. 347, 34–42.

Ohmoto, H., Lasaga, A.C., 1982. Kinetics of reactions between aqueous sulfates and sulfides in hydrothermal systems. Geochim. Cosmochim. Acta 46, 1727–1745.

Ohmoto, H., Rye, R.O., 1979. Isotopes of sulfur and carbon. In: Barnes, H.L. (Ed.), Geochemistry of Hydrothermal Ore Deposits, second ed. John Wiley & Sons, pp. 509–567.

Otake, T., Lasaga, C.A., Ohmoto, H., 2008. Ab initio calculations for equilibrium fractionations in multiple sulfur isotope systems. Chem. Geol. 249, 357–376.

Pinilla, C., Blanchard, M., Balan, E., Natarajan, S.K., Vuilleumier, R., Mauri, F., 2015. Equilibrium magnesium isotope fractionation between aqueous Mg²⁺ and carbonate minerals: Insights from path integral molecular dynamics. Geochim. Cosmochim. Acta 163, 126–139.

Pinilla, C., de Moya, A., Rabin, S., Morard, G., Roskosz, M., Blanchard, M., 2021. First-principles investigation of equilibrium iron isotope fractionation in Fe₁ₓSₓ alloys at Earth's core formation conditions. Earth Planet. Sci. Lett. 569, 117059.

Pokrovski, G.S., Dubessy, J., 2015. Stability and abundance of the trisulfur radical ion in hydrothermal fluids. Earth Planet. Sci. Lett. 411, 298–309.

Pokrovski, G.S., Dubrovinsky, L.S., 2011. The S₃⁻ ion is stable in geological fluids at elevated temperatures and pressures. Science 331, 1052–1054.

Pokrovski, G.S., Tagirov, B.R., Schott, J., Hazemann, J.-L., Proux, O., 2009. A new view on gold speciation in sulfur-bearing hydrothermal fluids from in situ X-ray absorption spectroscopy and quantum-chemical modeling. Geochim. Cosmochim. Acta 73, 5406–5427.
Pokrovski, G.S., Akinfiev, N.N., Borisova, A.Y., Zotov, A.V., Kouzmanov, K., 2014. Gold speciation and transport in geological fluids: Insights from experiments and physical-chemical modeling. In: Garofalo, P., Ripley, E. (Eds.), Gold-Transporting Fluids in the Earth's Crust, Geol Soc London Spec Publ, vol. 402, pp. 9–70. https://doi.org/10.1144/sp402.4.

Pokrovski, G.S., Kokh, M.A., Guillaume, D., Borisova, A.Y., Gisquet, P., Hazemann, J.-L., Lahera, E., Del Net, W., Proux, O., Testemale, D., Haigis, V., Jonchière, R., Seitsonen, A.P., Ferlat, G., Vuilleumier, R., Saitta, A.M., Boiron, M.-C., Dubessy, J., 2015. Sulfur radical species form gold deposits on Earth. Proc. Natl. Acad. Sci. 112 (44), 13484–13489.

Pokrovski, G.S., Kokh, M.A., Desmaele, E., Laskar, C., Bazarkina, E.F., Borisova, A.Y., Testemale, D., Hazemann, J.-L., Vuilleumier, R., Ferlat, G., Saitta, A.M., 2021a. The trisulfur radical ion S<sub>3</sub><sup>2-</sup> controls platinum transport by hydrothermal fluids. Proc. Natl. Acad. Sci. 118(34), e2109768118. doi: 10.1073/pnas.2109768118.

Pokrovski, G.S., Blanchard, M., Poitrasson, F., Saunier, G., 2021b. Mechanisms and rates of pyrite precipitation from hydrothermal fluids revealed by iron isotopes. Geochim. Cosmochim. Acta 304, 281–304.

Pokrovski, G.S., Desmaele, E., Laskar, C., Bazarkina, E.F., Testemale, D., Hazemann, J.-L., Vuilleumier, R., Seitsonen, A.P., Ferlat, G., Saitta, A.M., 2022. Gold speciation in hydrothermal fluids revealed by in situ high energy resolution X-ray absorption spectroscopy. Am. Mineral. 107, 369–376.

Polyakov, V.B., 1991. Quantum-statistical discussion of the isotopic bond number method (estimation of isotopic partition functions). Russ. J. Phys. Chem. 65, 694–698.

Rabin, S., Blanchard, M., Pinilla, C., Poitrasson, F., Grégoire, M., 2023. Iron and silicon isotope fractionation in silicate melts using first-principles molecular dynamics. Geochim. Cosmochim. Acta 343, 212–233.

Richet, P., Bottinga, Y., Javoy, M., 1977. A review of hydrogen, carbon, nitrogen, oxygen, and chlorine stable isotope fractionation among gaseous molecules. Annu. Rev. Earth Planet. Sci. 5, 65–110.

Schmidt, C., Seward, T.M., 2017. Raman spectroscopic quantification of sulfur species in aqueous fluids: Ratios of relative molar scattering factors of Raman bands of H₂S, HSO₃⁻, SO₂, HSO₄⁻, SO₄²⁻, S₂O₃²⁻, S₃ and H₂O at ambient conditions and information on changes with pressure and temperature. Chem. Geol. 467, 64–75.

Seward, T.M., 1973. Thio complexes of gold and the transport of gold in hydrothermal ore solutions. Geochim. Cosmochim. Acta 37, 379–399.

Syverson, D.D., Ono, S., Shanks, W.C., Seyfried Jr., W.E., 2015. Multiple sulfur isotope fractionation and mass transfer processes during pyrite precipitation and recrystallization: an experimental study at 300 and 350°C. Geochim. Cosmochim. Acta 165, 418–434.

Thode, H.G., Macnamara, J., Collins, C.B., 1949. Natural variations in the isotopic content of Sulphur and their significance. Can. J. Res. B27, 361–373.

Thode, H.G., Monster, J., Dunford, H.B., 1961. Sulphur isotope geochemistry. Geochim. Cosmochim. Acta 25, 159–174.

Tossell, J.A., 2012. Calculation of the properties of the S<sub>3</sub><sup>2-</sup> radical anion and its complexes with Cu<sup>2+</sup> in aqueous solution. Geochim. Cosmochim. Acta 95, 79–92.

Truche, L., Bazarkina, E.F., Barré, G., Thomassot, E., Berger, G., Dubessy, J., Robert, P., 2014. The role of S<sub>3</sub><sup>2-</sup> ion in thermochemical sulfate reduction: Geological and geochemical implications. Earth Planet. Sci. Lett. 396, 190–200.

VandeVondele, J., Hutter, J., 2007. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. J. Chem. Phys. 127, 114105.

VandeVondele, J., Krack, M., Mohamed, F., Parrinello, M., Chassaing, T., Hutter, J., 2005. Quickstep: fast and accurate density functional calculations using a mixed gaussian and plane waves approach. Comput. Phys. Commun. 167, 103–128.

Walters, J.B., Cruz-Uribe, A.M., Marschall, H.R., 2020. Sulfur loss from subducted altered oceanic crust and implications for mantle oxidation. Geochem. Persp. Let. 13, 36–41.