The sulfur isotope evolution of magmatic-hydrothermal fluids: insights into ore-forming processes

Metal-rich fluids that circulate in magmatic-hydrothermal environments form a wide array of economically significant ore deposits. Unravelling the origins and evolution of these fluids is crucial for understanding how Earth’s metal resources form and one of the most widely used tools for tracking these processes is sulfur isotopes. It is well established that S isotopes record valuable information about the source of the fluid, as well as its physical and chemical evolution (i.e. changing pH, redox and temperature), but it is often challenging to unravel which of these competing processes drives isotopic variability. Here we use thermodynamic models to predict S isotope fractionation for geologically realistic hydrothermal fluids and attempt to disentangle the effects of fluid sources, physico-chemical evolution and S mineral disequilibrium. By modelling a range of fluid compositions, we show that S isotope fingerprints are controlled by the ratio of oxidised to reduced S species (SO4 2 /H2S), and this is most strongly affected by changing temperature, fO2 and pH. We show that SO4 2 /H2S can change dramatically during cooling and our key insight is that S isotopes of individual sulfide or sulfate minerals can show large fractionations (up to 20‰) even when pH is constant and fO2 fixed to a specific mineral redox buffer. Importantly, while it is commonly assumed that SO4 2 /H2S is constant throughout fluid evolution, our analysis shows that this is unlikely to hold for most natural systems. We then compare our model predictions to S isotope data from porphyry and epithermal deposits, seafloor hydrothermal vents and alkaline igneous bodies. We find that our models accurately reproduce the S isotope evolution of porphyry and high sulfidation epithermal fluids, and that most require magmatic S sources between 0 and 5‰. The S isotopes of low sulfidation epithermal fluids and seafloor hydrothermal vents do not fit our model predictions and reflect disequilibrium between the reduced and oxidised S species and, for the latter, significant S input from seawater and biogenic sources. Alkaline igneous fluids match model predictions and confirm magmatic S sources and a wide range of temperature and redox conditions. Of all these different ore deposits, porphyry and alkaline igneous systems are particularly well-suited to S isotope investigation because they show relationships between redox, alteration and ore mineralogy that could be useful for exploration and prospecting. Ultimately, our examples demonstrate that S isotope forward models are powerful tools for identifying S sources, flagging disequilibrium processes, and validating hypotheses of magmatic fluid evolution. 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/ licenses/by/4.0/).


INTRODUCTION
Most magmatic-hydrothermal systems on Earth contain fluids with significant concentrations of sulfur. In these https://doi.org/10.1016/j.gca.2020.07.042 0016-7037/Ó 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). environments S exists in a variety of oxidation states (S 2À to S 6+ ) and phases (melt, minerals, fluids and gases), and participates in many chemical reactions. As a result, S often shows large isotopic fractionations between co-existing Sbearing minerals and is well suited to geochemical investigation (Seal, 2006;Schauble, 2008). At active volcanic centres, S isotopes have been used to provide critical insights into magma chamber processes and the dynamics of hydrothermal systems (Marini et al., 2011), while at extinct magmatic centres, S isotopes have been used to evaluate magma sources Penniston-Dorland et al., 2012) and understand fluid evolution (Rye, 2005;Hutchison et al., 2019).
S isotopes have been particularly beneficial in the study of Earth's metal resources because many ore deposits are derived from S-and metal-rich fluids that circulate above and around magmatic intrusions (Seal, 2006). Since most economically important metals are chalcophile, reduced S plays an essential role complexing and transporting these metals (Mountain and Seward, 2003;Stefánsson and Seward, 2004;Pokrovski et al., 2015) and may also control their precipitation as stable sulfide minerals (e.g. chalcopyrite, CuFeS 2 , in porphyry deposits). Thus, metal mobility in ore-forming environments is intimately linked to S speciation which is encoded in the isotopes of S-bearing minerals (Ohmoto, 1972;Rye, 2005;Seal, 2006). Key ore deposits where these techniques have been applied include: Porphyry deposits -major sources of Cu, Mo, W, Sn and Au, which are derived from hydrothermal fluids exsolved from calc-alkaline magmas and precipitated at depths of 1-6 km (Sillitoe, 2010;Wilkinson, 2013). Epithermal deposits -important sources of Au, Ag, Cu, Pb and Zn, which form in shallow magmatichydrothermal environments (generally <1 km) in a variety of tectonic settings (Sillitoe and Hedenquist, 2003) Volcanogenic massive sulfide (VMS) ore deposits and modern seafloor massive sulfide (SMS) deposits -sources of Fe, Cu, Zn, Pb, Au, and Ag, which are associated with seafloor hydrothermal systems and form when convecting hydrothermal fluids are expelled into the ocean and precipitate their metal load (Herzig and Hannington, 1995;Shanks, 2001;Piercey, 2011). Alkaline igneous deposits -key sources of rare earth and high field strength elements (REE and HFSE, including Nb, Ta, Zr, Hf and Sc), which often form from latestage magmatic-hydrothermal fluids that have undergone protracted fractional crystallization and differentiation (usually in an intracontinental rift setting, Marks and Markl, 2015;Dostal, 2017;. A compilation of d 34 S (a ratio of 34 S and 32 S) is shown for each of these metal deposits in Fig. 1. The histograms show that both reduced (sulfide) and oxidised (sulfate) minerals are present in these magmatic-hydrothermal deposits. Fig. 1. Comparison of sulfur isotope results (d 34 S) from porphyry (a), epithermal (b), modern seafloor massive sulfide (SMS) and alkaline igneous (d) ore systems. Histograms showing sulfide (red) and sulfate (blue) mineral d 34 S as filled and outlined bars, respectively. The two histograms in each subplot are binned separately according to the fraction of the total population (with N indicating the number of data in each population of sulfides and sulfates). Sulfide d 34 S is centred around typical magmatic values of 0 ± 5‰ (Chaussidon et al., 1989;Labidi et al., 2012Labidi et al., , 2013Labidi et al., , 2015. In each case the wide spread in d 34 S reflects variations in the fluid source (i.e. the proportion of magmatic to externally-derived fluids), and the offsets between sulfide and sulfate potentially encode the variable physico-chemical properties of the fluids (i.e. differences in fluid pressure, temperature, fO 2 and pH). Data and references for all samples are provided in the Supplementary Information (Tables S1-4). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) In general, sulfides cluster around 0‰ (overlapping the range of most mantle-derived rocks, i.e. 0 ± 5‰, Chaussidon et al., 1989;Labidi et al., 2012Labidi et al., , 2013Labidi et al., , 2015 while sulfates cluster at higher d 34 S (note that the heavy 34 S isotope forms stronger stiffer bonds with oxidised S species and therefore preferentially substitutes into the sulfate phase, Urey, 1947;Thode et al., 1953Thode et al., , 1954Schauble, 2008).
Although each of these deposits is precipitated by hot, metal-rich hydrothermal fluids advecting from or convecting around a deep magmatic system, our compilation ( Fig. 1) reveals clear variations in both the average and range of d 34 S values between the different deposit types. These contrasting S isotope values can be explained by a number of features. Firstly, the S source of the hydrothermal fluids could be variable; S could be magmatic in origin ($0‰), mixed with externally-derived S (such as seawater, $21‰, Rees et al., 1978;Tostevin et al., 2014) or leached from the local crust. Secondly, the physico-chemical parameters of the ore-forming fluids could vary. In this case differences in fluid pressure (P), temperature (T), pH and oxygen fugacity (fO 2 ) could lead to variable proportions of reduced to oxidised S and hence differences in S isotope fractionation between the magmatic-hydrothermal systems (Ohmoto, 1972). Finally, the sulfate and sulfide may have precipitated from different fluids such that the isotopes of these minerals record variable S sources and/or physicochemical evolution of different generations of fluid (this would lead to bulk rock sulfate-sulfide disequilibrium). In short, although S isotopes provide powerful insights into fluid sources, evolution and disequilibria in magmatichydrothermal deposits, the key challenge is disentangling these competing effects.
One new and innovative approach to this problem is to perform multiple S isotope analyses (i.e. 36 S, 34 S, 33 S and 32 S, Farquhar et al., 2000) and identify mass independent fractionations of S (MIF-S). Archean rocks show distinctive MIF-S and numerous studies have used this to fingerprint different S reservoirs that contribute to Archean ore deposits (particularly VMS, e.g. Golding et al., 2011;Jamieson et al., 2013;Sharman et al., 2015), as well as Proterozoic ore bodies that are derived, in part, from Archean S sources (Selvaraja et al., 2017, Bolhar et al., 2020. Multiple S isotope analyses have also enhanced our understanding of modern seafloor hydrothermal systems because small but detectable MIF-S can be used to confirm the input biogenic S (Ono et al., 2007;Peters et al., 2011;Aoyama et al., 2014;McDermott et al., 2015). These efforts have been bolstered by the development of new microanalytical techniques and reference materials to collect in situ (multiple) S isotope datasets at high spatial-resolution (using, Secondary Ion Mass Spectrometry and laser ablation ICP-MS, e.g. Mason et al., 2006;Mandeville et al., 2009;Brueckner et al., 2015;Lode et al., 2015;LaFlamme et al., 2016LaFlamme et al., , 2018aLaFlamme et al., , 2018bCaruso et al., 2018). Although these new techniques represent major advances, they still have their limitations. For example, MIF-S sources are unlikely to have been incorporated into most post-Archean ore systems, and since most magmatichydrothermal environments have temperatures )100°C this precludes biotic MIF-S processes. Even when high spatial-resolution analyses are applied there is still ambiguity in how to interpret d 34 S and whether the values are an intrinsic feature of the source or related to a magmatichydrothermal process.
Another widely used technique for visualizing d 34 S results and evaluating fluid source and evolution is to consider the d 34 S of co-existing sulfate and sulfide phases. d 34 S sulfate -d 34 S sulfide diagrams were first considered by Fifarek and Rye (2005) and have been summarised in subsequent reviews of S isotope geochemistry (e.g. Rye, 2005;Seal, 2006;Marini et al., 2011). An illustration of this approach is provided in Fig. 2  /H 2 S = 1). The red dashed lines represent an identical case but with a d 34 S P S of 5‰. Lines of equivalent temperature (grey) were calculated from equilibrium fractionation factors of Ohmoto and Lasaga (1982). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Although these diagrams provide a method for unravelling the S source and fluid redox (SO 4 2À /H 2 S), it is necessary to assume that both the d 34 S P S and SO 4 2À /H 2 S remained constant during cooling and that mineral d 34 S faithfully records the d 34 S of the fluid. Studies which have applied these methods to natural systems have come to varying conclusions about the usefulness of this approach. For example, some deposits define strong linear d 34 S sulfate -d 34 S sulfide arrays and appear to provide good evidence for isotopic equilibrium (Fifarek and Rye, 2005;Rye, 2005), while data from other deposits are much more scattered and are taken as evidence for sulfate-sulfide disequilibrium (Shelton and Rye, 1982;Eastoe, 1983). The assumption of a fixed SO 4 2À /H 2 S during fluid decompression and cooling is particularly problematic because few minerals actually record this value (e.g. Hettmann et al., 2012;Konecke et al., 2017) and, moreover, it is widely acknowledged that magmas, and the fluids they release, show significant variations in redox during ascent (Burgisser and Scaillet, 2007).
Ultimately, if we are to use S isotopes to understand the origin, precipitation, and alteration of metal deposits, we require robust thermodynamic constraints on the S isotope evolution of geologically realistic magmatic-hydrothermal fluids. Although there have been significant recent developments in analytical techniques (above), there have been comparatively few attempts to develop S isotope models of magmatic-hydrothermal systems (e.g. Janecky and Shanks, 1988). To date the most significant contribution to this issue is the study of Ohmoto (1972) who developed a series of equations to calculate fluid SO 4 2À /H 2 S and d 34 S for fixed temperature, pH and fO 2 values. Our paper develops this approach. By tying the fluid fO 2 to wellknown mineral redox buffers (e.g. nickel-nickel oxide and quartz-fayalite-magnetite, Frost, 1991) we calculate SO 4 2À / H 2 S and d 34 S over a wider range temperature, pressure, and pH conditions. We first evaluate the d 34 S evolution of a variety of generic hydrothermal fluids. We then consider the d 34 S systematics of typical porphyry, epithermal, SMS and alkaline ore forming fluids using appropriate P-T-fO 2 -pH parameters which compare to natural data. Finally, we explore whether d 34 S forward models can be used to identify magmatic-hydrothermal systems with anomalous non-magmatic S sources, flag disequilibrium processes, and test hypotheses of magmatic fluid evolution.

METHODS
Our modelling builds on the work of Ohmoto (1972), who used thermochemical data and isotopic fractionation factors to calculate S speciation and d 34 S in fluid and mineral phases. This approach treats the hydrothermal fluid as a concentrated alkali-chloride solution, a reasonable assumption given that chloride is the dominant ligand in fluid inclusions from virtually all ore-forming systems (e.g. Smith and Henderson, 2000;Seward et al., et al., 2014). The main aqueous S species considered by Ohmoto (1972) are H 2 S, HS -, S 2-, SO 4 2-, HSO 4 -, KSO 4 and NaSO 4 -, and their relative abundance is expressed via the following reactions: NaSO For a solution containing these various S species (i), the total sulfur isotope of the fluid (d 34 S P S ) can be related to the individual S isotope compositions (d 34 S i ) and mole fractions (X i ) through the equation: Expressing the d 34 S of each S species relative to d 34 S H 2 S , an isotope fractionation factor (D i ) can be used, such that: By combining equations (7) and (8), Ohmoto (1972) showed that: Thus, the isotopic composition of an individual S species in a hydrothermal fluid at equilibrium is controlled by: (1) the total sulfur isotope of the fluid (d 34 S P S ); (2) the mole fractions (X i ) of the various S species and (3) their specific fractionation factors (D i ). To solve Eq. (9), we use isotope fractionation factors of Ohmoto and Rye (1979), for sulfate species, and Ohmoto and Lasaga (1982), for all other S species, and calculate the mole fractions following Eqs. (17)-(24) from Ohmoto (1972), which show that: where T, P, fO 2 and I are the temperature, pressure, oxygen fugacity and ionic strength of the alkali-chloride solution, respectively. To calculate X i we need to know the equilibrium constants for reactions (1)-(6), as well as the activity coefficients for the aqueous species. Equilibrium constants are controlled by both temperature and pressure and in our study, we exported these from the SUPCRTBL thermodynamic dataset (Johnson et al., 1992;Zimmer et al., 2016). Activity coefficients of aqueous species are related to temperature, pressure and ionic strength of the fluids, and were calculated for simple S-bearing alkali-chloride solutions using an extended Debye-Hü ckel equation using the GEM-Selektor software (Wagner et al., 2012). It is important to note that Ohmoto (1972) did not consider the pressure sensitivity of these values in his models and so we have updated Eq. (10) accordingly.
Considering equations (9) and (10), then the d 34 S of each species i: d 34 S i ¼ f ðd 34 S P S ; T; P; fO 2 ; pH; IÞ ð 11Þ In the earlier paper, Ohmoto (1972) considered the absolute values of fO 2 and pH and compared various isotopic fractionations to mineral stability fields. Here we evaluate d 34 S of a cooling hydrothermal fluid with fO 2 tied to common mineral redox buffers: magnetite-hematite (MH), nickel-nickel oxide (NNO) and quartz-fayalite-magnetite (QFM), using the equations compiled by Frost (1991). The advantage of this approach is that we examine the d 34 S evolution of fluids at realistic redox conditions (constrained via petrological methods, Andersen et al., 1993), and predict (rather than assume) the ratio of oxidised to reduced S (SO 4 2À /H 2 S). As with all models we must make several simplifications. Firstly, our models assume that chemical and isotopic equilibria have been established between the aqueous and mineral S species. Ohmoto and Lasaga (1982) showed that temperature, sulfur concentration and, most importantly, pH control sulfate-sulfide equilibrium in hydrothermal fluids. For high temperature hydrothermal systems equilibria is likely achieved irrespective of pH (i.e. at 600°C the time to obtain equilibrium is on the order of minutes, while at 350°C it is hours-weeks). However, when temperatures are low (<300°C ) and pH > 3 equilibria between sulfates and sulfides are slow and unlikely to be attained (i.e. at 250°C equilibrium is achieved in $4 years, while at 200°C it is $100 years, Ohmoto and Lasaga, 1982). Thus, low temperature, nearneutral hydrothermal environments seldom achieve sulfatesulfide equilibrium and hence data from these systems are unlikely to match our equilibrium model predictions (this will be fully discussed in context of our case studies in Section 4).
Further requirements of our model are a constant d 34 S P S value throughout cooling. This assumes the amount of S precipitated by the fluids is negligible compared to the total S left. In accordance with most other studies, we simplify the sulfur system to a single sulfide (S 2-) and sulfate (S 6+ ) couple. Recent work by Binder and Keppler (2011) and Pokrovski et al. (2015) provide experimental evidence for S species (SO 3 and S 3 -) which are not included in our modelling (reactions (1)-(6)). These are excluded because thermochemical data are poorly constrained and fractionation factors unknown (although our model will be easy to adapt once such data become available). Finally, when calculating the d 34 S trajectory of natural hydrothermal fluids, we assume that the fO 2 derived from petrological investigations of the mineral assemblage correspond to the fO 2 of the hydrothermal fluid. This is a reasonable assumption and is supported, for example, by studies of seafloor hydrothermal fluids, where calculated fluid fO 2 can be directly compared to fO 2 from mineral assemblages, and in which the results are in good agreement (Kawasumi and Chiba, 2017).

RESULTS
In this section we describe the key results of our numerical modelling and the S isotope evolution of generic hydrothermal fluids (these models are then applied to porphyry and epithermal deposits, seafloor hydrothermal vents and alkaline igneous bodies in Section 4). Before exploring the d 34 S trajectory of a hydrothermal fluid, it is useful to understand how the position of the mineral redox buffers compare with S speciation predicted by our model. In the Supplementary Information we provide a full description of our S specifications results but here, for brevity, we evaluate S speciation by defining the sulfate-sulfide fence. This is defined when the mole fractions of reduced and oxidised S are equal. Considering the oxidised (S 6+ ) species, the sulfate-sulfide fence is given by: and is shown by the curved lines in Fig. 3 (above the line sulfate dominates and below the line sulfide dominates).
In this model fluid pressure and ionic strength are constant at 1 kbar and 1 M, respectively, and the results clearly show that the position of the sulfate-sulfide fence in T-fO 2 space is governed by fluid pH.
Comparing the location of mineral redox buffers and the sulfate-sulfide fence (Fig. 3), our key observation is that acidic and reduced conditions promote sulfide stability, while alkaline and oxidised conditions promote sulfate. It is important to note that the pH-controlled sulfate-sulfide fences and mineral redox buffers are closely overlapping. Hence, pH and fO 2 are the main controls on the ratio of oxidized to reduced S (SO 4 2À /H 2 S), and, broadly speaking, a change in pH by two units has a similar effect on S speciation to an fO 2 increase of one log unit (e.g. NNO to NNO + 1). Only for a very oxidized fluid at the MH buffer would we observe a sulfate dominated assemblage irrespective of changing pH (between 2 and 8, Fig. 3). Fig. 3. log fO 2 -temperature diagram comparing the location of common mineral redox buffers (from Frost, 1991) with the sulfatesulfide fence of a hydrothermal fluid (Eq. (12)). The sulfate-sulfide fence varies as a function of pH and is defined when the mole fraction of reduced to oxidized S species equals 1. Systems above the line will be sulfate-dominated, while those below the line will be sulfide-dominated. In this model fluid pressure and ionic strength are constant at 1 kbar and 1 M, respectively. Temperature, pH and fO 2 strongly influence the ratio of reduced to oxidised S in a hydrothermal solution and are the main controls on S speciation (and hence isotopic fractionation).
In Fig. 4 we show the S isotope evolution of a range of model fluids. Each model considers a fluid with a fixed bulk S isotope value (d 34 S P S = 0‰) at specified redox buffers (with key variables, pH, P and ionic strength, noted in the annotations). Overall, the models generate a wide spread of isotopic values and mostly fall between the vertical oxidised S line (SO 4 2À /H 2 S = 1) and the horizontal reduced S (SO 4 2À /H 2 S = 0). This projection confirms our conclusion that fluids above the MH buffer are the only ones dominated by a single (oxidised) S phase throughout cooling and follow the vertical SO 4 -dominated line irrespective of changing pH, P and ionic strength (Fig. 4). Fig. 4a-d examines the role of fluid pH and fO 2 on S isotope fractionation. Taking a fixed fO 2 buffer (e.g. NNO) we observe that increasing pH stabilises sulfate (Fig. 3) and causes the sulfide population to tend towards increasingly negative d 34 S (Fig. 4a-d). The sulfate population also tends towards lower d 34 S values with increasing pH and ultimately approach the d 34 S P S value of 0‰ (for NNO this occurs at pH 6, Fig. 4c). The effect of increasing fO 2 on d 34 S is also clear; for example, in Fig. 4a going from NNO to NNO + 1 to NNO + 2 causes the sulfide population to tend toward increasingly negative d 34 S, while the sulfate also decreases. A noticeable feature of our modelled isotope trajectories is that fluids of fixed fO 2 often deflect from reduced (H 2 S dominated) systems at high temperature to more oxidised (SO 4 -dominated) systems at intermediate temperatures, and in several cases switch back to reduced values at lower temperatures. This leads to a U-shaped d 34 S trajectory that swings between the reduced and oxidised S lines (e.g. NNO in Fig. 4b). The d 34 S trajectory is the result of the very subtle differences in curvature of the mineral fO 2 buffers and sulfate-sulfide fence (as a function of temperature, Fig. 3) and illustrates the sensitivity of SO 4 2À /H 2 S to both fluid pH and fO 2 . We also find that reduced, high-pH (>7) systems (e.g. QFM-2 in Fig /H 2 S ratio of the fluid. The bulk S isotope value (d 34 S P S ) is fixed at 0‰ and corresponds to the intersection of the three dashed lines at high temperature (c.f. Fifarek and Rye, 2005;Rye, 2005). The horizontal dashed line indicates cooling of an H 2 S-dominated fluid, the vertical line indicates cooling of an SO 4 -dominated fluid, while the slope of unity indicates equal proportions of SO 4 and H 2 S. Our model calculates the mole fraction of oxidised:reduced S species and S isotope fractionation for a given redox buffer (shown in the legend) at varying temperature, pressure and pH. S isotope trajectories between the horizontal and vertical lines are identified by coloured lines and an erroneous linear interpretation of d 34 S P S is shown by dashed lines in a and i. a-d examine the role of variable pH at fixed pressure. e-g evaluate the role of increasing pressure at fixed pH. i-k show increasing fluid pH at fixed pressure (note pH increases linearly at each 50°C temperature increment). In all models except h, the ionic strength of the fluid (I) is fixed at 1 M. By comparing f and h one can examine the minor impact that changing I has on S isotope fractionation. In general, oxidized fluids (at the MH buffer) plot along the SO 4 -dominated line, while fluids at or below QFM follow the H 2 S-dominated line. Fluids that follow the NNO buffers produce a wide variety of trends (see text for Section 4). have greater abundances of the S 2À species (rather than H 2 S and HSwhich dominate all other models). Since fractionation factors between aqueous S 2À and the precipitated sulfide are larger than for H 2 S and HS - (Ohmoto, 1972;Ohmoto and Rye, 1979;Ohmoto and Lasaga, 1982) reduced, high-pH fluids show sulfide d 34 S above the horizontal line (Fig. 4d).
Sulfur speciation and isotope fractionation are also a function of fluid P and ionic strength (Eq. (9)). Taking a fluid of fixed fO 2 (e.g. NNO + 1 in Fig. 4e-g), we find that changing P has a relatively minor effect, compared to other variables (note that as above, d 34 S P S = 0‰ in all models). Increasing P stabilises sulfate, causing d 34 S sulfate to approach 0‰ and d 34 S sulfide to become increasingly negative (i.e. the U-shaped trajectory is pulled more tightly to the SO 4 -dominated axis at higher pressures). This is driven by the equilibrium constants in our thermodynamic dataset (SUPCRTBL, Johnson et al., 1992;Zimmer et al., 2016) and is consistent with experimental studies (e.g. Binder and Keppler, 2011, who suggest oxidised S dominates in fluids where fO 2 is above NNO and pressures are above 1.5 kbar) and thermodynamic modelling (e.g. Connolly & Galvez, 2018, who predict sulfate, rather than sulfide, species in high pressure slab-derived fluids). Increasing pressures generally favour greater co-ordination numbers, and since S is more highly coordinated with O 2 in sulfate than other species this explains why oxidised S becomes the dominant aqueous species at high pressures. Varying the ionic strength of the fluid (from 1 to 3, Fig. 4f and h) modifies the activity coefficients of the aqueous species but has little impact on the S speciation and isotopic fractionation (in agreement with Ohmoto, 1972). In summary, while a number of variables impact S speciation and isotopes in detail, we find that T-pH-fO 2 are the dominant parameters.
Until now we examine the S isotope evolution of fluids with fixed P-I-pH-fO 2 but in real magmatic-hydrothermal systems these parameters are likely to vary during fluid evolution. While there are many potential fluid evolutionary pathways one could investigate, for brevity we consider a simple scenario where fluid pH increases during cooling. These models are shown in Fig. 4i-k and in each case an initially acidic (pH 1) and high-temperature (600°C) fluid linearly increases pH, to values of 3, 5 and 7, as it cools to 200°C. While this is clearly a simplification of an evolving hydrothermal fluid, these are reasonable model parameters because high-temperature fumarole condensate from typical arc volcanoes have pH of 0-2 (Symonds et al., 1990;Africano and Bernard, 2000), and interaction between these magmatic-hydrothermal fluids and the host rock progressively increase pH (Reed, 1997;Smith et al., 2017). The results ( Fig. 4i-k) show similarities to our earlier findings and confirm that acidic conditions stabilise fluids along the H 2 S-dominated line, while alkaline conditions stabilise fluids along the SO 4 -dominated line. Importantly, the competing effects of changing pH and fO 2 , and their impact on SO 4 2À /H 2 S, can generate a surprising variety of d 34 S trajectories. For example, in Fig. 4j the d 34 S sulfate at NNO is roughly constant between 400 and 200°C, while d 34 S sulfide decreases markedly.
A key observation from Fig. 4 is that magmatichydrothermal fluids with fO 2 between NNO and NNO + 2 seldom define linear d 34 S sulfate -d 34 S sulfide arrays. This emphasises that SO 4 2À /H 2 S is rarely constant in an evolving hydrothermal fluid and has important implications for interpreting natural data. Firstly, individual sulfide or sulfate minerals can show large shifts in d 34 S during cooling (up to $20‰) even when pH is constant and fO 2 fixed to a specific mineral redox buffer. Secondly, the common assumption that SO 4 2À /H 2 S is constant during cooling of a hydrothermal fluid is very rarely achieved. If these diagrams are used to reconstruct d 34 S P S by fitting a linear regression to an array of d 34 S sulfate -d 34 S sulfide data over a limited temperature interval (c.f. Shanks, 2013) then one could easily extract erroneous parameters. As an example, assuming the fluid at NNO + 1 in Fig. 4a, i only precipitated and/or preserved S minerals between 400 and 600°C then a line through these data would be steeply dipping and point to an elevated d 34 S P S of $10‰. This is an erroneous interpretation because, as we stated previously, all models have a fixed d 34 S P S of 0‰. Thus, one of the most significant results from our study is that unless one can demonstrate sulfatesulfide equilibrium and fixed SO 4 2À /H 2 S over a wide range of temperatures (ideally )200°C) then it is inappropriate to use a linear interpretation of the d 34 S sulfate -d 34 S sulfide results to reconstruct d 34 S of the fluid source.

DISCUSSION
Having evaluated the S isotope evolution of generic hydrothermal fluids (above) we now discuss the d 34 S evolution of fluids that form porphyry, epithermal, SMS and alkaline igneous ore deposits. For each deposit type we compiled independent physico-chemical constraints on their formation (summarised in Table 1) and use these to evaluate the typical T, P, pH and redox evolution (or TPpR trajectory) of their ore forming fluids. We use the TPpR trajectory to generate forward models of d 34 S and, rather than focus on single case studies, we evaluate the breadth of d 34 S sulfate -d 34 S sulfide expected for each deposit. We then compare these values to natural data in order to validate our model and identify systems where fluid d 34 S P S departs significantly from magmatic values ($0‰) or where there is disequilibrium between co-existing S-bearing minerals.
Porphyry deposits usually show multiple overlapping stages of alteration and mineralization (Sillitoe, 2010). This starts with a phase of potassic alteration (biotite, Kfeldspar, magnetite assemblages) formed at high temperatures (700-350°C) and near-neutral pH, which is usually associated with ore deposition. The potassic phase is commonly followed and overprinted by phyllic alteration (quartz, muscovite (sericite), pyrite, chlorite assemblages) indicative of lower temperatures ($350°C) and moderately acidic pH. The latest stages are characterized by advanced argillic alteration (quartz, kaolinite, alunite pyrophyllite, dickite assemblages) which form at low temperatures ($300°C) and highly acidic conditions. Ultimately, this temporal evolution reflects the progressive cooling of the magmatic fluids and the fact that acid volatiles (especially SO 2 ) are fully associated at high temperatures but disproportionate at lower temperature to form sulfuric acid and H 2 S (Giggenbach, 1997).
Our models (Fig. 5a-c) replicate this TPpR evolution. We consider fluid pressures of 1-2 kbar (Table 1) which are appropriate for the typical depths of porphyry mineralisation (above) and also encompass upper crustal magma storage conditions (e.g. 1.8-2 kbar at Pinatubo, Scaillet and Evans, 1999). Temperatures span 800-350°C and are consistent with fluid inclusion and isotope geothermometers from porphyry deposits (Roedder, 1971;Eastoe, 1978;Landtwing et al., 2005;Rye, 2005), as well as experimental constraints on Pinatubo magmas (Scaillet and Evans, 1999). Redox conditions fall between NNO and NNO + 2 in agreement with constraints from petrology and gas geochemistry (Gerlach and Casadevall, 1986;Symonds et al., 1990;Scaillet and Evans, 1999;Audétat and Pettke, 2006). We assign pH intervals representative of near-neutral (5-6), moderately acidic (4-5) and highly acidic conditions (2-4), which we tentatively link to the three main phases of alteration (potassic, phyllic and advanced argillic, Fig. 5a-c). As many variables can potentially alter fluid pH these intervals are considered indicative rather than absolute and are used to garner a broad understanding of the TPpR and d 34 S evolution of porphyry fluids. Finally, our porphyry models (Fig. 5a-c) are setup with a bulk S isotope value (d 34 S P S ) of 5‰. We selected this value because numerous studies of volcanic rocks and melt inclusions from calc-alkaline volcanoes suggest elevated d 34 S P S (e.g. Pinatubo = 5-6.5‰, McKibben et al., 1996; El Chichó n = 5.8‰, Luhr et al., 1984 and Krakatoa = 1.5-4‰, Mandeville et al., 1998) and this allows us to test whether genetically-linked porphyry fluids also require enriched S sources.
Model results (Fig. 5a-c) clearly show the dual roles of pH and fO 2 in controlling S speciation and isotope fractionation. While there are a wide range of model outputs that fall between the reduced (SO 4 2À /H 2 S = 0) and oxidised (SO 4 2À /H 2 S=1) lines we predict that fluid evolution (increasing acidity Fig. 5a-c) would result in and initially vertical d 34 S array (high SO 4 2À /H 2 S) being driven to the upper right (low SO 4 2À /H 2 S) region of the plot (i.e. decreasing d 34 -S sulfide and increasing d 34 S sulfate ). In Fig. 5d, we overlay d 34 S data from volcanic and plutonic samples typical of calcalkaline porphyry environments. Most sulfate-sulfide pairs suggest temperatures of 800-300°C and we note that those that do show temperatures <300°C (e.g. Far Southeast, Hedenquist et al., 2017) were usually sampled from the margins of the porphyry bodies. Comparing the d 34 S data Table 1 Summary of the ore deposits considered in this study and the physico-chemical parameters (temperature-pressure-pH-fO 2 : TPpR) used to calculate S speciation and isotope fractionation in our model hydrothermal fluids. Note abbreviation SMS: seafloor massive sulfide. See text and references for further information on our chosen physico-chemical parameters.

Ore deposit Tectonic setting Magmatic composition
Modelled physico-chemical parameters (TPpR) T (°C) P (kbar) pH log fO 2 Key references Porphyry Subduction zone Calc-alkaline 800-350 1-2 2-6 NNO to NNO + 2 Luhr et al., 1984;Gerlach and Casadevall, 1986;Reed, 1997;Scaillet and Evans, 1999;Rye, 2005;Audétat and Pettke, 2006;Sillitoe, 2010 Epithermal Subduction zone Calc-alkaline 350-150 0.2 0-2 NNO to NNO + 3 Casadevall & Ohmoto, 1977;White and Hedenquist, 1990;Stoffregen, 1987;Fifarek and Rye, 2005 Braunger et al., 2018;Friel and Ulmer, 1974;Markl and Baumgartner,2002;Hettmann et al., 2014;Graser and Markl, 2008;Schö nenberger and Markl, 2008;Mitchell and Krouse, 1975 Carbonatite 800-200 1-2 5-7 QFM to QFM + 4 for different complexes, shows some display linear d 34 S arrays (e.g. El Salvador and Gaspe) while others (e.g. Butte and Frieda) are more scattered. Our modelled TPpR trajectories, with an enriched S source (d 34 S P S = 5‰‰), capture the isotopic variation in the natural samples (Fig. 5d). For comparison, a typical MORB like d 34 S P S of 0‰ is shown in Fig. 5e and is notable for the data points above the reduced S line. These outliers cannot be explained by the typical TPpR trajectory of a porphyry fluid with d 34 S P S of 0‰. Thus, an elevated d 34 S P S (Fig. 5d) is required to explain the S isotope systematics of several porphyry systems and supports observations from calc-alkaline volcanic rocks which suggest high d 34 S P S in subduction zone settings (Luhr et al., 1984;McKibben et al., 1996;Mandeville et al., 1998). Elevated d 34 S P S may be linked to a variety of processes including slab contributions to the mantle (De Hoog et al., 2001); underplating and degassing phenomena (McKibben et al., 1996) and crustal assimilation (Field et al., 2005). Not all deposits necessarily require d 34 S P S > 0‰, for example, Panguna, which shows a constant low SO 4 2À /H 2 S over a wide temperature range (Eastoe, 1983). Again, several S isotope studies of calc-alkaline magmas agree with this, for example tephra and scoria from Mt Mazama system suggest d 34 S P S values of $0‰ (Mandeville et al., 2009). We emphasise that porphyry fluids, like arc magmas, show a range of d 34 S P S (mostly between 0 and 5‰). Although our models do not allow us to assess whether the enriched d 34 S P S is a genuine feature of arc mantle sources or reflects similar magmatic processing, it is worth noting that primitive melt inclusions from the Lesser Antilles Arc (Bouvier et al., 2008) and metamorphic sulfides from exhumed subduction zones (Walters et al., 2019) have been used to estimate a wide range of d 34 S values for slab fluids /H 2 S ratio of the fluid. In our models (a-c) data symbols correspond to the redox buffer, annotations describe fixed model parameters and a bulk S isotope value (d 34 S P S ) of 5‰ was. In a-c, decreasing fO 2 and pH stabilise reduced S species and drive d 34 S values to the upper right portion of the plot. Arrows show the average d 34 S trajectory for the NNO buffer changing with increasing acidity. In d and e, model results (shown by the shaded regions) are compared to observations from natural systems. Data symbols correspond to the sulfate-sulfide mineral pair while data colours reflect the specific complex. Note that correct usage of d 34 S sulfate -d 34 S sulfide diagrams requires all sulfide mineral d 34 S to be adjusted to aqueous H 2 S d 34 S values (Rye, 2005) and for completeness we updated these values using fractionation factors of Ohmoto and Rye (1979) and Ohmoto and Lasaga (1982). In d d 34 S P S = 5‰, while in e d 34 S P S = 0‰. The former provides a better fit to the natural data and suggests that many calc-alkaline hydrothermal systems are enriched in d 34 S due to crustal assimilation, degassing processes and/or a slab-influenced mantle source. In d and e coloured dashed lines show our interpretation of SO 4 2À /H 2 S. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) (-11 to +8‰). A final point regarding S sources is that because porphyry ore deposits (Fig. 5d) exhibit variations in d 34 S P S this implies that the vital step in ore formation is not the magmatic and deep mantle processes which modify d 34 S P S but rather the magmatic-hydrothermal processes which locally concentrate S and metal-rich fluids.
In previous reviews and compilations of porphyry d 34 S (Rye, 2005;Seal, 2006;Shanks, 2013) it was suggested that linear interpolations of d 34 S sulfate -d 34 S sulfide arrays could be used to reconstruct d 34 S P S . However, our new modelling (Section 3, Fig. 4) emphasises that these approaches often yield unreliable results, and that the slope of natural d 34 S sulfate -d 34 S sulfide arrays should only be used to evaluate redox evolution (i.e. SO 4 2À /H 2 S). Porphyry systems which define linear arrays over a temperature range >200°C are highlighted in Fig. 5d,e. For the two volcanic systems, Pinatubo and El Chichó n, our interpretation of SO 4 2À / H 2 S is in broad agreement with their calculated magma redox (i.e. Pinatubo magmas are known to be more oxidised than those of El Chichó n, Scaillet et al., 1998). This is consistent with d 34 S values being controlled by SO 4 2À /H 2 S and thus closely linked to the oxidation state of the parent magmas. For porphyry bodies between NNO and NNO + 1 (Fig. 5a-c) we predict that if isotopic measurements were made over a range of temperatures, one would observe a transition in d 34 S values from an initially vertical (SO 4 -dominated) array to a more horizontal (H 2 S-dominated) array at lower temperatures reflecting increasing acidification of the fluids (Fig. 5a-c). The important point is that the oxidised (and mineralized) potassic core of the porphyry body should display the lowest d 34 S sulfide values. Although few systems have enough paired sulfate-sulfide data to examine this thoroughly (Fig. 5d), relationships between d 34 S sulfide and alteration that match our model predictions have been reported at alkalic porphyry deposits (e.g. Wilson et al., 2007). Thus, our models ( Fig. 5a-c) explain S isotope zonation in porphyry systems and highlight the potential usefulness in fingerprinting different phases of alteration and mineralization.
Our interpretations of SO 4 2À /H 2 S ratios (Fig. 5d, e) also provide insights into ore forming processes and it is notable that no single SO 4 2À /H 2 S typifies economic mineralisation. In a long-lived magmatic-hydrothermal environment, typical of a porphyry system (Sillitoe, 2010), we would expect significant spatial and temporal variations in SO 4 2À /H 2 S. We have shown that S speciation can change dramatically during the TPpR evolution of a single fluid pulse (Fig. 4) and if multiple fluid injections, with slightly different TPpR trajectories, occurred through time then this would yield a scattered d 34 S sulfate -d 34 S sulfide distribution (Fig. 5a-c). It is therefore very surprising that several porphyry systems display linear d 34 S sulfate -d 34 S sulfide arrays over a wide range of temperatures (Fig. 5d, e). Although relatively rare, linear d 34 S sulfate -d 34 S sulfide arrays were generated in a few of our simulations (e.g. NNO + 1 in Fig. 4a, i) and occur when the changing fluid TPpR parameters balance out any variations in SO 4 2À /H 2 S. While it is possible that the linear trends observed in some porphyries reflect steady cooling of a single fluid with TPpR parameters that effectively held SO 4 2À /H 2 S constant, we believe that most magmatic-hydrothermal environments would be subject to repeated injections of new melt and fluid (e.g. Taran et al., 2013), which would cause large changes in TPpR and hence significant variations in SO 4 2À /H 2 S and S isotopes. An alternative explanation is that linear d 34 S trends record rapid precipitation of sulfide and sulfate minerals from a hydrothermal fluid with approximately fixed SO 4 2À /H 2 S. For example, if the NNO + 1 hydrothermal fluid shown in Fig. 4a, i was at equilibrium over a temperature range of $600-300°C, and was suddenly decompressed then this would trigger S mineral precipitation and preserve a linear d 34 S sulfate -d 34 S sulfide array. Importantly, numerous authors have suggested that abrupt physical perturbations (i.e. decompression caused by faulting or magmatic fluid injections) can initiate ore depositional events (Weis et al., 2012;Richards, 2018). Our modelling supports these scenarios because fixed SO 4 2À /H 2 S is unlikely to be maintained in a dynamic magmatic-hydrothermal system (although detailed textural analysis and S-isotope microanalysis, e.g. Peterson and Mavrogenes, 2014, complemented by TPpR modelling, would be the best approach for testing this further, Section 4.5).

Epithermal deposits
Epithermal deposits are also generated by metal-rich hydrothermal fluids circulating above magmatic bodies, but, in contrast to porphyry deposits, they form at shallower depths (generally <1 km) and occur in a wider variety of tectonic settings (Sillitoe and Hedenquist, 2003). Epithermal deposits are often divided into acidic, oxidized 'high sulfidation' systems and near-neutral, reduced 'low sulfidation' systems (White and Hedenquist, 1990). These divisions reflect differences in fluid chemistry, mineralogy and alteration, which are ultimately tied to tectonic setting, magmatic source and host rock composition (Table 1, Sillitoe and Hedenquist, 2003). High sulfidation deposits are usually found in zones of advanced argillic alteration (mentioned above) and are thought to represent the uppermost section of a calc-alkaline magmatichydrothermal system, and thus genetically related to deeper porphyry deposits (Hedenquist et al., 1998;Rye, 2005;Sillitoe, 2010). During formation of an epithermal deposit, fluid compositions and host rock interactions may change spatially and temporally (Smith et al., 2017), such that a low pH high sulfidation system may evolve to a nearneutral low sulfidation system (e.g. Summitville, Bethke et al., 2005).
Our d 34 S models investigate the different TPpR trajectory of high and low sulfidation deposits (Table 1). High sulfidation models (Fig. 6a, b) consider acidic conditions (pH 2, e.g. Stoffregen, 1987;Fifarek and Rye, 2005) while low sulfidation models (Fig. 6c, d) consider near-neutral conditions (pH 4-6, e.g. Casadevall and Ohmoto, 1977). Pressure is set at 200 bars consistent with shallow depths of formation and constraints from fluid inclusions (Casadevall and Ohmoto, 1977;White and Hedenquist, 1990). Temperatures range from 350 to 150°C which are also consistent with fluid inclusion and isotopic evidence Fig. 6. Models of equilibrium fractionation of sulfur isotopes between sulfide and sulfate in high (a, b) and low sulfidation (c, d) epithermal deposits compared to natural samples (e, f). Solid diagonal lines define isotherms while dashed lines define the SO 4 2À /H 2 S ratio of the fluid. In a-d, symbols correspond to the redox buffer and annotations describe fixed model parameters. Arrows in a and d show the dominant d 34 S trajectory predicted by our models. In e and f, model results (shown by the shaded regions) are compared to observations from natural systems. Data symbols correspond to the sulfate-sulfide mineral pair while data colours reflect the specific complex. It is important to note that all sulfide mineral d 34 S was adjusted to aqueous H 2 S d 34 S using fractionation factors of Ohmoto and Rye (1979) and Ohmoto and Lasaga (1982). d 34 S values from Creede caldera (Rye et al., 1988) and Crofoot-Lewis (CFL, Ebert and Rye, 1997) are shown by the grey dashed rectangles (note these data indicate sulfate-sulfide d 34 S disequilibrium, see text for Section 4). White stars show averaged d 34 S values of jarosite (a sulfate mineral that often forms from supergene alteration of hydrothermal sulfides, Stoffregen, 1993, Rye, 2005 from Creede, Crofoot-Lewis and Paradise Peak (PP, Rye and Alpers, 1997). In e and f, d 34 S P S of the shaded model is 5‰. In e the outline of the high sulfidation (low pH) models with d 34 S P S = 0‰ is shown by the dashed red line. For high sulfidation epithermal deposits (e), our models and data are in good agreement and suggest that magmatic S sources (0-5‰) typify most of these systems. Data from low sulfidation and transitional systems (i.e. those that show characteristics of both low and high sulfidation deposits) are shown in f and display a wide range of d 34 S values indicative of disequilibrium processes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) (e.g. Bethke et al., 2005;Wilkinson, 2001). Detailed information on the redox state of epithermal fluids is generally lacking, although it is reasonable to assume similar, if not greater, fO 2 than the porphyry fluids modelled previously (Table 1), and so we extend our analysis to NNO + 3 (c.f. Pokrovski et al., 2015). To maintain consistency with the porphyry models a d 34 S P S of 5‰ is used in Fig. 6a-d. As with the previous example we also evaluate d 34 S P S of 0‰ in the model-data comparisons (Fig. 6e).
Our results show that all low pH models (Fig. 6a, b), with the exception of those at NNO + 3, plot close to the reduced S line (SO 4 2À /H 2 S = 0). Near-neutral pH models (Fig. 6c, d) exhibit a wider range in d 34 S but are mostly focused on the oxidised S line (SO 4 2À /H 2 S = 1). In Fig. 6e and f data from high and low sulfidation epithermal systems are overlain on their respective TPpR models. Samples from Summitville and Pueblo Viejo are transitional between high and low sulfidation systems and included with the latter. At Summitville the hydrothermal system evolved from an acidic high sulfidation system to a near-neutral low sulfidation system (Bethke et al., 2005), while at Pueblo Viejo fluid pH was apparently 2-3 (Muntean et al., 1990) and so less acidic than typical high sulfidation systems we modelled.
Focusing first on the high sulfidation systems (Fig. 6e), we see good overlap between the TPpR models and the natural samples. Most data suggest fluid temperatures between 300 and 200°C, substantially lower temperatures than porphyry systems (Fig. 5d, e), and very typical of fluid inclusion trapping temperatures measured in epithermal systems (300-100°C, Wilkinson, 2001). Magmatic sources (d 34 S P S of 0-5‰) provide a good fit to the data. The only exception is Rodalquilar where one data point necessitates a d 34 S P greater than 5‰. In the case of Rodalquilar, Sr isotopes of sulfate minerals show that this non-magmatic S was most likely derived from seawater or leached from nearby marine sedimentary rocks, before being mixed with magmatic S in the deposit . Given that epithermal deposits form near surface (Hedenquist et al., 1998;Sillitoe, 2010) one might expect greater influx of external (non-magmatic) S, yet our models are able to confirm that magmatic S dominates virtually all high sulfidation epithermal systems.
Data from low sulfidation and transitional systems are plotted in Fig. 6f and show slightly lower d 34 S sulfide than high sulfidation (low pH) systems. Few low sulfidation systems show evidence for sulfate-sulfide equilibrium and so the limited data on Fig. 6f reflect a lack of d 34 S sulfate -d 34 S sulfide pairs rather than a paucity of case studies. Well-known examples showing isotope disequilibrium are Crofoot-Lewis (Ebert and Rye, 1997), which would imply unrealistically high temperatures (Fig. 6f), and Creede (Rye et al., 1988;Rye, 2005), which shows an incredibly wide spread of d 34 S sulfate and clear disequilibrium with co-existing sulfides. Note that in both these cases we plot the range of values in Fig. 6f because mineral pairs are lacking. The low temperatures (<250°C) and high pH (3-7) that typify low sulfidation epithermal fluids mean that timescales of sulfate-sulfide equilibrium are several years (Ohmoto and Lasaga, 1982) and unlikely to be achieved in these dynamic hydrothermal systems (Rye, 2005). Thus, while there is overlap between the data and a few of our models (and indeed the limited data are shifted to lower d 34 -S sulfide consistent with our predictions, Fig. 6f) we are cautious ascribing these values to the d 34 S sulfate -d 34 S sulfide of the fluid.
Our investigation of epithermal systems reveals some of the limitations of our models. Although TPpR trajectories predict that S isotopes may be used to distinguish high and low sulfidation epithermal deposits (Fig. 6a-d), the timescales required for S mineral equilibrium in low sulfidation systems are much longer than the timescales of fluid advection and mineral precipitation. Thus, sulfate-sulfide disequilibrium is the norm in low sulfidation systems. Another complicating factor in near-surface epithermal environments are supergene processes such as oxidative weathering. Supergene oxidation of primary sulfide minerals leads to the formation of hydrous sulphates, such as jarosite, which have been identified in various weathered epithermal deposits (e.g. Rye and Alpers, 1997;Rye, 2005). Bulk-rock comparisons of d 34 S sulfate and d 34 S sulfide from supergene environments clearly indicate sulfatesulfide disequilibrium (i.e. unrealistically high temperatures shown by white starts in Fig. 6f) and emphasise the need for careful textural and mineralogical investigations prior to isotope modelling (Section 4.5). Despite these complicating factors, we underscore that for high sulfidation deposits (that are unmodified by supergene processes, and where equilibrium is likely obtained in hours-weeks) a robust model-data comparison can be made. In this case our TPpR trajectories provide a good fit to the data (Fig. 6e) and clearly indicate a dominance of magmatic S sources.

Seafloor massive sulfide (SMS) deposits
Modern SMS deposits and their geological equivalents volcanogenic massive sulphides (VMS) form when convecting hydrothermal fluids are exhaled along seafloor vents (Herzig and Hannington, 1995). The ore-forming fluid is usually seawater that has undergone reaction with the host rock (Woodruff and Shanks, 1988;Shanks, 2001), although in a few modern hydrothermal fields there is evidence for a significant contribution of juvenile magmatic fluids and volatiles (Gamo et al., 1997).
Hydrothermal fluids associated with active seafloor vents can be sampled directly (e.g. Shanks, 2001;Ding et al., 2005;Kawasumi and Chiba, 2017) and so we focus our modelling on modern SMS rather than ancient VMS deposits (because the TPpR trajectory of the former can be more accurately constrained). While SMS deposits are found in a variety of settings (Shanks, 2001;Piercey, 2011) we focus our modelling on two endmembers: (1) Typical mid-ocean ridge (MOR) sulfide chimney deposits hosted in basaltic rocks. In this case, hydrothermal fluids are dominated by evolved seawater (Woodruff and Shanks, 1988;Shanks et al., 1995;Shanks, 2001) and are only slightly acidic (pH 4-5, Ding et al., 2005;Seyfried et al., 2011).
Our modelled TPpR trajectories (Fig. 7a-c) evaluate how the above variations in fluid pH may govern the d 34 S systematics of seafloor vent fluids (at temperatures of 350-200°C, fO 2 between QFM to NNO and fixed pressures of 300 bars, see Table 1 for key references). A typical MORB like d 34 S P S of 0‰ is used in modelling the TPpR trajectory and the results show that vent fluids with a pH between 4 and 5 (typical of most MOR systems) should generate a wide spread of d 34 S values depending on the fO 2 of the fluid (Fig. 7a). Low pH models (typical of back-arc calcalkaline settings, Fig. 7b, c) cluster around the reduced S line (SO 4 2À /H 2 S = 0) and produce similar d 34 S results to high sulfidation epithermal deposits (Fig. 6a, b).
When the d 34 S values of S minerals from modern SMS deposits are overlain on our TPpR trajectories (Fig. 7d, e) it is clear that only a couple of sites show any agreement. Most data plot as vertical arrays with d 34 S sulfate between 19 and 22‰ (equivalent to the d 34 S of modern seawater sulfate, Rees et al., 1978;Tostevin et al., 2014). Changing d 34 S P S does little to improve the model-data fit, and several datasets from the Mid-Atlantic Ridge and the Valu Fa back-arc Ridge suggest sulfate-sulfide formation temperatures that are unrealistically high (600-400°C, Chiba et al., 1998;Herzig et al., 1998). Thus, our equilibrium TPpR trajectories (Fig. 7e) fail to predict the d 34 S of natural seafloor hydrothermal settings and provide strong evidence that S phases were not formed in equilibrium (i.e. sulfates and sulfides were precipitated from different fluids and did not re-equilibrate, leading to bulk rock disequilibrium, Seal et al., 2000). Fig. 7. Models of equilibrium fractionation of sulfur isotopes between sulfide and sulfate in seafloor massive sulfide (SMS) deposits (a-c) compared to data from active seafloor vents (d, e). Solid diagonal lines define isotherms while dashed lines define the SO 4 2À /H 2 S ratio of the fluid. In our models (a-c) a bulk S isotope value (d 34 S P S ) of 0‰ was used and annotations above each plot describe fixed model parameters. Arrow in c shows the dominant d 34 S trajectory predicted by our models. In d and e, model results (shown by the shaded regions) are overlain by observations from natural systems. Data symbols correspond to the sulfate-sulfide mineral pair while data colours reflect the specific vent site (see key, lower right). The natural data show a poor fit to the models and indicate disequilibrium between the sulfate-sulfide mineral pairs. Sulfates mostly reflect the d 34 S of contemporaneous seawater (labelled in e). Sulfides that cluster around 0-5‰ are mainly derived from igneous S (0‰) with a small contribution from reduced seawater S (vertical black bar). At a few sites greater contributions from sulfate reduction or biogenic S are required to explain high or low d 34 S sulfide values (arrows). Anomalous low sulfide d 34 S at Susu Knolls and Hine Hina likely reflect a direct input of magmatic volatiles (see text for Section 4). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) These findings agree with most previous studies of seafloor hydrothermal systems which have consistently highlighted disequilibria and distinct origins for sulfate and sulfide minerals (Shanks, 2001). At seafloor pressures sulfate exhibits retrograde solubility and so downwelling seawater precipitates anhydrite at temperatures above 150°C (Shanks et al.1981;Sleep, 1991). Since most d 34 S sulfate measurements overlap the composition of seawater (Fig. 7e), this supports near quantitative stripping of seawater sulfate during hydrothermal convection. Sulfide minerals may be derived from a variety of sources, including: (1) S leached from the volcanic host rocks ($0‰); (2) S from reduced seawater (with d 34 S typically 15-19‰ lower than modern seawater, $21‰, Huston, 1997;Rees et al., 1978;Tostevin et al., 2014); and (3) S leached from biogenic materials in seafloor sediments (which averages -25‰ but ranges from 10 to -55‰, Canfield & Farquhar, 2009). Most sulfides have d 34 S between 1 and 5‰ (Fig. 7e), which is usually taken as evidence that ore-forming fluids are dominated by a volcanic host rock d 34 S signature, with a small contribution from seawater sulfate reduced to sulfide in the reaction zone or chimney environment (Woodruff and Shanks, 1988;Janecky and Shanks, 1988;Alt, 1995;Alt and Shanks, 2011). In general, sulfides with anomalously high d 34 S (>5‰, Fig. 5) require greater contributions from reduced seawater S and/or reduction of early-formed anhydrite (e.g. TAG, Chiba et al., 1998), while those with anomalously low d 34 S (<0‰) require biogenic S sourced from the overlying sediment (e.g. Guaymas, Peter and Shanks, 1992). It is important to recognize that d 34 S on its own does not provide a unique fingerprint of these sources (i.e. the 1-5‰ range observed in most sulfides could be explained via a mixture of reduced seawater S and biogenic S sources, rather than a predominance of volcanic host rock S) and we emphasise that textural studies and multiple S isotope investigations are usually required to unpick these processes (discussed in Section 4.5).
The only data which overlap our TPpR trajectories are from Hine Hina and Susu Knolls (Fig. 7e) which are both located in back-arc settings and hosted in calc-alkaline rocks. Their S isotope values place them in the low pH region of our d 34 S sulfate -d 34 S sulfide plots and is consistent with their acidic vent chemistries (pH 1-3, Herzig et al., 1998;Kim et al., 2004). Although sulfate-sulfide equilibrium is unlikely (Herzig et al., 1998), their isotopic values suggest a d 34 S P S value closer to $0‰ (Fig. 7e) and hence a much greater input of acidic magmatic volatiles at these sites (which also supports interpretations from vent fluid chemistry, Reeves et al., 2011).
S isotopes of SMS deposits are a valuable counterpoint to the porphyry and epithermal systems investigated in Sections 4.1 and 4.2. While the latter are dominated by magmatic S ($0‰), seafloor hydrothermal fluids incorporate a wide range of isotopically variable S reservoirs and, for the vast majority of systems, the direct contribution from magmatic S is extremely limited. Sulfate-sulfide disequilibrium in SMS systems reflects their wide range of S sources as well as the strong thermal and chemical gradients which typify these reaction zones (Lowell et al., 1995), and our equilibrium models of a cooling hydrothermal fluid with d 34 S P S of 0‰ clearly do not capture these complexities. Despite this, the results from SMS deposits do provide a useful example of how our forward models can be used to flag bulk rock d 34 S disequilibrium in magmatichydrothermal systems.
Finally, though fossil VMS deposits were not the focus of our modelling it is important to note that similarities are often drawn between these and active SMS environments (Herzig and Hannington, 1995). VMS span a much wider range of d 34 S than modern SMS. Individual VMS deposits show extreme isotopic variability (e.g. sulfides from Dry Creek VMS span 71‰, Slack et al., 2019) and there are well-documented secular variations in VMS d 34 S throughout the geological record related to the evolution of the atmosphere, biosphere and hydrosphere (e.g. Huston, 1997;Farquhar et al., 2010, -50 to +50‰ in sulfides and 5-50‰ for sulfates). VMS d 34 S primarily captures the long-term evolution of exogenous S reservoirs (e.g. seawater d 34 S) rather than the evolution of their contemporaneous magmatic system. Thus, model-data mismatch observed in contemporary SMS (Fig. 7e) is also expected for fossil VMS.

Alkaline igneous deposits
We now consider the d 34 S systematics of alkaline magmatic fluids which, by comparison to well-studied examples above, have few detailed investigations of their TPpR evolution and less developed models of ore-forming processes. Alkaline magmatic fluids are associated with carbonatites and alkaline silicate intrusions (Verplanck and Van Gosen, 2011;Jones et al., 2013), and represent some of the most chemically evolved magmatic fluids on Earth. They often exhibit extreme enrichments in REE and HFSE which are increasingly sought for high technology applications (Zhou et al., 2017). The mineralized bodies associated with alkaline igneous systems are located within the vicinity of the magmatic intrusion, usually within the roof zone (Marks and Markl, 2015; or within late-stage dykes radiating from the complex (Le Bas, 1977). These deposits host REE and HFSE in a variety of minerals including silicates, oxides, carbonates and phosphates.
Our review of the TPpR parameters that typify alkaline magmatic fluids (Table 1) suggests that mineralisation can occur over a wide range of temperatures (800-200°C) and oxidation states (from QFM-4 to QFM + 3, Marks and Markl, 2015;Braunger et al., 2018Braunger et al., , 2020. These parameters vary significantly from complex to complex and so our models compare two endmember redox trends for alkaline fluids (a reduced alkaline silicate and an oxidised carbonatite, Table 1).
Reduced models (Fig. 8a, b) are based on alkaline magmatism in the Gardar Province (Greenland), and in particular the agpaitic nepheline syenites found in the roof of the Ilímaussaq complex (Marks and Markl, 2015). The TPpR trajectory of these REE-rich fluids have been constrained by mineralogy, phase equilibria and fluid inclusions (Table 1) and are notable for their high pH (>7, Markl and Baumgartner, 2002) and extremely low fO 2 (QFM to QFM-8, Markl et al., 2001;Hettmann et al., 2014). Our models (Fig. 8a) predict that d 34 S values for these highly reduced fluids would closely follow the horizontal sulfidedominated line (SO 4 2À /H 2 S = 0). It is important to note that their high pH and reduced nature means S 2is the dominant S phase and would lead to elevated d 34 S in the precipitated sulfide minerals (see Section 4 and Fig. 4d). A key insight from previous studies of the Gardar province is that when late-stage alkaline fluids are injected into country rocks, at the roof and margins of the complex, they evolve to more oxidised conditions (around the MH buffer) at temperatures of 300-200°C (Graser and Markl, 2008). This shift in redox reflects the fact that fluids within the intrusion are initially buffered by reduced alkaline rocks but once injected into the surrounding country rock they are mixed with externally derived oxidised fluids (Hutchison et al., 2019). Thus, a realistic evolutionary trajectory for initially reduced Ilímaussaq-type fluids would involve a rapid change in redox at low temperatures and generates a clockwise rotation in d 34 S (combining Fig. 8a and b).
Oxidised models (Fig. 8c) are based on the Kaiserstuhl Volcanic Complex (SW Germany), which comprises a variety of alkaline silicate rocks and carbonatites (Keller, 1981). Detailed petrological investigations of Braunger et al. (2018) have recently developed a thorough understanding of the conditions of magma storage (Table 1) and, in contrast to the examples above, have recognized relatively oxidised fO 2 values between QFM + 1 to QFM + 3. Unfortunately, pH values for oxidised carbonatitic fluids are poorly constrained at present and so we examine a Fig. 8. Models of equilibrium fractionation of sulfur isotopes between sulfide and sulfate in alkaline igneous deposits (a-c) compared to data from natural systems (d-f). Solid diagonal lines define isotherms while dashed lines define the SO 4 2À /H 2 S ratio of the fluid. In our models (a-c) a bulk S isotope value (d 34 S P S ) of 0‰ was used and annotations above each plot describe fixed model parameters. a and b consider the cooling and oxidation of an initially reduced alkaline silicate fluid, while c considers the cooling of a carbonatititc fluid (see text for Section 4).
Arrow in a-c shows the dominant d 34 S trajectory predicted in each model. In d we compare the d 34 S of sulfate-sulfide mineral pairs from alkaline rocks to the fields defined by porphyry and epithermal systems (grey shaded regions). Note data symbols correspond to the mineral pair while data colours reflect the specific complex (see key at bottom of plot). In e model results (shown by the coloured shaded regions) are overlain by observations from natural systems. Data in e have been corrected to a d 34 S P S of 0‰ using values of d 34 S P S from magmatic cumulate (this allows a more direct comparison of the data and models since both have d 34 S P S = 0‰). In f data symbols correspond to mineralogy of the main REE-bearing phase. Data from alkaline and carbonatitic fluids suggest a wide range of temperature and redox values, and we envisage a variety of TPpR trajectories (shown schematically by the arrows in e and f). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) similar range of values to the reduced case (note that these values are consistent with the pH used in previous modelling, c.f. Drü ppel et al., 2006). The oxidised model (Fig. 8c) shows that most d 34 S values plot along the sulfate dominated line (SO 4 2À /H 2 S = 1), only QFM models record a transition from reduced to oxidised conditions at temperatures above 600°C.
S isotope data from alkaline igneous systems are shown in Fig. 8d. Virtually all S isotopic measurements have been made on carbonatites, the only exceptions are the alkaline silicate rocks of the Gardar (Ilímaussaq, Motzfeldt and Ivigtut, Hutchison et al., 2019) and the unusual alkaline-carbonatite kimberlites of the Udachnaya-East pipe (Siberia, Kitayama et al., 2017). An important caveat is that most Gardar fluids were sufficiently reduced (at temperatures above 300°C) such that no sulfates were present. To plot isotopic data for these intrusions in Fig. 8d we used sulfide mineral pairs to estimate temperatures and calculate the theoretical d 34 S sulfate that would be in equilibrium with these minerals. Since our models make predictions about the d 34 S of sulfate-sulfide equilibrium pairs irrespective of mineral stability, this approach allows us to make comparisons between the models and data (Fig. 8e). Note this approach should only be applied when there is robust petrographic evidence for a completely reduced or oxidised S mineral assemblage (such as the Gardar, Hettmann et al., 2012;Hutchison et al., 2019).
d 34 S data from alkaline igneous systems show a remarkably wide range of values (Fig. 8d). Alkaline silicate fluids extend along the reduced S line (SO 4 2À /H 2 S = 0), while carbonatite fluids tend to show lower d 34 S indicative of more oxidised conditions. Compared to porphyry and epithermal systems, d 34 S data from alkaline systems clearly show that mineralisation take places over a much wider range of temperatures, and that redox conditions of most carbonatitic fluids are more oxidised and sulfate dominated (Fig. 8d) than fluids circulating calc-alkaline environments.
Another important feature of alkaline d 34 S data ( Fig. 8d) is that they appear to show both positive d 34 S P S (e.g. Gardar, Ilímaussaq and Ivigtut) and negative d 34 S P S (e.g. Brazil, Saltire and Tapira). Hutchison et al. (2019) showed that these variations in d 34 S P S are unrelated to changes in melt fractionation and fluid evolution and appear to represent genuine characteristics of the magma source. Using the d 34 S P S values calculated by Hutchison et al. (2019), mostly from magmatic cumulate, we correct all data in Fig. 8e to a d 34 S P S of 0‰ so that we can more directly compare our TPpR trajectories to the data. We find that our reduced and oxidised forward models encompass the wide d 34 S range of natural samples. Alkaline silicates closely match the modelled d 34 S trajectory for Ilímaussaqtype fluids, while most carbonatites have d 34 S trajectories between the reduced and oxidised endmembers (with an intermediate example shown by the blue arrow in Fig. 8e).
It is clear from Fig. 8e that d 34 S data from individual alkaline igneous systems are sparse and this makes it challenging to unravel their T and redox evolution. We find that a useful alternative is to consider the key ore mineral at each locality, and three obvious groupings emerged: (1) silicates; (2) oxides and (3) carbonates and phosphates. Labelling the data in this way (Fig. 8f) reveals a relationship between S isotopes, redox (i.e. SO 4 2À /H 2 S) and ore mineralogy. Under very reduced conditions silicate minerals (e.g. eudialyte and steenstrupine) form the main REE-bearing minerals, while at high fO 2 carbonate and phosphate minerals prevail (e.g. monazite, bastnäsite and apatite). In between these two extremes, oxide minerals (e.g. pyrochlore) are the most common hosts for REE.
The relationships in Fig. 8e emphasise that S isotopes and ore mineralogy are governed by similar factors (i.e. the fluid composition and its redox evolution). These data suggest that d 34 S measurements could be used to make inferences about the ore mineralogy and we suggest that this relationship could benefit mineral exploration at poorly exposed or unexplored alkaline igneous systems. For example, by analysing a few samples containing S bearing minerals one could evaluate the d 34 S trajectory of the alkaline magmatic-hydrothermal fluid, place constraints on the redox evolution, and hence the REE-bearing ore minerals that are likely present. Sulfate-sulfide mineral pairs are known to occur at many alkaline complexes (Mitchell and Krouse, 1975;Mäkelä and Vartiainen, 1978;Gomide et al., 2013;Song et al., 2016;Hutchison et al., 2019) and while this approach would not provide unambiguous proof of mineralisation it could help guide exploration, and the strategy for expensive remote sensing surveys (since the spectral, magnetic, and radiometric response of these different REE-bearing ore minerals is variable, Neave et al., 2016;Zimmermann et al., 2016).
A final point is that all the alkaline igneous ore minerals highlighted above can be altered by low temperature hydrothermal processes (e.g. Mitchell and Liferovich, 2006;Salvi and Williams-Jones, 2006;Borst et al., 2016). Although HFSE are generally considered immobile (Lipin and McKay, 1989) there is increasing evidence that hydrothermal alteration causes variable remobilisation of HFSE in alkaline igneous systems (Salvi et al., 2000;Vasyukova et al., 2016;van de Ven et al., 2019). For example, peralkaline granitic systems such as Strange Lake show evidence for extensive hydrothermal remobilization of HFSE (Salvi and Williams-Jones, 1996), while peralkaline syenitic systems such as Thor Lake (Basal Zone deposit) and Ilímaussaq (Kringlerne deposit) show very minor remobilization of HFSE by hydrothermal fluids (despite significant alteration of primary ore minerals, Mö ller and Williams-Jones, 2017;Borst et al., 2016;van de Ven et al., 2019). Our d 34 S compilation demonstrates that the TPpR of alkaline igneous systems is more variable than any other ore-forming system and that there are significant variations in the TPpR trajectories of different alkaline complexes. These data imply large variations in pH and fO 2 during hydrothermal activity and may go some way to explaining why hydrothermal alteration is so variable between complexes and, importantly, why large variations in HFSE remobilization are observed at otherwise comparable deposit types.

Model caveats and future opportunities
The above examples illustrate how one can take independent constraints on TPpR of a magmatichydrothermal system, model S isotope evolution and carefully examine model-data fit to understand S sources, fluid evolution and disequilibrium processes. It is important to reiterate that all our models assume a fixed d 34 S P S as well as chemical and isotopic equilibrium in the evolving fluid. The models reproduce well the d 34 S evolution of porphyry, high sulfidation epithermal and alkaline igneous fluids but struggle for SMS and low sulfidation epithermal systems (where numerous S reservoirs with different d 34 S P S are involved and/or sulfate-sulfide equilibrium is not obtained).
Our forward models offer a new means of validating and testing hypotheses of magmatic fluid evolution and ore formation. When applying these models to other systems it is important to note the following caveats. Firstly, systems where d 34 S P S has changed significantly during fluid evolution should be avoided. The most obvious examples are those that have undergone a significant phase separation (i.e. loss of an immiscible fluid or gas). In such cases the impact on d 34 S P S will depend on the composition and redox state of both the immiscible phase and the hydrothermal fluid, as well as the fraction of total S that is separated. Phase separation is often to linked to increased oxidation of the hydrothermal fluid (by removal of reduced phases, e.g. H 2 , CO 2 and CH 4 , Drummond and Ohmoto, 1985;Palin and Xu, 2000). With respect to the S species, this would cause aqueous H 2 S to be progressively oxidized and would lead to partitioning of 34 S into SO 2 and SO 4 2and concentrate 32 S in the remaining H 2 S fluid phase. This scenario was modelled by McKibben and Eldridge (1990) who showed that if these oxidised species were subsequently removed then the residual S pool could be depleted by 3‰ at 10% S removal, 9‰ at 25% S removal and 21‰ at 50% S removal (at temperatures of 200°C, typical of epithermal systems). Although these values suggest that substantial depletions in d 34 S P S are possible they will only occur when S-bearing minerals are precipitated from a small finite S reservoir. Our study found no evidence of this reservoir effect, although we recognize that it could be important in small volume hydrothermal veins which are unconnected to or evolve separately from the gross magmatic-hydrothermal system. Systems that have undergone phase separation (e.g. boiling or rapid depressurisation) are easily identified by porous sulfide textures (Román et al., 2019) as well as cyclic trace element and d 34 S zonation in individual sulfide grains (McKibben and Eldridge, 1990;Peterson & Mavrogenes, 2014). Importantly, most natural samples that show evidence of phase separation and extreme d 34 S variation are associated with vein networks and, on this basis, we caution the use of our S isotope model for small volume hydrothermal veins.
As with all models it is also important to carefully consider the input parameters (TPpR) and whether these accurately capture the evolution of the magmatic-hydrothermal system. An understanding of fluid redox is critical when applying these models. In our study, we assumed that petrological estimates of fO 2 correspond to the fO 2 of the co-existing fluid (Kawasumi and Chiba, 2017) and used this value to evaluate S speciation. A key assumption when using the calculated fO 2 to define redox state is that O 2 must be considered a 'perfectly mobile component' (Korzhinskii, 1959). While this assumption likely holds for systems with a high fluid/rock ratio (similar to those studied here) it may breakdown at greater depths where fluid/rock ratios are low and most redox reactions take place between solid oxides and silicates (i.e. O 2 is no longer a mobile species because O is mainly bonded to inert components such as FeO and Fe 2 O 3 ). This is the case for subduction melanges formed at the slab-mantle interface, where oxygen molar quantity (nO 2 ) becomes the best descriptor of redox conditions (Tumiati et al., 2015). Another vital point regarding fO 2 is that while we model S isotope evolution at fixed fO 2 buffers (e.g. QFM) it is highly unlikely that any natural system would follow a single buffer throughout its evolution (Frost, 1991;Anenburg and O'Neill, 2019). Most magmatic-hydrothermal systems show significant redox heterogeneity, and this is particularly pronounced in hydrothermal veins (as fluids fluid transition between being rock buffered and fluid buffered, Peterson and Mavrogenes, 2014). S isotope evolution should not be considered along fixed redox buffers; many systems exhibit rapid shifts in redox over small temperature intervals and our models should be used to explore these (as shown by our example of Ilímaussaq-type fluids, Section 4.3, Fig. 8a, b).
A final point is that isotope modelling should always be accompanied by a detailed textural and mineralogical study to carefully evaluate sulfate-sulfide equilibrium and avoid spurious interpretations (e.g. associating supergene S phases, like jarosite, with magmatic-hydrothermal process, Fig. 4f). New analytical techniques and the advent of multiple S isotopes (Mason et al., 2006;Mandeville et al., 2009;Brueckner et al., 2015;Lode et al., 2015;Caruso et al., 2018;LaFlamme et al., 2016LaFlamme et al., , 2018aLaFlamme et al., , 2018b are particularly useful for understanding how of S-bearing minerals were precipitated and for fingerprinting S sources in settings where sulfate-sulfide disequilibrium is common (e.g. SMS, Ono et al., 2007;Peters et al., 2011;Aoyama et al., 2014, McDermott et al., 2015. The purpose of our modelling is to complement these approaches by allowing new d 34 S data to be understood in the context of fluid TPpR evolution. Future work could explore the S speciation outputs from the model (see Supplementary Information) because these make important predictions about the modal abundance of S-bearing minerals which could also be compared to natural observations. Further, since reduced S is often directly involved in complexing and mobilizing economically important (chalcophile) metals, future work should also attempt to delineate mobile and immobile fields for particular metals on the d 34 S sulfate -d 34 S sulfide diagrams. This would greatly improve our understanding of the links between S isotopes, metal transport and ore mineral precipitation.

CONCLUSIONS
Our study has updated and expanded the classic Ohmoto (1972) models for S isotope fractionation in magmatic-hydrothermal fluids. Our new models (Fig. 4) are in good agreement with his earlier study and confirm that the S isotope trajectory of an ore-forming hydrothermal fluid is controlled by the ratio of oxidised to reduced S, which is itself governed by physico-chemical parameters of the fluid (temperature, fO 2 and pH being of greatest importance). By fixing fluid fO 2 to a specific mineral redox buffer this allows us to directly calculate SO 4 2À /H 2 S and isotope fractionation for geologically realistic magmatichydrothermal fluids. We find that unless these fluids are highly reduced ((QFM) or oxidised (!MH), SO 4 2À /H 2 S and d 34 S will vary significantly during cooling (Fig. 4). We show that d 34 S P S can be reconstructed by combining our modelling with d 34 S sulfate -d 34 S sulfide diagrams but emphasize that previous attempts to calculate d 34 S P S (which assume a linear extrapolation at constant SO 4 2À /H 2 S, Fifarek and Rye, 2005) are unlikely to yield reliable results.
Comparing our equilibrium models to S isotope data from different ore systems is instructive and allows us to evaluate the effects of variable TPpR trajectories, S sources and S mineral disequilibrium. Porphyry fluids match the modelled TPpR trajectory well and require magmatic S sources of $5‰ (values which are equivalent to the d 34 S P S measured in calc-alkaline volcanic rocks, McKibben et al., 1996;Mandeville et al., 1998Mandeville et al., , 2009). Our models (e.g. NNO Fig. 5a-c) predict a temporal evolution of S isotopes from low d 34 S sulfide in zones of potassic alteration (associated with near-neutral mineralizing fluids) to high d 34 S sulfide in zones of phyllic and argillic alteration (associated with increasingly acidic fluids). Such trends may not be apparent at all porphyry deposits but, importantly, have been identified at several complexes (e.g. Porphyry Au-Cu bodies in SE Australia, Wilson et al., 2007). Epithermal deposits suggest similar S sources to porphyry fluids (0-5‰) and data from high sulfidation (acidic) systems match equilibrium model predictions (Fig. 6f). While our epithermal models predict marked differences in d 34 S sulfate -d 34 S sulfide between high and low sulfidation systems, the high pH (4-6) of the latter requires years to acquire sulfate-sulfide equilibrium and thus timescales of fluid advection and mineral precipitation become the main control on d 34 S in low sulfidation systems (Ohmoto and Lasaga, 1982;Rye, 2005). Active seafloor vents (SMS deposits) generally show little compliance with our equilibrium TPpR trajectory (Fig. 6d). They are a good example of how forward models can be used to flag disequilibrium and, in the case of SMS deposits, the very different S sources for sulfide and sulfate minerals. Finally, the origin of alkaline igneous ore deposits has been difficult to interpret because of their mineralogical complexity and variable hydrothermal overprint (e.g. Smith et al., 2015). Our models and data (Fig. 8) show that alkaline igneous fluids display extreme variability in their TPpR trajectory consistent with the wide range of measured physico-chemical parameters (Table 1). Compared to other ore-forming environments, alkaline systems show the greatest variability in TPpR and it is here that our model-data comparison provides new insights into why these systems are so variable in terms of their hydrothermal alteration and metal remobilization (Section 4.4).
New analytical techniques are transforming our ability to measure (multiple) S isotopes and are providing detailed textural information on the origins and evolution of S in active and extinct magmatic-hydrothermal systems (e.g. LaFlamme et al., 2016LaFlamme et al., , 2018aLaFlamme et al., , 2018bBolhar et al., 2020). The modelling framework developed herein is highly complementary and will allow these exciting new data to be carefully interrogated and, hopefully, understood.

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.

ACKNOWLEDGMENTS
This project was supported by the European Union's Horizon 2020 research and innovation programme under grant agreement No. 689909. W.H. also acknowledges support from a UKRI Future Leaders Fellowship (MR/ S033505/1). A.J.B. is funded by the NERC National Environment Isotope Facility award (NE/S011587/1) and the Scottish Universities Environmental Research Centre. Constructive and thorough reviews by C. LaFlamme, J.B. Walters, and an anonymous reviewer greatly improved this manuscript. A.C. Simon is also thanked for editorial handling. This manuscript stemmed from discussions with experts in porphyry and epithermal systems at a HiTech AlkCarb consortium meeting in Czech Republic. W.H. acknowledges K. Smith, F. Wall, J. Kynický who organised this event.

RESEARCH DATA FOR THIS ARTICLE
The isotopic data compiled for this study are provided in the Electronic Annex. These include separate spreadsheets for porphyry, epithermal, SMS and alkaline igneous systems. The MATLAB code written for this study is available from the corresponding author.