Theoretical estimates of sulfoxyanion triple-oxygen equilibrium isotope effects and their implications

Triple-oxygen isotope ( d 18 O and D


Introduction
Sulfur can exist in a range of redox states from S(-II) to S(+VI).As such, its oxidation and reduction represent major electron fluxes into and out of Earth's biosphere; these fluxes regulate atmospheric O 2 content and the Earth-surface redox state over multi-million year timescales (Berner, 2001).Today, dissolved marine sulfate [S(+VI) redox state] constitutes one of the largest oxidant reservoirs on Earth's surface (Blättler et al., 2018), whereas sulfide minerals such as pyrite [FeS 2 ; S(-I)] contained in marine sediments and sedimentary rocks constitute one of the largest reductant reservoirs (Berner, 1984).
Several biotic and abiotic processes can transfer sulfur between these oxidized and reduced forms, often via sulfoxyanion species that exist at intermediate redox states.For example, microbial sulfate reduction (MSR) is a metabolism that gains energy in the absence of O 2 by reducing dissolved sulfate to hydrogen sulfide [S(-II)]-the precursor to pyrite-via the intermediate species sulfite [S(+IV)] (Fike et al., 2015).Related metabolisms gain energy by simultaneously oxidizing and reducing sulfur to sulfate and sulfide via the disproportionation of intermediate redox compounds such as sulfite, zero-valent sulfur [S(0)], or the mixed-valence species thiosulfate [S(-I)/S(+V); see Vairavamurthy et al. (1993) for atom-specific redox state determination] (Jørgensen, 1990;Fike et al., 2015).Similarly, the (a)biotic oxidative weathering of pyrite in exhumed rocks must occur via sulfite and thiosulfate intermediates, although the exact mechanism is complex and not fully constrained (e.g., Balci et al., 2007;Schoonen et al., 2010;Kohl and Bao, 2011).Furthermore, sulfoxyanion species existing in intermediate redox states are generally short-lived at Earth-surface conditions.For example, the hydrolysis and oxidation of sulfur dioxide gas [S(+IV)] to sulfate in atmospheric water-which represents a major pathway of acid rain formation-occurs on the order of microseconds (Brandt and van Eldik, 1995).The residence time of intracellular sulfite produced during MSR is similarly estimated to be on the order of microseconds before it is either fully reduced to hydrogen sulfide or reoxidized to sulfate (Bertran et al., 2020).Thus, even though the overall abundance of intermediate sulfoxyanions on Earth's surface at any point in time is low, nearly all sulfur-cycle processes require their transient production and consumption (e.g., Jørgensen, 1990).
One method to assess the relative importance of various sulfurcycle processes is by analyzing the sulfur and oxygen isotope compositions of sulfate ( 33 S/ 32 S, 34 S/ 32 S, 36 S/ 32 S, 17 O/ 16 O, and 18 O/ 16 O; reported as D 33 S, d 34 S, D 36 S, D 0 17 O, and d 18 O, respectively; see Section 2.1).For example, Wing and Halevy (2014) and Bertran et al. (2020) showed that 34 S and 18 O fractionation during MSR represents an intermediate between kinetic and equilibrium fractionation factors; MSR rates in marine sediments can thus be predicted using the d 34 S and d 18 O values of residual sulfate.Similarly, d 34 S, d 18 O, and D 0 17 O values of dissolved sulfate in rivers has been used to estimate the relative importance of anaerobic pyrite weathering, aerobic pyrite weathering, and evaporite dissolution on land (e.g., Turchyn et al., 2013;Burke et al., 2018;Killingsworth et al., 2018;Hemingway et al., 2020;Burt et al., 2021;Kemeny et al., 2021).Furthermore, sulfate is preserved in minerals such as anhydrite/gypsum (CaSO 4 ), barite (BaSO 4 ), and carbonate (as carbonate-associated sulfate, or CAS); the sulfur and oxygen isotope compositions of these minerals can thus be used to reconstruct sulfur-cycle processes through time.
Despite the utility of sulfate isotopes as geologic tracers, their proper interpretation requires knowledge of fractionation factors for each step of each sulfur-cycle process.Although sulfurisotope fractionation has been studied extensively (e.g., Eldridge et al., 2016;Eldridge et al., 2021), less is known about sulfoxyanion oxygen-isotope fractionation.Specific to d 18 O and D 0 17 O, intermediate sulfoxyanions can rapidly exchange oxygen atoms with surrounding water, potentially exhibiting redox state-and speciesor isomer-specific equilibrium fractionation factors (e.g., Pryor and Tonellato, 1967;Betts and Voss, 1970;Müller et al., 2013a;Wankel et al., 2014).Thus, the final oxygen isotope composition of sulfate produced or consumed by any process should depend strongly on (i) the isotope composition of water in which it formed and (ii) the specific intermediate sulfoxyanion species involved.Although some empirical 18 O fractionation estimates exist for sulfate, sulfite, and thiosulfate (Lloyd, 1968;Chiba et al., 1981;Müller et al., 2013a;Wankel et al., 2014, this study), many sulfoxyanion oxygen-isotope fractionation factors-especially those for 17 O-remain unknown due to the difficulty of experimentally measuring these short-lived compounds.This lack of fractionation factor constraints hinders our ability to interpret sulfate d 18 O and D 0 17 O values, both today and in the geologic past.
In the absence of experimental constraints, quantum-chemistry computational methods have been shown to yield acceptable triple-oxygen isotope fractionation factor estimates for a range of oxygen-containing anions and minerals (e.g., Cao and Liu, 2011;Hayles et al., 2018;Schauble and Young, 2021;Yeung and Hayles, 2021).Importantly, D 0 17 O estimates are particularly robust since any biases (e.g., arising from methodological inadequacies in potential energy surface calculations) largely cancel due to the mass-dependent nature of equilibrium 17 O fractionation (Cao and Liu, 2011).Thus, if computational 18 O fractionation factors can be shown to reasonably match experimental constraints, then D 0 17 O is expected to be accurate within analytical uncertainty.
The remainder of this study is organized as follows: (i) first, we outline the necessary notation and quantum mechanical theory to computationally estimate triple-oxygen isotope fractionation factors (Section 2); (ii) next, describe the computational methods used, including those to estimate methodological scaling factors and the importance of anharmonic zero-point energy (ZPE) corrections (Section 3); (iii) we then report predicted fractionation factors and compare to available experimental results from the literature, including new experimental results included as part of this study (Section 4 and Supplementary Material); and finally (iv) we interpret these fractionation factors within the context of several sulfur-cycle processes and compare predictions to environmental data (Section 5).This work-combined with Eldridge et al. (2016) and Eldridge et al. (2021)-yields equilibrium fractionation factor estimates of all major and minor isotopes ( 33 S/ 32 S, 34 S/ 32 S, 36 S/ 32 S, 17 O/ 16 O, and 18 O/ 16 O) for several important sulfoxyanion species.

Notation
The oxygen-isotope composition of a given compound ''A" can be written as where Ã R denotes the ⁄ O/ 16 O ratio, ''*" indicates the heavy isotope 17 O or 18 O, and VSMOW is the Vienna Standard Mean Ocean Water international reference standard.Here we report results in units of ''permil" (‰) by multiplying Eq. 1 by 1000.For any fractionation process (kinetic or equilibrium), the oxygen-isotope composition of product ''A" and reactant ''B" are related by the fractionation factor: While not explicitly stated here for notational convenience, all equilibrium fractionation is temperature dependent.When considering all three oxygen isotopes, 18 a A=B and 17 a A=B are related by the mass law for a given process: Although each process is described by a unique mass law, fractionation is deemed ''mass dependent" if where 17 h RL is the mass law for the ''reference line" and prime (0) indicates the use of logarithmic d Ã O values.Like d Ã O A , here we report D 0 17 O A in units of permil by multiplying Eq. 4 by 1000.
Although the choice of reference line is arbitrary, here we set 17 h RL ¼ 0:5305 since this corresponds to the high-temperature limit of equilibrium fractionation; i.e., the ratio of reduced masses between isotopes of atomic oxygen (Young et al., 2002).Finally, the temperature-dependent D 0 17 O offset between product ''A" and reactant ''B" for any fractionation process is defined as Our goal is to theoretically estimate 18 a A=H 2 Oðliq:Þ and DD 0 17 O A=H 2 Oðliq:Þ over a range of environmentally relevant temperatures, where ''A" is any sulfoxyanion species of interest.

The Bigeleisen-Goeppert Mayer-Urey equation
To estimate fractionation factors, we utilize the Bigeleisen-Goeppert Mayer-Urey (B-GM-U) equation, which predicts the equilibrium constant of isotope substitution using translational, rotational, and vibrational reduced partition function ratios (RPFRs) for each isotopologue of a given molecule (Bigeleisen and Goeppert, 1947;Urey, 1947).It is subject to four main approximations: (i) a molecule's rotational motion can be treated as a rigid rotor, (ii) its vibrational motion can be treated as a harmonic oscillator, (iii) the motion of electrons and nuclei are decoupled such that isotopic substitution has no effect on electronic potential energy surface and molecular structure (the so-called Born-Oppenheimer approximation; Born and Oppenheimer, 1927), and (iv) the ratio of the moments of inertia of two isotopically substituted molecules depends only on their masses and the product of their vibrational frequency ratios (the so-called Teller-Redlich product rule; Redlich, 1935;Wilson et al., 1955).3 , ðHSÞSO À 3 , and S2O2ðOHÞ À .Reported bond lengths and angles refer to those calculated for a 30 Á H2O water droplet cluster using the B3LYP/6-31+G(d,p) method (-OH bond lengths and angles are not shown but are included in Table S3); water molecules are omitted for visual clarity.Atom colors refer to: yellow = sulfur, red = oxygen, white = hydrogen.
By following these approximations and utilizing a statistical mechanical treatment of molecular motion the B-GM-U equation can be written as: h is Planck's constant, c is the speed of light, x i is the harmonic normal mode frequency for degree-of-freedom i; k B is Boltzmann's constant, T is temperature in Kelvin, n is the number of atoms in the compound of interest, and x ¼ 5 for linear molecules or x ¼ 6 for nonlinear molecules.As above, ''*" denotes terms related to the compound containing the heavy isotope, whereas the subscript ''h" indicates the pure harmonic approximation result.We have amended Eq. 6 with the subscripts ''TR", ''ZPE", and ''EXC" to denote contributions from partition functions for translation/rotation, zero-point energy of vibration, and excited vibrational states, respectively.
To tabulate fractionation factors for molecules containing several non-equivalent isotopically substitutable atoms (as is the case for oxygen atoms in sulfoxyanions considered here), Richet et al. (1977) defined the Ã b factor, which represents bulk isotope fractionation between the molecule of interest and an idealized mono-atomic, non-interacting gas (e.g., an O atom).Ignoring multiply substituted isotopologues, Ã b can be approximated as the geometric mean of the RPFRs for all singly-substituted isotopomers of a given molecule: where m is the number of atoms of the element of interest (e.g., oxygen) contained within the molecule (e.g., Richet et al., 1977;Liu et al., 2010).Eq. 8 states that Ã b % Ã RPFR only for the special case where all isotopically substitutable atoms are equivalent (e.g., oxygen-isotope substitution in symmetric sulfoxy gases; Fig. S1).
Oxygen-18 fractionation between two compounds ''A" and ''B" can then be written as the ratio of their 18 b values: Following Cao and Liu (2011), we extend this nomenclature for all three oxygen isotopes by defining: which similarly represents the equilibrium mass law between the compound of interest and an idealized, mono-atomic, noninteracting O atom.For any compound of interest, it follows that lim T!1 where 16 M, 17 M, and 18 M are the atomic masses of 16 O, 17 O, and 18 O (see Matsuhisa et al., 1978;Young et al., 2002, for derivation).That is, 17 j should approach the high-temperature limit of equilibrium fractionation.Combining Eqs. 3, 9, and 10, the mass law for equilibrium fractionation between two compounds ''A" and ''B" can be written as a function of their 18 b and 17 j values: Similar to previous studies (e.g., Cao and Liu, 2011;Hayles et al., 2018;Schauble and Young, 2021), we report all results as 18 b and 17 j rather than reporting 18 a A=B and 17 h A=B directly.

Corrections to the B-GM-U equation
The approximations required to derive the B-GM-U equation can lead to large deviations between predicted and experimental equilibrium fractionation factors; several studies have thus proposed RPFR correction terms to reduce this inaccuracy (e.g., Richet et al., 1977;Liu et al., 2010;Zhang and Liu, 2018;Schauble and Young, 2021).Theoretical corrections generally fall into one of three categories: (i) deviations from the rigid rotor approximation, including quantum mechanical rotation corrections, centrifugal distortion, and torsion effects (Liu et al., 2010); (ii) deviations from the harmonic oscillator approximation, including anharmonic vibrational energy-state corrections and doublewell potentials (Liu et al., 2010;Schauble and Young, 2021); and (iii) coupling, including vibration-rotation coupling and electronnuclear coupling (i.e., deviations from the Born-Oppenheimer approximation; Born and Huang, 1954;Liu et al., 2010;Zhang and Liu, 2018).
As has been done previously, we ignore corrections of categories (i) and (iii) throughout this study.Specifically, Liu et al. (2010) have shown using a suite of gaseous compounds that all deviations from the rigid rotor approximation for rotational RPFRs, as well as vibration-rotation coupling, negligibly impact Ã b estimates for non-hydrogen elements, including oxygen.Similarly, while inaccuracies due to deviations from the Born-Oppenheimer approximation can become significant at ultra-cold temperatures and for hydrogen-isotope fractionation, their impacts on 18 b and 17 b are likely small at temperatures relevant to Earth-surface conditions.For example, Zhang and Liu (2018) predict that the Born-Oppenheimer approximation leads to an over-estimate of 18 a between gas-phase H 2 O and SO 3 of 1.7‰ at À25 °C and that this offset decreases with increasing temperature.While not negligible, this inaccuracy is well within the experimental 18 a uncertainty for many species studied here (cf., Lloyd, 1968;Chiba et al., 1981;Müller et al., 2013a;Wankel et al., 2014, this study).
For category (ii) errors, Liu et al. (2010) predict that anharmonic corrections to the zero-point energy (ZPE) of vibration may exhibit an appreciable impact on Ã RPFR and thus Ã b estimates for nonhydrogen elements, whereas anharmonic corrections to excited vibrational states do not.Furthermore, Schauble and Young (2021) recently concluded that double-well vibrational potentials-which are present in hydrogen-bonding solutions such as water-do not influence D 017 O by more than 0.01‰.We therefore ignore corrections for anharmonic excited states and double-well potentials, but we do estimate anharmonic ZPE effects for gaseous species using a modified B-GM-U equation.Specifically, if the anharmonic ZPE can be quantified directly, then the ZPE partition functions can be removed from the RPFR product and Eq. 6 can be rewritten as f is the total ZPE (including harmonic and anharmonic contributions), and the subscript ''AnZPE" indicates that the result includes anharmonic ZPE corrections.Combining Eq. 6 and 13 yields where is the anharmonic correction to the ZPE partition function ratio following the nomenclature of Liu et al. (2010).
In practice, direct determination of f values for aqueous sulfoxyanions is not feasible due to computational constraints, whereas empirical ZPE correction factors likely carry prohibitively large bias and uncertainty (e.g., Irikura et al., 2009).Following current best practice (e.g., Liu et al., 2010;Eldridge et al., 2016;Liu and Liu, 2016), we therefore instead qualitatively assess the impact of anharmonicity in Section 4.1 using ZPEs determined for ideal gaseous (''in vacuo") sulfoxy species with and without anharmonic corrections.We specifically calculate Ã AnZPE for a suite of sulfoxy gases (Fig. S1) using Eq.16 and discuss the likely magnitude of anharmonic ZPE corrections on aqueous sulfoxyanions in light of these results.

Potential energy surface corrections
Computational constraints require that methods with a less accurate potential energy surface (PES) are used when analyzing large molecules-including explicitly solvated anions-resulting in x values that are systematically biased (Irikura et al., 2005).These biases can be partially corrected by empirically scaling sulfoxyanion x values to those determined using highly accurate PES methods (e.g., Scott and Radom, 1996;Merrick et al., 2007). 1   Here, we calculate PES scaling factors as the ratio of x values for gaseous sulfoxy species determined using several computational theories and basis sets; we report all scaled results with the subscript ''Sc".An analogous approach was used in a previous theoretical estimate of 18 O fractionation between SO 2À  4 and water; it was shown to reduce the misfit between theoretical experimental results by several permil, highlighting the utility of PES scaling factors when determining 18 b values (Zeebe, 2010). 2   Finally, because x and ZPE values are related to electron density, it is reasonable to hypothesize that correction factors depend on sulfur-atom(s) redox state(s).This hypothesis is supported by Eldridge et al. (2016), who observed a larger anharmonic ZPE correction for in vacuo species in the S(+VI) redox state relative to more reduced species.As detailed in Sections 3 and 4, we therefore explore the effect of redox-specific PES scaling factors on resulting Ã b h;Sc values for aqueous sulfoxyanions.

Methods
Descriptions of the experimental methods for all sulfite-water and thiosulfate-water equilibrium experiments included as part of this study are reported in the Supplementary Material.All ZPEs and harmonic frequencies were estimated using the computational chemistry software Gaussian 16 on the Research Computing cluster at Harvard University (Frisch et al., 2016).For all species, geometry optimizations were performed first, followed by isotopespecific frequency calculations.The masses of all atoms other than the oxygen atom of interest were assumed to equal their respective major isotopes (i.e., 1 H, 16 O, and 32 S).For species with several possible isotopomers (e.g., isotopic substitution at each O atom in SO 2À  4 ), frequencies and RPFRs were calculated individually for each isotopomer, and resulting Ã b values were subsequently calculated using Eq. 8.After each step, convergence was ensured by confirming that stationary points were found (i.e., that no imaginary frequencies exist).
All optimization and frequency calculations were performed using one of three methods of varying accuracy: (i) B3LYP/6-31 +G(d,p) (Lee et al., 1988;Becke, 1993), (ii) CCSD/aug-cc-pVTZ (Scuseria et al., 1988;Kendall et al., 1992), and (iii) MP2/aug-cc-pVTZ (Møller and Plesset, 1934;Frisch et al., 1990;Kendall et al., 1992).Method (i) is a low/moderate-complexity hybrid of Hartree-Fock (HF) and density functional theory (DFT) with a double-zeta Pople basis set, including diffuse and polarization functions.It has been used previously to estimate multiple-sulfur isotope fractionation factors of sulfoxyanions and polysulfur compounds (Eldridge et al., 2016;Eldridge et al., 2021); similar methods have additionally been used to estimate triple-oxygen isotope fractionation factors for a range of gaseous, aqueous, and mineral species (Cao and Liu, 2011;Hill et al., 2014;Hayles et al., 2018;Guo and Zhou, 2019).Here, we applied method (i) to all gaseous and aqueous species, including anharmonic ZPE corrections for gaseous species.Method (ii) is a high-complexity coupled cluster HF theory with a correlation-consistent polarized triple-zeta basis set, augmented with diffuse functions.It is a highly accurate but computationally expensive method and was used here for gaseous species to estimate the effect of basis-set accuracy on fractionation factor results (i.e., to calculate PES scaling factors).Method (iii) utilizes second-order Møller-Plesset perturbation theory with the same basis set used in method (ii).Like method (ii), it is accurate yet computationally expensive; however, unlike method (ii), it additionally allows for anharmonic corrections because the analytic second derivatives can be calculated (Frisch et al., 1990;Frisch et al., 2016).This method was used here for gaseous species to estimate the effect of basis-set accuracy as well as anharmonic ZPE corrections on fractionation factor results (i.e., to calculate PES scaling factors and estimate the impact of ZPE anharmonicity), as has been done previously (Liu et al., 2010).
We note that the atom-specific formal oxidation states of ðSÞSO 2ðgÞ are assumed by comparison to aqueous thiosulfate (Vairavamurthy et al., 1993) and have not been determined experimentally.Similarly, the physical-chemical relationship between gas-phase ðSÞSO 2ðgÞ and aqueous-phase thiosulfate species has not been demonstrated directly.Rather, this is inferred based on theoretical results and by analogy to the known relationship between 1 This differs from the common practice of scaling harmonic to fundamental frequencies, which are inappropriate for the B-GM-U equation (see Liu et al., 2010, for details).

Liquid water computations
To obtain 17 j estimates, liquid water harmonic frequencies were calculated following the procedure of Hayles et al. (2018) using the B3LYP/6-31+G(d,p) method.First, water ''droplets" were generated by starting with 6 H 2 O molecules and iteratively performing a geometry optimization followed by the addition of 4 more H 2 O molecules until a final cluster size of 30 H 2 O molecules was achieved.Such iterative optimizations are commonly performed to ensure that water droplet geometry remains stable (e.g., Li et al., 2009;Zeebe, 2010;Li and Liu, 2011;Eldridge et al., 2016).Then, each H 2 O molecule was individually isotopically substituted and harmonic frequencies were determined, yielding n ¼ 30 sets of frequencies.Finally, to minimize the influence of water droplet cluster geometry on resulting Ã b values, the whole procedure was repeated 4 more times with arbitrary molecular geometries for a total of n ¼ 5 water droplets.

Sources of uncertainty and statistical methods
There exist three main potential sources of uncertainty and/or bias: (i) assumptions required to derive the B-GM-U equation (Section 2.2-2.3);(ii) approximations inherent to the particular computational theory and basis set used to model the molecular PES, including scaling factor uncertainty (Section 2.4); and (iii) variability due to water droplet geometry, for example due to the importance of water dimers and trimers at higher temperature (Eldridge et al., 2016;Hayles et al., 2018).To assess points (i) and (ii), we empirically calculate PES scaling factors and estimate possible bias arising from ZPE anharmonicity for a suite of gaseous sulfoxy compounds using several theories and basis sets (e.g., following Liu et al., 2010;Cao and Liu, 2011).To assess point (iii), we determine triplicate water droplet geometries for a subset of sulfoxyanion species [SO 3 (OH) À and S 2 O 2 ðOHÞ À ].We additionally calculate semi-empirical liquid water Ã b values following the method of Hayles et al. ( 2018) since theoretical results may not fully capture liquid water molecular interactions.
Finally, we assess the overall accuracy of our results by comparing predicted 18 a A=H 2 Oðliq:Þ with all existing theoretical and experimental results from the literature, where ''A" is the sulfoxyanion of interest (Lloyd, 1968;Chiba et al., 1981;Zeebe, 2010;Müller et al., 2013a;Wankel et al., 2014, this study).Unfortunately, no experimental equilibrium DD 0 17 O A=H 2 Oðliq:Þ estimates currently exist for sulfoxyanions.However, as mentioned previously, theoretical DD 0 17 O A=H 2 Oðliq:Þ error has been shown to be small-even if 18 a A=H 2 Oðliq:Þ is in error by several permil-due to mass-dependent error cancellation (Cao and Liu, 2012).We therefore assume that theoretical DD 0 17 O A=H 2 Oðliq:Þ predictions are accurate if the corresponding 18 a A=H 2 Oðliq:Þ agrees with experimental results to within several permil.
All regressions (e.g., when calculating PES scaling factors) were performed using orthogonal distance regression, which allows for uncertainty in both x and y variables.Regression results are reported with AE1r uncertainty.Statistical differences between populations were determined using a two-way analysis of variance (ANOVA) and evaluated at the p ¼ 0:05 level.Misfit between theoretical predictions and experimental results was determined as the root mean square error (RMSE); to estimate bulk solution fractionation for each experiment, ''mean" fractionation factors were calculated as the average of all species (including isomers where appropriate) weighted by their pH-specific relative abundance (see Millero et al., 1989;Eldridge et al., 2018, for speciation constants used here).

Results and interpretation
Detailed descriptions of the experimental results for all sulfitewater and thiosulfate-water equilibrium experiments included as part of this study are reported in the Supplementary Material.Predicted bond lengths and angles for all gaseous species using all methods are reported in Table S1, whereas calculated ZPEs and harmonic frequencies are reported in Table S2.Similarly, predicted bond lengths and angles for aqueous sulfoxyanions are reported in Table S3, whereas ZPEs and harmonic frequencies are reported in Table S4.All optimization and frequency results, as well as python scripts used to calculate all Ã RPFR; Ã b; 17 j; 18 a and D 0 17 O values, are included in the Supplementary Data.Finally, following common practice, seventh-order polynomial fits for calculating 18 b h;Sc and and all sulfoxyanions of interest as a function of temperature are reported in Tables 1and 2. 4.1.Gaseous species and scaling factors For gaseous species, all methods result in similar optimized geometries with slight differences in bond lengths and angles.Bond length differences between methods reach a maximum of 0.12Å (for S-S in S 2 O 2ðgÞ ) and average 0.03Å, with MP2/aug-cc-pVTZ consistently predicting the longest bonds and CCSD/aug-cc-pVTZ consistently predicting the shortest bonds.Bond angle differences between methods are similarly small, reaching a maximum of 6.1°(for O-S-S in S 2 O 2ðgÞ ) and averaging 1.3°with no clear trend between methods.
Here, we assess the impact of PES scaling and anharmonic ZPE correction factors and we compare all results to available experimental data.When considering all gaseous species together, x values predicted by the CCSD/aug-cc-pVTZ method are nearly identical to experimental results for all compounds reported in Johnson (2020), with an experimental vs. theoretical regression slope of 1:0062 AE 0:0002 (r 2 ¼ 0:9999; Fig. S2A).Similarly, x predictions using the MP2/aug-cc-pVTZ method are generally nearly identical to experimental values (noting that experimental x val-  ues are derived from observed fundamental frequencies based on a physical model; e.g., Irikura, 2007).However, all diatomic molecules (SO, S 2 , and O 2 ) deviate significantly from this trend, with predicted values consistently lower than experimental results.
Resulting 18 b h values determined using unscaled harmonic frequencies (Eq.6) for all gaseous sulfoxy species, H 2 O ðvap:Þ , and O 2ðgÞ are shown in Fig. S5.For all methods at a given temperature, sulfoxy species 18 b h values generally increase with increasing sulfur redox state.For example, 18 b h calculated at 25 °C using the B3LYP/6-31+G(d,p) method increases from 1:0738 for S 2 O 2ðgÞ to 1:0947 for SO 3ðgÞ ; H 2 O ðvap:Þ and O 2ðgÞ predictions generally lie between those of S(+II) and S(+IV) species.For a given compound, the CCSD/aug-cc-pVTZ method predicts the highest 18 b h at all temperatures, with B3LYP/6-31+G(d,p) and MP2/aug-cc-pVTZ both predicting lower-yet similar-values.For example, 18 b h calculated at 25 °C for SO 3ðgÞ ranges from 1:0947 using B3LYP/6-31+G(d,p) to 1:1065 using CCSD/aug-cc-pVTZ.For all compounds using all methods, 18 b h values approach unity with increasing temperature, as expected (i.e., no fractionation as T ! 1).
PES scaling factor corrections to B3LYP/6-31+G(d,p) results are shown in Fig. S6.When considering the all-compound average scaling factors, 18 b h;Sc values at 25 °C only differ from their unscaled counterparts by a maximum of 3.6‰ using CCSD/aug-cc-pVTZ scaling and 2.4‰ using MP2/aug-cc-pVTZ scaling [both maxima corre-spond to SO 3ðgÞ ].In contrast, redox-state specific PES scaling factors yield 18 b h;Sc corrections at 25 °C that range from 0.7‰ [H 2 O ðvap:Þ ] to 10.4‰ [SO 3ðgÞ ] using CCSD/aug-cc-pVTZ and from À2.0‰ [SO 2ðgÞ ] to 5.9‰ [ðSÞSO 2ðgÞ ] using MP2/aug-cc-pVTZ.Interestingly, SO 2ðgÞ scaling factors using the MP2/aug-cc-pVTZ method are consistently 6 1 for all T; in contrast, all other PES scaling factors are always P 1.Like 18 b h results, all PES scaling factors approach unity with increasing temperature.
Anharmonic ZPE corrections to 18 b h values (Eq.16) calculated using B3LYP/6-31+G(d,p) and MP2/aug-cc-pVTZ methods are shown in Fig. S7.All corrections range from À1.0‰ [O 2ðgÞ ] to À0.3‰ [SO ðgÞ ] for B3LYP/6-31+G(d,p) and from À1.3‰ [O 2ðgÞ ] to 0.1‰ [SO ðgÞ ] for MP2/aug-cc-pVTZ.Interestingly, S(+II) anharmonic ZPE corrections calculated using MP2/aug-cc-pVTZ are consistently P 1, unlike all other results.Similar to 18 b h and PES scaling factor results, all anharmonic ZPE corrections approach unity with increasing temperature.Across all methods tested here, PES x scaling factors yield consistently larger changes in 18 b than do AnZPE corrections; this is especially true for redox-state specific PES results.For example, the overall impact on resulting 18 b values of redox-state specific PES scaling factors (scaled to CCSD/aug-cc-pVTZ) and of AnZPE corrections [using B3LYP/6-31+G(d,p)] for all gaseous sulfoxy species, H 2 O ðvap:Þ , and O 2ðgÞ are shown in Fig. S8.For all sulfoxy species, PES scaling factors yield 18 b corrections that are P 5:8-fold larger than those from ZPE anharmonicity at 25 °C.
Resulting 17 j h values calculated using unscaled harmonic frequencies (Eq.6) for all gaseous sulfoxy species, H 2 O ðvap:Þ , and O 2ðgÞ are shown in Fig. S9.For all methods at all temperatures, H 2 O ðvap:Þ consistently displays the highest 17 j h values, whereas all other compounds cluster at lower values.For example, the B3LYP/6-31+G(d,p) method at 25 °C predicts a 17 j h value for H 2 O ðvap:Þ of 0:5300 with all other compounds ranging from 0:5282 [SO 2ðgÞ ] to 0:5283 [S 2 O 2ðgÞ ].For all compounds and all methods, 17 j h values approach 0:5305 with increasing temperature, as predicted by the high-temperature theoretical limit (Eq.11; Matsuhisa et al., 1978;Young et al., 2002).
Unlike for 18 b h;Sc predictions, 17 j h;Sc values calculated using PES x scaling factors do not greatly deviate from their unscaled counterparts (Fig. S10), as observed previously (Cao and Liu, 2011).For example, 17 j h;Sc offsets using the all-compound average scaling factors at 25 °C range from À5.9 Â 10 À5 However, anharmonic ZPE correction factors also lead to diverging predictions with increasing temperature.That is, 17 j AnZPE does not converge on 0:5305 as T ! 1, as is theoretically predicted (Eq.11; Matsuhisa et al., 1978;Young et al., 2002).As detailed in Cao and Liu (2011), this results from the 17 AnZPE= 18 AnZPE term, which itself does not converge on 0:5305 as T ! 1.
Based on these results, we choose to scale aqueous-phase x values determined using the B3LYP/6-31+G(d,p) method by redox-state specific CCSD/aug-cc-pVTZ PES scaling factors when calculating all fractionation factors.The reason for choosing CCSD/aug-cc-pVTZ rather than MP2/aug-cc-pVTZ scaling factors is twofold: (i) Eq. 9. We then calculate 17 a liq:=vap: as a function of 18 a liq:=vap: and 17 h liq:=vap: using Eq. 3.However, 17 h liq:=vap: has not yet been empirically determined with sufficient precision across the entire temperature range of interest (cf., Barkan and Luz, 2005).We therefore theoretically estimate 17 h liq:=vap: using 18 b h;Sc and 17 j h;Sc values calculated with the B3LYP/6-31+G(d,p) method for H 2 O ðliq:Þ and H 2 O ðvap:Þ following Eq.12. Resulting theoretical 17 h liq:=vap: agrees with available empirical results to within 2.0 Â 10 À4 , suggesting that our predictions are robust (Fig. S13; Barkan and Luz, 2005).
We then determine 17 b se for H 2 O ðliq:Þ as: 17 b se ¼ 18 a liq:=vap: where 17 b h;Sc here refers to the values for H 2 O ðvap:Þ calculated using the B3LYP/6-31+G(d,p) method with redox-state specific CCSD/ aug-cc-pVTZ scaling factors.Finally, 17 j se is determined following Eq.10.To test the influence of salinity on resulting fractionations, we additionally repeat this process using empirically measured 18 a liq:=vap: values for a 4 M NaCl solution from Horita et al. (1995) (their Eq. 10).
Resulting freshwater 18 b se calculated here only deviates from pure theoretical predictions by a maximum of 1.7‰ and from semi-empirical values reported in Hayles et al. (2018) by a maximum of 3.8‰ (Fig. S13A).Similarly, freshwater 17 j se agrees with theoretical results and with semi-empirical results of Hayles et al. (2018) to within 2.5 Â 10 À5 , indicating that isotope effects not captured by the water droplet method are small (Fig. S13B).Furthermore, salinity exhibits a negligible effect on both 18 b se and 17 j se , leading to maximum deviations from freshwater results of only 0.4‰ and <1 Â 10 À6 , respectively.We therefore use freshwater 18 b se and 17 j se when calculating all sulfoxyanion fractionation factors below (Section 4.3-4.6).Nonetheless, seventh-order polynomial fits for calculating theoretical, semi-empirical (freshwater), and semi-empirical (saline) 18 b and 17 j values as a function of temperature are reported in Tables 1,2.

Aqueous sulfate species
On average, sulfate species geometries are similar to those calculated by Eldridge et al. (2016), with bond lengths differing by a maximum of 0.01Å and angles differing by a maximum of 3.1°( Fig. 1; Table S3).Replicate SO 3 (OH) À geometries calculated here show nearly identical S-O bond lengths, differing by only 0.01Å.However, replicate S-(OH) and O-H bond lengths differ by up to 0.10Å and 0.09Å; similarly, O-S-O, O-S-(OH), and S-O-H bond angles differ by up to 4.0°, 4.2°, and 6.4°between replicates, likely due to the influence of H 2 O geometry on hydrogen bond length and angle.Despite these geometric differences, 18 b h;Sc and 17 j h;Sc are nearly identical across replicates (Fig. 2).This leads to standard deviations in 18 a SO 3 ðOHÞ À =H 2 Oðliq:Þ and DD 0 17 O SO 3 ðOHÞ À =H 2 Oðliq:Þ of only ±0.6‰ and ±0.002‰ at 25 °C (n ¼ 3; Fig. 3), which is near typical analytical precision.Seventh-order polynomial fits for calculating 18 b h;Sc and 17 j h;Sc values [replicate average results for SO 3 (OH) À ] are reported in Tables 1,2.
Resulting 18 a SO 2À 4 =H 2 Oðliq:Þ and 18 a SO 3 ðOHÞ À =H 2 Oðliq:Þ values are in close agreement with all previous experimental and theoretical data for dissolved sulfate and sulfate-bearing minerals (Table S5; Lloyd, 1968;Kusakabe and Robinson, 1977;Chiba et al., 1981;Zeebe, 2010;Schauble and Young, 2021).Comparing dissolved sulfate experimental results to SO 2À 4 predictions yields an RMSE of 4.5‰, much larger than expected analytical precision.However, this RMSE decreases substantially to 1.6‰ when comparing to SO 3 (-OH) À , rather than SO 2À  4 , predictions (Fig. 3A).As pointed out by Mizutani and Rafter (1969) and Zeebe (2010), dissolved sulfate was actually present as SO 3 (OH) À under the experimental conditions used in Lloyd (1968), consistent with our predictions.Anhydrite results of Chiba et al. (1981) are similarly better described by 18 a SO 3 ðOHÞ À =H 2 Oðliq:Þ predictions than by 18 a SO 2À 4 =H 2 Oðliq:Þ predictions.As discussed in Chiba et al. (1981), this likely results from the smaller ionic radius of the Ca 2+ cation in anhydrite relative to the Ba 2+ cation used to precipitate barite for dissolved sulfate isotope analysis; smaller ionic radius and lower cation mass has been shown to lead to larger oxygen isotope fractionation (O'Neil et al., 1969).Chiba et al. (1981) therefore hypothesized that the agreement between dissolved-phase SO 3 (OH) À and solid-phase CaSO 4 fractionation results from similar bonding environments between these two species.This is confirmed by the barite experimental results of Kusakabe and Robinson (1977), which show smaller 18 a BaSO 4 =H 2 Oðliq:Þ relative to the 18 a CaSO 4 =H 2 Oðliq:Þ of Chiba et al. (1981) at a given temperature.Similarly, lattice-dynamic theoretical results of Schauble and Young (2021) predict that BaSO 4 equilibrium fractionation is % 3 to 5‰ smaller than that of CaSO 4 across all temperatures, similar to the observed experimental offset and that seen here between 18 a SO 2À 4 =H 2 Oðliq:Þ and 18 a SO 3 ðOHÞ À =H 2 Oðliq:Þ (Fig. 3A), further supporting the interpretation of Chiba et al. (1981).
Importantly, 18 a SO 2À 4 =H 2 Oðliq:Þ values calculated here between 0 °C and 150 °C are in close agreement with predictions from Zeebe (2010) (their Eq. 5), exhibiting an RMSE of 0.5‰ (Fig. 3A).This agreement occurs despite the conclusion by Zeebe (2010) that MP2 or B3LYP functionals yield unstable geometries and inaccurate scaling factors for hydrated SO 2À 4 .However, the basis sets tested by Zeebe (2010) for B3LYP and MP2 functionals did not include hydrogen atom polarization functions, which are necessary to minimize basis set superposition error when using DFT methods for systems with hydrogen bonds (Novoa and Sosa, 1995).By including hydrogen atom polarization functions, our results never yielded unstable geometries (i.e., imaginary frequencies) for any sulfate species at any water droplet cluster size (Supplementary Data).Furthermore, the agreement between 18 a values predicted here with past experimental and theoretical results confirms our choice of PES scaling factors.In contrast, the use of unscaled 18  4 =H 2 Oðliq:Þ and DD 0 17 O SO 3 ðOHÞ À =H 2 Oðliq:Þ values are nearly identical across all temperatures, differing by a maximum of 0.009‰ at 0 °C (Fig. 3B).Interestingly, both species yield large negative DD 0 17 O predictions (i.e., D 0 17 O values less than that of surrounding water), reaching values as low as À0.199‰ at 0 °C.Predicted frac-tionations for both 18 O and 17 O decrease with increasing temperature, as expected.These results are again nearly identical to the solid-phase BaSO 4 and, especially, CaSO 4 predictions of Schauble and Young (2021), which were calculated using a similar theoretical method but with independent input data (i.e., lattice dynamics; Fig. 3B).

Aqueous sulfite species
Similar to sulfate, sulfite species geometries calculated here are in close agreement with those reported in Eldridge et al. (2016),  4 ; SO3ðOHÞ À , n ¼ 3 replicates] and ''semi empirical" liquid water using 18 b h;Sc and 17 j h;Sc values calculated by Tables 1 and 2. Also shown in panel (A): experimental results for dissolved HSO À 4 over a range of pH conditions (Lloyd, 1968), experimental results for barite mineral (BaSO 4 ) in saline hydrothermal conditions (Kusakabe and Robinson, 1977), experimental results for anhydrite mineral (CaSO 4 ) at high-pressure hydrothermal conditions (Chiba et al., 1981), and ab initio predictions for SO 2À 4 calculated using the Hartree Fock method (Zeebe, 2010).Also shown in both panels: theoretical predictions for BaSO 4 and CaSO 4 calculated using the lattice dynamics method (Schauble and Young, 2021).Panel (B) values correspond to 17 hRL ¼ 0:5305 in the definition of D 0 17 O.
with bond lengths differing by a maximum of 0.03Å and angles differing by a maximum of 2.5°(Fig.1; Table S3).However, unlike for sulfate species, there exist large differences in 18 b h;Sc and 17 j h;Sc between sulfite species (Fig. 2).Across all temperatures, ðHSÞO À 3 consistently exhibits the highest 18 b h;Sc values whereas SO 2À 3 consistently exhibits the lowest values.In contrast, SO 2À 3 displays the highest 17 j h;Sc values whereas SO 2ðaq:Þ displays the lowest values.Seventh-order polynomial fits for calculating 18 b h;Sc and 17 j h;Sc values for all sulfite species are reported in Tables 1,2.
These differences in 18 b h;Sc and 17 j h;Sc lead to large differences in 18 a and DD 0 17 O predictions between different sulfite species and water (Fig. 4).For example, our results predict that ðHSÞO À 3 is 18.2‰ more enriched than SO 2À 3 when both are in equilibrium with water at 25 °C.Similarly, estimated D 0 17 O for SO 2À 3 is 0.090‰ higher than that for SO 2ðaq:Þ when both are in equilibrium with H 2 O at 25 °C; still, all species display D 0 17 O values lower than that of equilibrated water across all temperatures.
Above these temperatures, we predict that SO 2À 3 and SO 2 ðOHÞ À are more depleted in 18 O relative to H 2 O. Similar crossovers have been observed previously for other oxygen-bearing species (e.g., Hayles et al., 2018).In contrast, all other species exhibit 18 a values as high as 26.4‰ at 0 °C and do not display crossover points.Unlike for 18 a; DD 0 17 O predictions never exhibit crossover points and instead trend toward zero at high temperature for all species, as expected (Young et al., 2002).
All existing experimental sulfite 18 a results are in relatively close agreement with predictions calculated here (Fig. 4A, Table S5; Müller et al., 2013a;Wankel et al., 2014, this study).
Experiments were performed at a range of pH values from % 2 to % 10, leading to large differences in isomer relative abundances between experimental conditions.We therefore compare experimental results to predicted ''bulk solution" fractionation; i.e., the temperature-and pH-specific abundance-weighted average frac-Fig.4. Triple-oxygen fractionation factors between liquid water and sulfite species.Predicted (A) 1000 Â lnð 18 aÞ and (B) DD 0 17 O between each sulfite species [SO 2À 3 , ðHSÞO À 3 , SO2ðOHÞ À , and SO2ðaq:Þ] and ''semi empirical" liquid water using 18 b h;Sc and 17 j h;Sc values calculated by Tables 1 and 2. Predicted ''bulk solution" (C) 1000 Â lnð 18 aÞ and (D) DD 0 17 O as a function of pH for the four temperatures in which experimental data exist: 2°C (black), 25°C (blue), 50 °C (red), and 95 °C (yellow).Temperature-and pHspecific relative abundances of each species were calculated using the speciation constants and isomerization ratio of Millero et al. (1989) and Eldridge et al. (2018) (Fig. S14A).Also shown in panel (A) and (C): experimental results for dissolved sulfite over a range of temperature and pH conditions from three laboratories (Müller et al., 2013a;Wankel et al., 2014, this study).Experiments performed at 4 °C and 22 °C are plotted as 2 °C and 25 °C, respectively.Panel (B) and (D) values correspond to 17 hRL ¼ 0:5305 in the definition of D 0 17 O.tionation factors.Speciation diagrams (Fig. S14A) were determined using the dissociation constants of Millero et al. (1989) (based on data from Goldberg and Parker, 1985) and the bisulfite isomerization ratio of Eldridge et al. (2018) assuming a dilute solution (i.e., pure water; see also Horner and Connick, 1986;Littlejohn et al., 1992;Damian et al., 2007, for further discussion on isomerization ratio estimates).This yields an overall experimental vs. predicted bulk solution RMSE of 3.7‰.
Experimental vs. predicted agreement is closest at circumneutral pH when SO 2À 3 and total bisulfite occur in roughly equal proportion, as expected given that experimental results consistently lie between SO 2À 3 and SO 2 ðOHÞ À predictions (Fig. 4A).In contrast, bulk solution theoretical results over-predict low-pH experiment18 a values by as much as 6‰ and under-predict high-pH experiment 18 a values by as much as 10‰ (Fig. 4C; Fig. S14B) due to the fact that experimental results are relatively insensitive to pH.For example, this can be seen by comparing all experimental results to predicted 18 a SO 2 ðOHÞ À =H 2 Oðliq:Þ values; this yields an RMSE of 2.2‰, much lower than for the bulk solution, despite the fact that SO 2 ðOHÞ À alone never dominates the bulk sulfite solution.As discussed at length in Müller et al. (2013a), retaining accurate isotope signatures during experimental conversion to SO 2À 3 and subsequent precipitation as BaSO 3 is challenging due to (i) potential kinetic fractionation effects (ii), incorporation of H 2 O molecules into the hygroscopic salt, and (iii) potential biases toward a particular species or isomer [e.g., SO 2 ðOHÞ À ] during the precipitation reaction.It therefore remains unclear if theory vs. experiment offsets are largely due to (i) inadequacies in the theoretical predictions, (ii) the of unknown (non-equilibrium) experimental isotope effects, or (iii) some combination of these factors.This discrepancy may be reconciled in the future by measuring oxygen-isotope compositions of intact aqueous sulfoxyanions using high-resolution mass spectrometry techniques (e.g., Neubauer et al., 2020).

Aqueous sulfoxylate species
Sulfoxylate species geometries again agree closely with those calculated in Eldridge et al. (2016); bond lengths differ by a maximum of 0.02Å and geometries differ by a maximum of 2.7°(Fig.1; Table S3).Like for sulfite species, there exist large differences in Interestingly, 18 a results for all sulfoxylate species either yield a crossover point [for ðHSÞO À 2 , SO(OH) À , (HS)O(OH), and S(OH) 2 ] or predict 18 O depletion relative to H 2 O across the entire temperature range considered [for SO 2À  2 ].Unlike for sulfoxyanion species at all other redox states, this additionally leads to crossovers in DD 0 17 O; specifically, SO 2À 2 reaches D 0 17 O values that are 0.003‰ higher than that of equilibrated water at 350 °C.Because little is known about the role of sulfoxylate species in the global sulfur cycle, there exist no experimental isotope fractionation results with which we can compare our predictions.Nonetheless, seventh-order polynomial fits for calculating 18 b h;Sc and 17 j h;Sc values for all sulfoxylate species are reported in Tables 1,2.

Aqueous thiosulfate species
Finally, S 2 O 2À 3 also displays a similar geometry to that calculated in Eldridge et al. ( 2016) (all other thiosulfate species were not included in their study); bond lengths differ by a maximum of 0.01Å and geometries differ by a maximum of 0.4°(Fig.1; Table S3).Like SO 3 (OH) À , we additionally calculated triplicate S 2 O 2 ðOHÞ À geometries and fractionation factors.Replicate geometries show similar bond lengths, with differences reaching 0.02Å for S-O bonds, 0.03Å for S-S and S-(OH) bonds, and 0.04Å for O-H bonds.Bond angle differences between replicates reach 2.2°for O-S-O, 3.2°for S-O-H, 3.6°for S-S-(OH), 3.8°for S-S-O, and 5.2°f or O-S-(OH).Despite these geometric differences, 18 b h;Sc and 17 j h;Sc are again nearly identical across replicates (Fig. 2), leading to standard deviations in 18 a S 2 O 2 ðOHÞ À =H 2 Oðliq:Þ and DD 0 17 O S 2 O 2 ðOHÞ À =H 2 Oðliq:Þ of only ±0.2‰ and ±0.001‰ at 25 °C (n ¼ 3; Fig. 6), well within analytical precision.
Few reliable experimental equilibrium exchange fractionation factor estimates exist for thiosulfate species.We compare our predictions to experimental results also presented in this study; however, it is likely that oxygen isotope equilibrium was not reached under some experimental conditions.We therefore exclude from our comparison any experimental results that are statistically identical to the Na 2 S 2 O 3 starting material d 18 O value (p < 0:01; two-tailed t test; Table S5).This leads to an RMSE between our predictions and retained experimental data of 2.2‰ when comparing to 18 a S 2 O 2À 3 =H 2 Oðliq:Þ , 3.2‰ when comparing to 18 a ðHSÞSO À 3 =H 2 Oðliq:Þ , and 3.6‰ when comparing to 18 a S 2 O 2 ðOHÞ À =H 2 Oðliq:Þ .Interestingly, unlike for sulfate and sulfite species, thiosulfate RMSE is highest when comparing to the S 2 O 2 ðOHÞ À isomer, although the exact proportion of protonated thiosulfate at our experimental conditions remains unknown.However, this RMSE should be considered a maximum estimate since it remains possible that isotope exchange remained incomplete in these experiments.Like all other sulfoxyanion species, no experimental DD 0 17 O data exist with which we can compare our theoretical results.

Discussion and implications
We discuss how these equilibrium fractionation factors update our understanding of several sulfur-cycle processes-including pyrite oxidation, MSR, thiosulfate disproportionation, and hydrothermal anhydrite precipitation-that represent the major sulfur fluxes on Earth's surface.For each process, we assess whether equilibrium predictions support or refute a certain mechanistic pathway.We focus specifically on the possible incorporation of atmospheric O 2 into sulfate, as this has been previously invoked to explain fluvial, lacustrine, and marine sulfate D 0 17 O values (e.g., Bao et al., 2008;Crockford et al., 2018;Killingsworth et al., 2018;Crockford et al., 2019).

Pyrite oxidation
Sulfate dissolved in modern rivers and preserved in ancient mineral deposits often displays negative D 0 17 O values (see Crockford et al., 2019, for compilation).This result is canonically interpreted to reflect incorporation of oxygen sourced from a mixture of water and dissolved O 2 -which carries a negative massindependent 17 O signal (Thiemens and Lin, 2021)-into sulfate during pyrite oxidation (Fig. 7A; e.g., Bao et al., 2008;Crockford et al., 2018;Killingsworth et al., 2018).However, dissolved sulfate in modern rivers draining pyrite-rich lithologies has recently been shown to exhibit D 0 17 O values equal to or slightly higher than those of concomitant water, questioning this mechanistic interpretation (Hemingway et al., 2020).
While several aspects of the pyrite oxidation mechanism remain unknown or underconstrained (e.g., Schoonen et al., 2010), it is generally accepted that pyrite sulfur is oxidized via a multi-step electron transfer process (the so-called ''semiconductor" model; Williamson and Rimstidt, 1994).Accordingly, pyrite sulfur acts as an anode that iteratively donates electrons to cathodic iron atoms.Electropositive sulfur is subsequently subject to nucleophilic attack by H 2 O or OH À , forming sulfoxy species and releasing H + to solution.Direct O 2 incorporation into sulfate is thus inconsistent with the semi-conductor model, although the importance of alternative, isotopically unique nucleophiles such as H 2 O 2 remains unknown (Schoonen et al., 2010;Hemingway et al., 2020).3 ; ðHSÞSO À 3 ; and S2O2ðOHÞ À , n ¼ 3 replicates] and ''semi empirical" liquid water using 18 b h;Sc and 17 j h;Sc values calculated by Tables 1 and 2 Here, we instead hypothesize that dissolved sulfate d 18 O and D 0 17 O values can reflect intermediate sulfoxyanion oxygenisotope equilibrium with water and subsequent (possibly microbially mediated) dissolved-phase oxidation, either during initial pyrite oxidation or downstream redox cycling.We test this hypothesis using recently reported triple-oxygen isotope values for a time-series of Mississippi River sulfate collected at Baton Rouge, Louisiana, USA (Killingsworth et al., 2018).Pyrite oxidation-derived sulfur is released to solution either as sulfite, thiosulfate, or sulfate depending on pH (e.g., Rimstidt and Vaughan, 2003).Furthermore, pyrite surfaces have been shown to catalyze thiosulfate oxidation to sulfite (Xu and Schoonen, 1995), which exhibits rapid oxygen-isotope exchange under circumneutral to acidic pH values (Betts and Voss, 1970;Wankel et al., 2014).If we assume the final oxygen atom is derived from water with a negligible kinetic isotope effect (cf., Müller et al., 2013b;Cao and Bao, 2021), then Mississippi River sulfate isotope compositions can be explained by sulfite-water equilibrium isotope exchange at pH 7 followed by terminal oxidation to sulfate (Fig. 7B).In contrast, these data are not consistent with a terminal oxygen atom derived from dissolved O 2 nor with thiosulfate isotope equilibrium followed by disproportionation (Fig. 7B), implicating sulfite as the critical intermediate sulfoxyanion.
Interestingly, this result is strongly dependent on solution pH.For example, low-pH equilibrium expected to occur in acid-mine drainage settings would exhibit much lower D 0 17 O near À0.175‰ at 25 °C (Fig. 4D).In contrast, alkaline solutions such as those expected during carbonate-or silicate-buffered rock weathering would display smaller D 0 17 O offsets from water at the same temperature.Thus, in addition to temperature, weathering-solution pH likely exhibits a strong control on resulting dissolved sulfate oxygen-isotope compositions.
Still, large uncertainties persist.For example, both Rimstidt and Vaughan (2003) and Kohl and Bao (2011) showed that the majority of released S can occur as thiosulfate in circumneutral to alkaline aerobic and anaerobic pyrite oxidation experiments, even after several weeks.Similarly, Hemingway et al. (2020) showed that pyrite oxidation-derived sulfate can retain anomalously positive D 0 17 O values-possibly sourced from atmospheric H 2 O 2 -although this signal is overprinted by downstream processes (e.g., biogeochemical sulfate recycling).While the mechanism of H 2 O 2 signal Fig. 7. New and canonical interpretations of sulfate isotope compositions in the oxidative sulfur cycle.Predicted sulfate triple-oxygen isotope compositions following: (A) canonical interpretations of experimentally observed net pyrite oxidation isotope effects, (B) sulfite-water isotope equilibrium at 10 °C and pH 7 followed by oxidation to sulfate, (C) thiosulfate-water isotope equilibrium at 10 °C followed by disproportionation to sulfate and hydrogen sulfide (Eq.19).Mississippi River dissolved sulfate isotope compositions are used as an example dataset to test these predictions.Markers common to all panels are as follows: blue squares = average Mississippi River water at Arkansas City, AR, USA (1984-1987;n ¼ 10;Coplen and Kendall, 2000;Killingsworth et al., 2018); red diamonds = atmospheric O 2 (Sharp and Wostbrock, 2021); red squares = dissolved oxygen in equilibrium with atmospheric O 2 at 10 °C (Benson and Krause, 1984;Luz and Barkan, 2009); gray circles = individual Mississippi River dissolved sulfate samples (2009-2014, n ¼ 38;Killingsworth et al., 2018); white squares = Mississippi River dissolved sulfate average composition; gray lines = fractionation trajectories; black lines = mixing trajectories.Originally reported Mississippi River sulfate D 0 17 O values have been shifted up by 0.07‰ as recommended by Cao and Bao (2021) to place results closer to the SMOW-SLAP calibration scale.For panel (A), 18 a values of net pyrite oxidation with all O from either H 2 O or O 2 are taken from Balci et al. (2007); corresponding DD 0 17 O values have not been measured and are assumed here to follow kinetic fractionation lines with slopes defined by the reduced masses of reactants (cf., Cao and Bao, 2021, their (Balci et al., 2007;Kohl and Bao, 2011).For panel (B), we assume T-and pH-specific bulk solution fractionation factors (see main text) and for panel (C) we assume fractionation with S2O 2À 3 .We additionally assume that sulfite/thiosulfate is quantitatively oxidized/disproportionated such that any kinetic fractionation is not expressed; fractionation of H 2 O during oxidation/disproportionation is not constrained but is thought to be of minor importance (cf., Müller et al., 2013b).All D 0 17 O values correspond to 17 hRL ¼ 0:5305 and are reported on the SMOW-SLAP calibration scale whenever possible (Sharp and Wostbrock, 2021).Importantly, only sulfite equilibrium followed by oxidation with terminal oxygen from H 2 O can explain observed Mississippi River results.incorporation into sulfate remains unknown, Fe(II)-mediated H 2 O 2 degradation has recently been shown to follow a mass law of 0:5347 AE 0:0006 (Sutherland et al., 2022); this is larger than 17 h RL and thus leads to 17 O enrichment in degradation products relative to reactant H 2 O 2 .It is thus possible that the incorporation of Fe(II)mediated H 2 O 2 degradation products into sulfate could also explain observed positive D 0 17 O values.In contrast, Cao and Bao (2021) reinterpreted positive D 0 17 O values in riverine sulfate as reflecting kinetic isotope fractionation rather than H 2 O 2 incorporation, although this interpretation relies on knowledge of kinetic triple-oxygen isotope fractionation factors, which remain unconstrained.Finally, the mechanism proposed here cannot explain D 0 17 O values as low as À1.0‰ observed in Neoproterozoic sulfate deposits (Bao et al., 2008).Future laboratory-and field-based work is therefore crucial to constrain in situ environmental parameters such as pH, thiosulfate and sulfite concentrations and isotope compositions, and H 2 O 2 concentrations at the site of pyrite oxidation.

Microbial sulfate reduction
Whereas pyrite oxidation is the dominant source of sulfate to Earth's surface, MSR and subsequent pyrite formation in marine sediments represents the dominant sink.For the purpose of tracing oxygen isotopes, MSR can be interpreted as following the simplified intracellular reaction network: where APS refers to the adenosine phosphosulfate enzyme complex; we exclude terminal reduction of SO 2À 3 to H 2 S since this does not involve oxygen exchange (Zeebe, 2010;Wankel et al., 2014;Wing and Halevy, 2014;Bertran et al., 2020).Furthermore, the oxidative back-reaction from APS to SO 2À  4 can occur either enzymatically or abiotically (Bertran et al., 2020;Benkovic and Hevey, 1970).In the enzymatic case, sulfate is released via nulcleophilic attack on the APS phosphorus atom; one of four oxygen atoms in resulting sulfate is thus derived from the phosphate group of adenosine monophosphate (AMP; Brunner et al., 2012).In the abiotic case however, sulfate is nearly quantitatively released by a unimolecular elimination reaction, leading to sulfate with three oxygen atoms derived from sulfite and one oxygen atom directly sourced from water (Benkovic and Hevey, 1970).At Earth-surface conditions, neither sulfate nor APS are known to exchange oxygen atoms with ambient water (Chiba and Sakai, 1985;Kohl et al., 2012), implicating sulfite and AMP phosphate as the primary species by which oxygen-isotope exchange can occur.
Depending on thermodynamic conditions and substrate concentrations-and thus sulfate reduction rates-MSR will exist between purely kinetic (i.e., unidirectional forward fluxes in Eq. 18) and and equilibrium (i.e., equal forward and backward fluxes in Eq. 18) limits (Wing and Halevy, 2014).At the kinetic limit, all generated SO 2À 3 is completely reduced to H 2 S within the cell; at the equilibrium limit however, isotopically equilibrated SO 2À 3 can back-react to SO 2À  4 .The exact position between these limits therefore determines the degree to which equilibrium isotope exchange during MSR (including AMP phosphate-derived oxygen for the enzymatic oxidation reaction) impacts marine sulfate oxygen isotope compositions (Bertran et al., 2020).Using a thermodynamic model, Wing and Halevy (2014) estimated that MSR in marine sediments likely approaches the equilibrium limit, even in coastal regions with high sulfate reduction rates.This has since been confirmed in several field localities by tracking porewater sulfate 33 S and 34 S evolution with sediment depth; results are inconsistent with MSR operating in the kinetic regime (Masterson et al., 2018;Masterson et al., 2022).Several aspects of equilibrium oxygen iso-tope exchange during MSR remain unknown or underconstrained, particularly regarding the timescale of AMP phosphate oxygen exchange (Brunner et al., 2012;Chang and Blake, 2015).Nevertheless, here we use our theoretical predictions to estimate two endmember scenarios that may prove useful for interpreting field results: (i) Full expression of the sulfate-water equilibrium fractionation factor.At the equilibrium limit, Bertran et al. (2020) estimated that the abiotic elimination pathway likely dominates the oxidative APS back-reaction.If true, this implies that SO 2À 4 (at circumneutral pH) regeneration involves an activated transition state resembling AMP Á Á Á SO 3 Á Á Á H 2 O, which may rapidly exchange oxygen atoms with water (Benkovic and Hevey, 1970).Such a direct sulfate-water exchange mechanism has been invoked previously to explain the similarity between observed and theoretical sulfate d 18 O predictions (Zeebe, 2010); however, this similarity may be fortuitous rather than mechanistic (Brunner et al., 2012).Nevertheless, our results predict direct sulfate-water equilibrium during MSR would push sulfate to higher d 18 O and lower D 0 17 O with a 17 h SO 2À 4 =H 2 O value ranging from 0.5240 to 0.5243 between 0 and 25 °C.For the modern ocean, this equates to an equilibrium MSR sulfate composition with d 18 O ranging from 29.1 to 24.8‰ and D 0 17 O ranging from À0.187 to À0.154‰ between 0 and 25 °C.
(ii) A weighted-average of the sulfite-water (75%) and AMP phosphate-water (25%) equilibrium fractionation factors.If we assume that the residence times of sulfite and AMP phosphate are long enough such that both reach isotopic equilibrium with water, the overall expressed sulfate-water fractionation should represent a weighted average of these two species (Brunner et al., 2012).Here, we ignore any additional fractionation between

SO 2À
3 and APS and between APS and SO 2À 4 since their triple-oxygen fractionation factors remain unknown, although these steps may prove important.Unfortunately, AMP phosphate triple-oxygen isotope fractionation factors have also not been measured or theoretically predicted.We instead use theoretical fluoroapatite predictions from Schauble and Young (2021); these are in close agreement with experimentally derived dissolved phosphate 18 O fractionations and should thus serve as a useful first approximation (Chang and Blake, 2015).Assuming bulk-solution sulfite fractionation at pH 7 (Fig. S14A; Millero et al., 1989;Eldridge et al., 2018), this results in an expressed 17 h ranging from 0.5176 to 0.5194 for temperatures between 0 and 25 °C.For the modern ocean, this leads to equilibrium MSR-derived sulfate d 18 O ranging from 12.8‰ to 9.8‰ and D 0 17 O ranging from À0.122‰ to À0.095‰.However, these d 18 O values are up to 15‰ lower than porewater observations, suggesting additional fractionation during reoxidation back-reactions (Zeebe, 2010;Brunner et al., 2012).Waldeck et al. (2019) recently reported that marine sulfate is described by d 18 O and D 0 17 O values of 8.67±0.21‰and À0.016 ± 0. 017‰, respectively.For both scenarios described here, this is consistent with a mixture between MSR recycling and an 18 Odepleted, 17 O-enriched input, presumably dissolved riverine sulfate (Hemingway et al., 2020).To estimate the relative importance of MSR and riverine inputs in setting observed marine signals, Waldeck et al. (2019) hypothesized a d 18 O value for the MSR end member ranging from 24‰ to 27‰, within our ''full sulfate expression" scenario range.However, that study also assumed an MSR 17 h value between 0:527 and 0:5305, substantially higher than the range estimated here for both scenarios (i.e., total range of 0:5176 to 0:5243).The 17 h values of Waldeck et al. (2019) would lead to a less negative MSR D 0 17 O prediction and thus-all else being equal-a greater contribution of this end member to observed marine signals relative to that predicted using the fractionation factors proposed here.Still, several aspects of the con-trols on marine sulfate D 0 17 O values remain unknown, particularly globally integrated riverine input values.Future studies are therefore needed to further constrain these signals.
Although several aspects of MSR triple-oxygen isotope fractionation remain under-constrained, one common feature of all predictions derived here is the generation of sulfate with D 0 17 O values lower than those of ambient water.It thus becomes clear that observed slightly negative D 0 17 O values of fluvial and marine sulfate (Killingsworth et al., 2018;Waldeck et al., 2019) can be explained purely by sulfate-water equilibrium during MSR operating near the equilibrium limit without the need to invoke incorporation of anomalously 17 O-depleted dissolved O 2 .

Thiosulfate disproportionation
Thiosulfate is produced by the partial oxidation of sulfide at redox boundaries, for example when aerated water penetrates into anoxic hot springs (Xu et al., 1998) or marine sediments (Jørgensen, 1990).At circumneutral pH, bulk thiosulfate isotopically exchanges oxygen with surrounding water with a half life on the order of hours (the presence of SO 2 ðOHÞ À can also act as a minor catalyst for this exchange; Pryor and Tonellato, 1967).Thiosulfate can then disproportionate either biologically or abiotically, but these two mechanisms follow unique reaction pathways.Microbial thiosulfate disproportionation has been shown to proceed as (Jørgensen, 1990;Finster et al., 1998): In contrast, abiotic disproportionation is thought to involve a suite of reactions with the rate-limiting step described by the bimolecular reaction (Johnston and McAmish, 1973;Xu and Schoonen, 1995;Xu et al., 1998): Assuming (i) reactant S 2 O 2À 3 reaches isotopic equilibrium with water, (ii) microbial disproportionation proceeds unidirectionally (i.e., no back-reaction in Eq. 19), and (iii) the final oxygen atom is derived from water with a negligible kinetic isotope effect (cf., Cao and Bao, 2021), then microbial disproportionation will yield product sulfate with d 18 O and D 0 17 O values %20‰ higher and %0.15‰ lower than surrounding water, respectively (Fig. 7C).While each of these assumptions must be thoroughly validated (e.g., the residence time of thiosulfate in a given environment may be too short to reach isotopic equilibrium), this calculation nevertheless provides a useful end-member scenario to interpret environmental data.
In contrast, abiotic thiosulfate disproportionation produces sulfite, which will itself rapidly exchange oxygen isotopes with surrounding water (Fig. 7B).The original thiosulfate isotope signature is therefore overprinted by sulfite oxygen exchange, independent of the degree to which reactant thiosulfate and water reached isotope equilibrium.Subsequent oxidation will yield sulfate with an isotope signature that reflects a mixing between equilibrated sulfite and the final oxygen atom source (in addition to any kinetic effects; Cao and Bao, 2021).Still, the importance of pH and temperature on thiosulfate oxygen-isotope equilibrium warrants further experimental study (cf., Eldridge et al., 2021, for related sulfur-isotope discussion).

Hydrothermal oxygen isotope exchange
Finally, we briefly consider hydrothermal oxygen-isotope exchange between water and sulfate.Although exchange is negligible at Earth-surface conditions-even over billion-year timescales-exchange rates increase drastically at elevated temperatures characteristic of hydrothermal vents (Chiba and Sakai, 1985).Measured d 18 O values of laboratory hydrothermally precipitated anhydrite (CaSO 4 ) are consistent with our SO 3 (OH) À fractionation predictions (Fig. 3A; see Sec. 4.3 and Chiba et al., 1981;O'Neil et al., 1969, for discussion on the comparison between SO 3 (OH) À and CaSO 4 fractionation factors), supporting equilibrium exchange under these conditions, as has been suggested previously (Lloyd, 1968;Chiba et al., 1981;Zeebe, 2010).However, no hydrothermal anhydrite D 0 17 O records currently exist to our knowledge.Nevertheless, anhydrite d 18 O has been used as a proxy for alteration fluid temperature at the time of mineral precipitation, although this requires that alteration fluid oxygen-isotope composition is accurately known (e.g., Teagle et al., 1998) and that anhydrite-water equilibrium is reached in natural vent settings (cf., Chiba et al., 1998).Still, if the temperature of isotope exchange can be constrained independently (e.g., by traditional or ''clumped" isotope thermometry of co-existing carbonates; Weinzierl et al., 2018) and equilibrium exchange can be ensured (e.g., by additionally measuring d 34 S; Chiba et al., 1998), then hydrothermal anhydrite veins preserved in oceanic crust and obducted ophiolites may record the triple-oxygen isotope composition of alteration fluid in the geologic past.This hypothesis remains speculative but warrants further study.

Conclusions
The triple-oxygen isotope composition of sulfate (d 18 O and D 0 17 O)-both in modern aquatic systems and in geologically preserved sulfate-bearing minerals-is becoming a common tool to constrain sulfur-cycle processes.However, equilibrium oxygenisotope fractionation factors between water and intermediate sulfoxyanion species remain largely unknown.Here, we estimate fractionation factors for a suite of sulfoxyanions-including several protonated species-using a quantum-chemistry computational approach; our results are in relatively good agreement with all available experimental constraints, and we discuss possible experimental and theoretical means of improving agreement, especially for the sulfite system.We highlight the potential importance of short-lived thiosulfate and, especially, sulfite species in setting d 18 O and D 0 17 O values of sulfate produced by several abiotic and biological processes (e.g., pyrite oxidation, MSR, thiosulfate disproportionation, anhydrite precipitation).Importantly, when equilibrium sulfite or thiosulfate fractionation factors are expressed, resulting sulfate can exhibit D 0 17 O values up to %0.20‰ lower than equilibrated water.Slightly negative D 0 17 O values thus do not require incorporation of alternative oxygen sources such as dissolved O 2 , as has been previously assumed.This result carries implications for the interpretation of isotope signals recorded in geologic sulfate archives.

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.

Table 1
Seventh-order polynomial fit parameters to the equation: ln 18 b h;Sc À Á ¼ P p½i T i , where T is in Kelvin.n.a.= not applicable; RMSE = root mean square error between ''true" values and polynomial fits; redox = sulfur atom redox state.Polynomial fits were determined over the temperature range 0 °C to 375 °C.All results are calculated using B3LYP/6-31+G(d,p) harmonic frequencies scaled to redox state-specific CCSD/aug-cc-pVTZ harmonic frequncies (Section 3).

Table 2
Seventh-order polynomial fit parameters to the equation: 17 j h;Sc ¼ P p½i T i , where T is in Kelvin.n.a.= not applicable; RMSE = root mean square error between ''true" values and polynomial fits; redox = sulfur atom redox state.Polynomial fits were determined over the temperature range 0 °C to 375 °C.All results are calculated using B3LYP/6-31+G(d,p) harmonic frequencies scaled to redox state-specific CCSD/aug-cc-pVTZ harmonic frequncies (Section 3)..D. Hemingway, M.L. Goldberg, K.M.Sutherland et al.

Table 1 )
. We assume final sulfate contains 75% O from H 2 O and 25% O from O 2 , consistent with previous interpretations