The behaviour of sulfur in magmas is complex because it dissolves as both sulfide (S2−) and sulfate (S6+) in silicate melt. Interesting aspects of the behaviour of sulfur are the solubility minimum (SSmin) and maxima (SSmax) observed with varying oxygen fugacity (fO2). We use a simple ternary model (silicate–S2–O2) to explore the varying fO2 paths where these phenomena occur. Both SSmin and SSmax occur when S2− and S6+ are present in the silicate melt in similar quantities owing to the differing solubility mechanisms of melt species containing these oxidation states of sulfur. At constant T, a minimum in dissolved total S content in vapour-saturated silicate melt (wSTm) occurs along paths of increasing fO2 and either constant fS2 or constant P. For paths on which wSTm is held constant with increasing fO2, the SSmin is expressed as a maximum in P. The SSmin occurs when the fraction of S6+ in the melt ([S6+/ST]m) is 0.25 for constant fS2 and [S6+/ST]m 0.75 for constant wSTm and P. A minimum in wSTm is not encountered during closed- or open-system depressurization in the simple system we modelled. However, the SSmin marks a change from reduction to oxidation during degassing. Various SSmax occur when the silicate melt is multiply saturated with at least two phases: vapour, sulfide melt, and/or anhydrite. The SSmin and SSmax are potentially important features of magmatic process involving S, such as mantle melting, magma mixing, and degassing. These concepts influence calculations of the pressures of vapour-saturation, fO2, and SO2 emissions using melt inclusions.

Supplementary material: Additional information and data used to create the figures are available at https://doi.org/10.6084/m9.figshare.c.6274527

Thematic collection: This article is part of the Sulfur in the Earth system collection available at: https://www.lyellcollection.org/topic/collections/sulfur-in-the-earth-system

It is widely accepted that there is a minimum in the solubility of sulfur in silicate melts (SSmin = sulfur solubility minimum; i.e. in the concentration of dissolved S in a silicate melt coexisting with an S-bearing vapour) as a function of oxygen fugacity (fO2). This has been proposed to occur when the speciation of S dissolved in the silicate melt changes from dominantly sulfide-bearing (S2−) to dominantly sulfate-bearing (SO42 or S6+) species, based on experimental studies of the concentration of S in vapour-saturated silicate melts spanning a range in fO2 (e.g. Fincham and Richardson 1954; Katsura and Nagashima 1974; Carroll and Rutherford 1985; Backnaes and Deubener 2011; Lesne et al. 2015; Matjuschkin et al. 2016; Nash et al. 2019) and on thermodynamic modelling (e.g. Moretti et al. 2003; Moretti and Ottonello 2005; Baker and Moretti 2011; Cicconi et al. 2020; Papale et al. 2022). Such a minimum in S solubility has implications for magmatic and volcanic processes. For example, any process where fO2 progressively changes and becomes closer to the fO2 of SSmin (e.g. mixing, progressive reduction or oxidation, degassing, etc.) will result in a decrease in the S solubility. This minimum has been used as evidence for the presence of additional, low-solubility, potentially unquenchable, S-bearing species in silicate melts that could be important for metal transport in arc systems (Matjuschkin et al. 2016). Also, understanding the thermodynamic basis for this feature is critical for calculating the pressure of vapour-saturation of S-bearing magmas using the volatile concentrations of quenched glasses (e.g. Lesne et al. 2015) and the composition of vapour and silicate melt during degassing (e.g. Moretti et al. 2003; Gaillard and Scaillet 2014; Liggins et al. 2020; Aiuppa et al. 2022).

Despite its potential importance, the existence of an SSmin with varying fO2 has been somewhat mischaracterized in the literature. As emphasized by O'Neill (2021), this is at least in part due to a lack of clarity regarding the independent variables and the path followed by sulfur fugacity (fS2) with increasing fO2 for specific natural or experimental processes (e.g. Moretti et al. 2003) and the number and identity of additional S-bearing phase(s) with which the silicate melt is saturated (e.g. Jugo et al. 2005). Thermodynamic modelling has shown that there is a minimum in the calculated S contents of melts with increasing fO2 if fSO2 (or the fugacity of another S-bearing species) is kept constant (e.g. Moretti et al. 2003; Moretti and Ottonello 2005; Baker and Moretti 2011; Cicconi et al. 2020; Papale et al. 2022). In this paper, we expand upon these previous studies to explore the conditions and paths for which an SSmin occurs and some of the implications for magmatic and volcanic processes. We also expand upon the work of Jugo (2009) regarding a S solubility maximum (referred to as an SSmax) for silicate melts that are multiply saturated with sulfide melt + anhydrite ± vapour. Although silicate melt and vapour in most natural systems contain H, C, halogens, metals, etc. in addition to S and O, in this paper we limit ourselves to a simple system in which S and O are the only volatile components in the silicate melt (such a system may be appropriate for Jupiter's moon Io; e.g. Zolotov and Fegley 2000). By limiting our treatment to this simple end-member system, the factors leading to an SSmin or an SSmax can be more easily isolated and understood.

Equilibria between silicate melt, vapour, sulfide melt and anhydrite

Sulfur occurs in several phases in magmatic systems, including dissolved S-bearing species in silicate melt, gaseous species in vapour, immiscible sulfide melts and various sulfate phases (e.g. reviews by Parat et al. 2011; Wallace and Edmonds 2011). In addition to silicate melt (silm) and vapour (v), we consider pure Fe-sulfide melt (FeS; sulfm) and anhydrite (CaSO4; anh) (Fig. 1a and b). For our treatment, the silicate melt end-member can be compositionally simple (e.g. SiO2, CaMgSi2O6, NaAlSi3O8) or complex (e.g. a natural basalt), provided it is fixed in composition. We assume that the silicate melt is insoluble in the vapour.

The vapour is assumed to be constrained to the S2–O2 binary subsystem and to contain only three species (S2, O2 and SO2; Fig. 1b), hence
xS2v+xO2v+xSO2v=1
(1)
where xiv is the mole fraction of species i in the vapour (v). Other species are present in an S2–O2 vapour (e.g. SO3, SO, S polymers, etc.) and could be added to our treatment. However, the three species in equation (1) are the most significant (e.g. Oppenheimer et al. 2011; Renggli et al. 2017; Henley and Seward 2018; Henley and Fischer 2021) and are sufficient to illustrate the salient points.
The composition of the silicate melt is confined to a ternary system with the following three components: silicate, S2 and O2. Sulfide (S2−) and sulfate (S6+) are assumed to be the only significant S-bearing species dissolved in natural silicate melts based on XANES measurements and solubility experiments (e.g. Fincham and Richardson 1954; Paris et al. 2001; Métrich et al. 2009; Wilke et al. 2011). Intermediate S-bearing species have been observed or inferred (e.g. S4+, S0, molecular SO2: Clemente et al. 2004; Métrich and Wallace 2008; Burgisser et al. 2015; Lesne et al. 2015; Matjuschkin et al. 2016), but these species are not thought to be significant in natural (especially in Fe-bearing) silicate melts (e.g. Jugo et al. 2010; Klimm et al. 2012). Hence, the silicate melt is assumed to contain only two S-bearing species (S2− and S6+; Fig. 1a), such that
wS2m+wS6+m=wSTm
(2)
where wim is the weight fraction of species i in the silicate melt (m), and ST refers to the total dissolved S content.
Three independent reactions control the coexisting compositions of silicate melt and vapour in this system. The first reaction describes a homogeneous equilibrium that governs the speciation of the vapour:
0.5S2(v)+O2(v)=SO2(v)
(3rma)
which is governed at equilibrium by
K3(T)=fSO2(fS2)0.5fO2
(3rmb)
where K3 is the equilibrium constant for reaction (3a), which depends on temperature (T), and fi is the fugacity of species i in the vapour. If the treatment were to include other S ± O-bearing vapour species (such as SO3, SO, etc.; see above), a statement of homogeneous equilibrium would have to be added for each additional vapour species.
The second reaction describes a heterogeneous equilibrium between silicate melt and vapour that governs the dissolution of sulfur from the vapour as S2− in the silicate melt (e.g. Fincham and Richardson 1954; Moretti and Ottonello 2003, 2005; Moretti and Papale 2004; Gaillard and Scaillet 2009, 2014; Baker and Moretti 2011; Gaillard et al. 2011, 2013, 2015; Baumgartner et al. 2017; Moretti 2021), as given by the following three relations:
0.5S2(v)+O2(m)=0.5O2(v)+S2(m)
(4rma)
K4(P,T)=aS2maO2m(fO2fS2)0.5xS2mxO2m(fO2fS2)0.5
(4rmb)
and
CS2=wS2m(fO2fS2)0.5
(4rmc)
where aim, xim and wim are the activity, mole fraction and weight fraction, respectively, in the silicate melt of the species i (either S2− or oxide [O2−], in this case); and CS2 is defined as the ‘sulfide capacity’ (e.g. Fincham and Richardson 1954; O'Neill 2021). The final term of equation (4b) makes the approximation that aim can be replaced by xim, and this approximation is adopted throughout. Given this approximation, CS2 is simply related to K4, the equilibrium constant for reaction (4a), although by convention it is defined in terms of wim rather than xim. Alternatively, the activity coefficients can be incorporated into CS2 (O'Neill 2021). Finally, we assume throughout that xO2m (i.e. the O2− in the silicate melt that can react with S2 to produce S2− via reaction (4a)) can be approximated as constant (i.e. xS2mxO2m).
A third reaction describes an additional heterogeneous equilibrium between silicate melt and vapour that governs the dissolution of sulfur from the vapour as SO42 in the silicate melt (e.g. Fincham and Richardson 1954; Moretti and Ottonello 2003, 2005; Moretti and Papale 2004; Baker and Moretti 2011; Moretti 2021):
0.5S2(v)+1.5O2(v)+O2(m)=SO42(m)
(5rma)
K5(P,T)=aS6+maO2m(fS2)0.5(fO2)1.5xS6+mxO2m(fS2)0.5(fO2)1.5
(5rmb)
and
CS6+=wS6+m(fS2)0.5(fO2)1.5
(5rmc)
where aS6+m and xS6+m are the activity and mole fraction, respectively, of sulfate dissolved in the silicate melt, and CS6+ is referred to as the ‘sulfate capacity’ (Fincham and Richardson 1954; Boulliung and Wood 2022; O'Neill and Mavrogenes 2022). Again, CS6+ is simply related to K5, the equilibrium constant for reaction (5a), but using weight fraction instead of mole fraction and assuming xO2m is constant.
An alternative heterogeneous equilibrium between silicate melt and vapour could be used instead of either reaction (4a) or (5a) to describe the conversion of S2− to S6+ in the silicate melt (e.g. Wallace and Carmichael 1994; Matthews et al. 1999; Moretti and Ottonello 2005; Métrich et al. 2009; Jugo et al. 2010; Baumgartner et al. 2017; Moretti 2021):
S2(m)+2O2(v)=SO42(m)
(6rma)
K6(P,T)=aS6+maS2m(fO2)2xS6+mxS2m(fO2)2
(6rmb)
and
xS6+mxS2m=[S6+S2]m=wS6+mwS2m=CS6+CS2(fO2)2.
(6rmc)
Reaction (6a) can be obtained by subtracting reaction (4a) from reaction (5a) and rearranging. Reaction (6a) is useful as it emphasizes that the oxidation state of S dissolved in the silicate melt (i.e. [S6+/S2−]m or [S6+/ST]m) at a given T and P is controlled only by fO2 provided that CS2 and CS6+ are known. However, CS2 and CS6+ are highly dependent on the composition of the silicate melt (e.g. O'Neill and Mavrogenes 2002, 2022; Moretti and Ottonello 2005; Nash et al. 2019; O'Neill 2021; Boulliung and Wood 2022; Moretti 2021). Additionally, these capacities will depend on T and P (because K4 and K5 must depend on T and P unless the standard state enthalpy and volume changes of the reactions are zero); moreover, they are likely to depend on the speciation of other multivalent elements in the silicate melt (e.g. [Fe3+/Fe2+]m; O'Neill 2021). Such factors could thus lead indirectly to variations in [S6+/S2−]m at constant fO2. However, these are expected to be minor effects for most for the examples considered here owing to the small amount of Fe3+ in most magmas (i.e. [Fe3+/FeT]m < 0.3; e.g. Gaborieau et al. 2020) and the relatively low P considered here. Therefore, we assume CS2 and CS6+ depend only on T and on the composition of the silicate component (and only on total Fe, not [Fe3+/Fe2+]m).
Given equation (1), if an S–O vapour is present in the system (i.e. if the silicate melt is vapour-saturated), the sum of the partial pressures (pi) of the species in the vapour equals the total pressure of the system (P):
P=pO2+pS2+pSO2
(7)
where the partial pressures are related to fugacity and mole fraction in the vapour through fugacity coefficients (γi):
fi=γipi=γixivP.
(8)
When the silicate melt is saturated with sulfide melt, the chemical potential of FeS in the silicate melt (μFeSm) and the sulfide melt are equal. This equilibrium condition has been parameterized by the ‘sulfide content at sulfide-saturation’ (S2−CSS), the dissolved S2− concentration in the silicate melt in equilibrium with sulfide melt (e.g. Shima and Naldrett 1975; O'Neill and Mavrogenes 2002; Smythe et al. 2017; O'Neill 2021):
wS2m=wS2CSSm.
(9)
The total S content of a silicate melt that is saturated with sulfide melt (STCSS) is then given by equation (9) in combination with equations (2) and (6c):
wSTm=wSTCSSm=[1+(CS6+CS2)(fO2)2]wS2CSSm.
(10)
Likewise, when the silicate melt is saturated with anhydrite, the chemical potential of CaSO4 in the silicate melt (μCaSO4m) and anhydrite are equal. This equilibrium condition has been parameterized by the ‘sulfate content at anhydrite-saturation’ (S6+CAS), the dissolved S6+ concentration in the silicate melt in equilibrium with anhydrite (e.g. Baker and Moretti 2011; Chowdhury and Dasgupta 2019; Zajacz and Tsay 2019):
wS6+m=wS6+CASm.
(11)
The total S content of a silicate melt that is saturated with anhydrite (STCAS) is then given by equation (11) in combination with equations (2) and (6c):
wSTm=wSTCASm=[1+(CS2CS6+)(fO2)2]wS6+CASm.
(12)
Iron is the only other multivalent element considered in our system, which can be dissolved in the melt as Fe2+ (FeO) or Fe3+ (FeO1.5), where [Fe3+/Fe2+]m is related to fO2 (e.g. Kress and Carmichael 1991).

The importance of the phase rule in our treatment of S-solubility in silicate melt

In its simplest form, the phase rule relates the number of components (c) and phases (φ) in a system to the variance (or the degrees of freedom, F) of the assemblage:
F=c+2φ.
(13)
For most of the calculations presented here, our system has three components (silicate, S2, and O2; c = 3) and two phases (silm + v; φ = 2) (Fig. 1a and b). This silm + v assemblage is thus trivariant (φ = 2, F = 3). Therefore, for the silm + v assemblage, if any three linearly independent intensive variables are chosen as independent variables, the state of each phase in the system is fully defined, and all other intensive variables in the coexisting phases are dependent variables (e.g. Prigogine and Defay 1954). Potential independent variables include P, T and μi (or other intensive variables in the coexisting phases such as fugacities or partial pressures (fi, pi), phase compositions (xiv, xim), etc.). Consequently, the values of three independent intensive variables must be given to specify completely the states of coexisting silicate melt and vapour in the system. The values of all other intensive variables in these two coexisting phases can then be calculated given the values chosen for the independent variables and knowledge of the thermochemical properties of the silicate melt and vapour phases. The phase rule alone does not allow the masses or mass fractions of the vapour and silicate melt in the system to be determined. We also present calculations of closed- and open-system degassing as a function of decreasing P at constant T. The thermodynamic basis for these calculations differs slightly from those described here and are explained in the section ‘Isothermal, decompression-induced degassing’.

The important point here is a restatement of the cautionary note of O'Neill (2021) and previous studies (e.g. Moretti et al. 2003; Moretti and Ottonello 2005; Baker and Moretti 2011; Cicconi et al. 2020; Papale et al. 2022) about the SSmin when fO2 is an independent variable; that is, assuming T is constant, the variation of wSTm is not uniquely defined if only fO2 is independently varied for a silm + v assemblage. The behaviour of a third variable must also be specified for the state of the silm + v assemblage to be defined at each point on a path of varying fO2. Only then can the variation in wSTm (including the nature of any minimum or maximum) as a function of fO2 be uniquely characterized.

As emphasized above, there are a variety of intensive variables that could be chosen as independent variables in our model system. In the calculations presented below, we use combinations of the following variables as independent variables: T, P, fO2, fS2, wSTm, μFeS, μCaSO4 and bulk composition. There are other possible choices of independent intensive variables that we have not considered but that could be of interest. The chosen variables must be linearly independent and thus some combinations of variables are not allowed. For example, because the vapour is assumed to consist only of S and O, xSTv+xOTv=1, so only one of these two variables can be chosen independently. Details on how the values of dependent intensive variables of interest in the silm + v system are solved for given the various choices of three independent variables are given in the Supplementary material.

Details of our choices of thermodynamic parameters for vapour and silicate melt (in our case a Hawaiian tholeiite) can be found in the Supplementary material. The compositions of sulfide melt and anhydrite do not plot in the plane of our chosen three-component system (Fig. 1b). Despite this, as long as the silicate melt is confined to the model ternary system, we can still model the effects of sulfide melt- and/or anhydrite-saturation (i.e. fixed μFeS and/or μCaSO4) on the properties of the silicate melt and vapour phases using parameterizations of S2−CSS and/or S6+CAS as described in the section ‘Variations in wSTm along isothermal and isobaric paths of increasing fO2’.

We use this conceptual background to explore trends in and interrelationships among various choices of independent and dependent variables, focusing on the implications for the SSmin and the SSmax. For all of our calculations, T and the base composition of the silicate component are held constant. Because we assume that CS2 and CS6+ depend only on T and the composition of the silicate component, both these parameters are the same in all calculations presented here. In particular, CS2 and CS6+ are assumed independent of P and [Fe3+/FeT]m and therefore of fO2. It is important to emphasize that the results based on our chosen parameters only describe the representative behaviour of the particular Hawaiian tholeiite melt composition used in our modelling. Although we are confident that the trends and insights derived from this choice are robust, the exact behaviour of other melt compositions will depend on the chosen values of CS2 and CS6+. Consequently, the specific values of various variables, including the precise values of fO2 where shifts in behaviour are predicted to occur, will probably vary with T, the composition of the silicate component, and the CS2 and CS6+ parameterization used (e.g. O'Neill and Mavrogenes 2002, 2022; Moretti and Ottonello 2005; Nash et al. 2019; O'Neill 2021; Boulliung and Wood 2022; Moretti 2021).

In this section we calculate the state of the system by choosing T (1200°C) and a value of fO2 (−5 < ΔFMQ < +5, where FMQ is the Fayalite–Magnetite–Quartz buffer; FMβQ in the study by Frost 1991) as two of the independent variables. The other independent variables considered are fS2, wSTm, P, μFeS or μCaSO4. For calculations in this section, the bulk compositions of the coexisting phases are not independent variables and hence are not constant. As a result [S6+/ST]m and [Fe3+/FeT]m vary with fO2 and therefore wO2m also varies (wO2m is the amount of O2 in the silicate melt in excess of the amount in the silicate component), as does the S/O ratio of the vapour. For any given state of the system (i.e. where the values of all intensive variables are defined in the coexisting phases), the results will be the same regardless of which three variables are chosen as independent. Therefore, the following figures convey the same results, but with different variables plotted on the axes and contours. For most of the P range shown in the following figures, γi1, hence fipi using equation (8) (e.g. for log10[P, bar] < 3, γi<1.25 for all species). The section ‘Silicate melt + vapour’ describes results for the two-phase silm + v assemblages, and the sections ‘Multiply saturated silicate melt with vapour, sulfide melt and/or anhydrite' and ‘Variations in wSTm along isothermal and isobaric paths of increasing fO2’ describe assemblages of silm + two other phases (v and/or sulfm and/or anh).

Silicate melt + vapour

For the silm + v assemblages, φ = 2 and F = 3; therefore, in addition to T and fO2, only one other independent variable is needed to specify fully the state of the system. This third independent variable will be referred to as the ‘Y’ variable: we first choose Y=fS2, then Y=wSTm, and finally Y = P. For each of these choices of the independent variables, we solve the system of equations (3b), (7), (8) and two out of (4c), (5c) and (6c) (details are given in the Supplementary material).

S speciation in the silicate melt and vapour

Figure 2 shows how the compositions of the silicate melt and vapour change with varying fO2 at constant T and either constant fS2 (Fig. 2a and b), constant wSTm (Fig. 2c and d) or constant P (Fig. 2e and f) given that the silicate melt is vapour-saturated. As fO2 changes, the speciation of sulfur in the silicate melt ([S6+/ST]m) and in the vapour (xSO2v) changes. This leads to changes with fO2 in the speciation of S in the silicate melt and vapour, which is indicated by the coloured horizontal bars in Figure 2. For the vapour, S2 is the dominant species at low fO2 (xSO2v<0.1, the fO2 range of the horizontal purple bar in each column), SO2 is the dominant species at high fO2 (xSO2v>0.9, turquoise–green–yellow) and between these there is a transition from dominantly S2 to SO2 (0.1xSO2v0.9, blue). At sufficiently high fO2, O2 becomes more abundant than S2 in the vapour (i.e. to the right of the points labelled α in the yellow bars in Fig. 2b, d and f), but both are much less abundant than SO2 in the range shown. For the silicate melt, S2− is the dominant S-bearing species at low fO2 ([S6+/ST]m < 0.1, purple–blue–turquoise), S6+ is dominant at high fO2 ([S6+/ST]m > 0.9, yellow) and between these there is a transition from dominantly S2− to S6+ (0.1 ≤ [S6+/ST]m ≤ 0.9, green).

At constant T, the effects of varying the third independent variable (Y=fS2, wSTm or P) in addition to fO2 on the dependent variables (which are referred to as ‘Z’ variables) are shown using contour plots in Figures 35. The colour scheme from Figure 2 showing the dominant silicate melt and vapour species and their changes is also shown in part (d) of Figures 35. Comparison of these figure parts shows that the topology of the silicate melt and vapour speciation (indicated by the coloured fields) is similar for all three of these choices of Y. The regions in Figures 35 are separated by black curves, where dashed curves are isopleths of vapour speciation (xSO2v=0.1 and 0.9) and continuous curves are isopleths of silicate melt speciation ([S6+/ST]m = 0.1 and 0.9).

A key point of part (d) in Figures 35 is that regardless of the choice of the Y variable, at low fO2 (ΔFMQ ≲ +0.7) nearly all dissolved S is present as S2− (purple–blue–turquoise). Under these conditions wSTmwS2m, and wSTm is controlled by equation (4c), giving wS2m as a simple function of fO2 and fS2. Likewise, at sufficiently high fO2 (ΔFMQ ≳ +1.7), nearly all S is dissolved as S6+ (yellow). Here wSTmwS6+m, and wSTm is instead controlled by equation (5c), which gives wS6+m as another simple function of fO2 and fS2. At intermediate fO2 (+0.7 ≲ ΔFMQ ≲ +1.7), the silicate melt changes from being dominated by S2− to being dominated by S6+ (green). The green region is relatively narrow in these figures because the S speciation in the silicate melt is defined by equation (6c), which depends on (fO2)2 (e.g. the difference in log10[fO2] between [S6+/ST]m = 0.1 and 0.9 is log10[9] ∼ 0.95). The silicate melt isopleths are not precisely vertical owing to the P dependence of fO2 on the FMQ buffer (Frost 1991). When the value of Y is large, P is large, which leads to pi deviating from fi, causing the silicate melt isopleths to deviate from vertical (i.e. to curve to lower fO2). They would, however, be vertical if plotted against log10[fO2] without normalization to a buffer, given that CS2 and CS6+ are assumed to be P-independent.

Unlike the speciation of S dissolved in the silicate melt, which depends only on fO2 via equation (6c), vapour speciation (i.e. xSO2v) depends on fS2 in addition to fO2 from equations (3b) and (8). Therefore, the boundaries separating where the vapour phase is dominated by S2 (purple), both S2 + SO2 (blue) and SO2 (turquoise–green–yellow) depend on fO2 (see the xSO2v isopleths in Figs 35). However, the slopes of these boundaries depend on the choice of the Y variable and its value, as well as the silicate melt speciation (e.g. see the change in slope of xSO2v=0.9 in Fig. 4 where it enters the green region with increasing fO2). The change in vapour speciation from dominantly S2 to SO2 occurs over a wider range of fO2 than the silicate melt speciation (compare the widths of the blue and green bands in Fig. 2 and in part (d) of Figs 35) because fSO2 (and xSO2v via equation (8)) depends on (fO2)1 from equation (3b), in contrast to [S6+/S2−]m, which depends on (fO2)2 from equation (6c).

Combining changes in silicate melt and vapour speciation, there are three main regions plus three transitional regions in log10[fO2]log10[Y] space (Figs 35): at low fO2 and high Y, the vapour is S2-dominated and the silicate melt is S2−-dominated (purple); at higher fO2 and lower Y, the vapour contains both S2 and SO2 in similar concentrations, whereas the silicate melt is still S2−-dominated (blue); at higher fO2 and lower Y, the vapour is now SO2-dominated and the silicate melt is still S2−-dominated (turquoise); at higher fO2 and most Y values shown, the vapour remains SO2-dominated but the silicate melt contains both S2− and S6+ in similar concentrations (green); at high fO2 and all Y values shown, the vapour is still SO2-dominated but the silicate melt is S6+-dominated (yellow); and there is also a small region where the vapour contains both S2 and SO2, and the silicate melt both S2− and S6+, in similar concentrations (orange region labelled [S2− + S6]m + [S2 + SO2]v in Figs 35: it is not intersected in Fig. 2). At higher Y values than shown in Figures 35, there are three additional regions: at intermediate fO2, the vapour is S2-dominated and the silicate melt contains both S2− and S6+ in similar concentrations; at high fO2 but intermediate Y, the vapour contains S2 and SO2 in similar concentrations and the silicate melt is S6+-dominated; and at higher Y, the vapour is S2-dominated and the silicate melt is still S6+-dominated (a sketch of this topology is shown in Supplementary material, Fig. S3).

Regions with a single dominant species in the silicate melt and vapour

For all three of our choices of a constant independent Y variable, curves (Fig. 2) and contours (Figs 35) of the dependent variables (Z) shown in the figures have constant slopes when there is a single dominant species in the silicate melt and a single dominant species in the vapour. Hence, when T is constant and there is only one dominant species in the silicate melt and one dominant species in the vapour,
Z(fO2)a(Y)b
(14)
for all choices of Z and Y given here. Therefore, contours of constant Z in log10(Y)log10(fO2) plots have slopes (σ) ∼ –a/b (Figs 35); and, when Y is constant, curves of different Z variables have slopes (ς) ∼ a (Fig. 2). Hence, the values of the slopes of the curves and contours of Z depend on the specific choice of the Y variable (see Table 1). It should be noted that ς and σ are approximately, rather than exactly, equal to these values because, (1) the x-axis is log10(fO2) relative to FMQ rather than strictly log10(fO2), (2) partial pressure rather than fugacity is sometimes plotted (e.g. Fig. 2), and these are related through equation (8), and (3) although one species is dominant, the concentrations of the other species are not zero.
At constant T, the relationships between fO2, Y (fS2, wSTm and P) and Z (fS2, fSO2, P, wS2m, wS6+m and wSTm) when the silicate melt is S2−-dominated and the vapour is S2-dominated (purple), the silicate melt is S2−-dominated but the vapour is SO2-dominated (turquoise), and the silicate melt is S6+-dominated and the vapour is SO2-dominated (yellow) can be derived by substituting
wSTmwS2m
(15rma)
when the silicate melt is S2−-dominated (purple–blue–turquoise), but
wSTmwS6+m
(15rmb)
when S6+-dominated (yellow). Similarly,
PpS2fS2
(16rma)
applies when the vapour is S2-dominated (purple), but
PpSO2fSO2
(16rmb)
when SO2-dominated (turquoise–green–yellow). Equations (15) and (16) are then substituted into equations (3)–(8) and rearranged into the form of equation (14). These relationships are derived in the Supplementary material, summarized (including ς and σ values) in Table 1 and labelled in Figures 25.

Regions with mixed speciation in the silicate melt or vapour

Changes in the slopes of curves (ς, Fig. 2) and contours (σ, Figs 35) occur for dependent variables (Z) when their slopes in log10(Y)log10(fO2) space differ in regions dominated by different species (i.e. S2 (purple) v. SO2 (turquoise) in the vapour and S2− (turquoise) v. S6+ (yellow) in the silicate melt). The regions in which the slopes change from one essentially constant value to another are where there is mixed speciation in the vapour (i.e. both S2 and SO2 in similar concentrations: blue) or in the silicate melt (i.e. both S2− and S6+ in similar concentrations: green).

If the log–log slope of a Z variable with respect to fO2 has the same sign in the low- and high-fO2 regions on either side of one of the mixed species regions, the slope simply becomes steeper or shallower along a constant Y path (Figs 25). For example, with increasing fO2 and Y = constant wSTm=400ppm, the slope in Figure 2d of log10(P) v. log10(fO2) changes at approximately the point labelled β. The slope is ς +1.0 where the silicate melt is S2−-dominated and the vapour is S2-dominated (in the purple region), changing to ς +1.5 where the vapour is SO2-dominated (in the turquoise region), but both slopes are positive (Table 1). For the contour plot in Figure 4c, this is shown as σ decreasing from −0.5 to −1.5, but both being negative. This is because when the silicate melt is S2−-dominated, the log–log slope for Z=pSO2 is steeper than for Z=pS2 (ς +1.5 v. +1.0 or σ −1.5 v. −0.5 in the purple–blue–turquoise regions in Figs 2d and 4a–c). Therefore, the curves cross when pS2=pSO2 (at the point labelled β in the blue region of Fig. 2d). Hence, along this path P is essentially equal to pS2 when S2 dominates the vapour (ς +1.0 or σ −0.5, purple), steepens continuously in the region where the concentration of SO2 in the vapour increases (blue), and is then well approximated by pSO2 when SO2 dominates the vapour (ς +1.5 or σ −1.5, turquoise).

If the log–log slope of a Z variable with respect to fO2 has opposite signs on the two sides of a mixed speciation region, a maximum (positive to negative slope with increasing fO2) or minimum (negative to positive slope with increasing fO2) in Z occurs (Figs 25). For example, for Y = constant fS2, there is a minimum in Z=wSTm at [S6+/ST]m = 0.25 where the silicate melt changes from being S2−- to S6+-dominated (γ in Fig. 2a; derivation of the [S6+/ST]m value at the wSTm minimum is given in the Supplementary material). This corresponds to maxima in Z=log10[wSTm] contours (Fig. 3c). This is an example of an SSmin on a specific path of increasing fO2 (i.e. fS2 and T are constant). A similar feature has been described for Y = constant fSO2 (e.g. Moretti et al. 2003; Moretti and Ottonello 2005; Baker and Moretti 2011; Cicconi et al. 2020; Papale et al. 2022), but the minimum will occur at [S6+/ST]m = 0.75 (derivation is given in the Supplementary material). However, there are paths with monotonically increasing fO2 and variable fS2 for which a minimum is not encountered (e.g. all linear paths in Fig. 3c with σ ≥ +1 or σ ≤ −3). There is also an SSmin when Y = constant P (δ in Fig. 2e) and the corresponding maxima in the log10[wSTm] contours (Fig. 5c). This minimum occurs at [S6+/ST]m 0.75, similar to constant fSO2 because SO2 is the main species in the vapour (i.e. in the green region; derivation is given in the Supplementary material). Similarly, there are paths with monotonically increasing fO2 and variable P that do not encounter the minimum (e.g. all log10(fO2)–log10(P) paths in Figure 5c with σ ≥ +1.5 and σ ≤ −0.5).

These SSmin occur because curves and contours of wSTm have opposite slopes at lower fO2 where S2− dominates the silicate melt compared with higher fO2 where S6+ dominates (Figs 2a, e, 3c and 5c). These slope changes occur because O2 in the vapour is on the product side of reaction (4a) when sulfur dissolves dominantly as S2− in the silicate melt, but on the reactant side when sulfur dissolves dominantly as S6+ in the silicate melt in reaction (5a). This is true whether the reactions describing silm + v heterogeneous equilibria are written with S2 as the vapour species (e.g. equations 4a and 5a) or with SO2 as the vapour species (not shown). Thus, for any path of constant fS2 or P (i.e. horizontal slices in Fig. 3c or Fig. 5c, respectively), wSTm reaches a minimum at [S6+/ST]m 0.25 or 0.75, respectively.

When Y=wSTm and is held constant as fO2 increases, the SSmin manifests itself as a maximum in Z = P (ε in Fig. 2d) and corresponding minima in Z = P contours (green region in Fig. 4c). This is because at the SSmin, a higher P is required to maintain the same S content in the vapour-saturated silicate melt. There are also maxima in Z=pSO2 or pS2 (ε and ζ in Fig. 2d, respectively) and corresponding minima in Z=fSO2 or fS2 contours (Fig. 4b and a, respectively). The change in sign of the slope of pSO2 and pS2 with increasing fO2 (e.g. from ς ∼ +1.5 to −0.5 for pSO2 in Fig. 2d) is due to the crossover of S2− to S6+ as the dominant species in the silicate melt (see the section ‘Regions with a single dominant species in the silicate melt and vapour’). Therefore, the maxima (and corresponding minima in the equivalent contours) in pSO2, pS2 and P all occur where the silicate melt speciation is mixed. Because SO2 is the dominant vapour species on both sides of the maximum in P, pS2 has negligible influence on the maximum in P (i.e. pS2pSO2 and therefore PpSO2). Hence, the maxima in pSO2 and P essentially coincide (ε) at [S6+/ST]m 0.75, but the pS2 maximum (ζ) is at [S6+/ST]m 0.25 (Figs 2d and 4a–c; derivation is given in the Supplementary material). However, the maximum in pSO2 (and P) is not exactly at [S6+/ST]m = 0.75 (or 0.25 for pS2) because (1) γSO2 (and γS2) are not one (causing a slight offset between fi and pi) and (2) because of the contribution of pS2 (and pO2) to P. These factors cause a very small offset between fSO2, pSO2 and P; however, [S6+/ST]m = 0.75 is a very good approximation for the position of the maximum.

We refer to the fO2 where the wSTm minimum or P maximum occurs as the SSminfO2; it can be calculated from [S6+/ST]m (i.e. 0.25 or 0.75) using equation (6c) and is dependent on CS2 and CS6+ (and therefore on melt composition and T). However, it can only be defined for silm + v at a given T and for a choice of Y=fS2, wSTm or P. The SSmin always occurs when both S2− and S6+ are dissolved in the silicate melt in similar concentrations, although it occurs at [S6+/ST]m = 0.25 when Y=fS2, but at [S6+/ST]m 0.75 when Y=wSTm or P (and SO2 dominates the vapour). Moreover, the minimum is not symmetric with changing fO2 (see Figs 2a, d, e and 35c) because the stoichiometries of the reactions for sulfur dissolution as S2− and S6+ differ. The different manifestations of the SSmin are related but not identical, reflecting the different fS2fO2 paths for different choices of Y. For example, there is a P maximum when Y=wSTm (Figs 2d and 4c), but no P maximum when Y=fS2 (Figs 2b and 3b) or P (Fig. 2f, by definition). Similarly, there is a wSTm minimum when Y=fS2 (Figs 2a and 3c) and P (Figs 2e and 5c), but not when Y=wSTm (Fig. 2c, again by definition). This highlights the importance of the choice of the third independent variable for silm + v in understanding the occurrence of the SSmin, which is required to specify the path followed by the system with changing fO2.

Multiply saturated silicate melt with vapour, sulfide melt and/or anhydrite

The trivariant two-phase silm + v assemblage is stable at Y values below the curves indicating sulfide melt-saturation (the grey continuous curves labelled ‘sulfm’) and anhydrite-saturation (the grey dashed curves labelled ‘anh’) in Figures 35. Above these grey curves, the silm + v assemblages are metastable, but at equilibrium the silicate melt is saturated only with sulfide melt and/or anhydrite, which will be discussed in the section ‘Variations in wSTm along isothermal and isobaric paths of increasing fO2’.

In this section, we first describe the behaviour of the system on the grey curves in Figures 35, where silm + v in the model ternary system is stable and saturated with sulfide melt and/or anhydrite. In our model ternary system, there are still only two phases present (silm + v) and therefore F = 3. Therefore, we assume that the silicate melt and vapour can still be described compositionally and thermodynamically by the ternary plane even though they are saturated with sulfide melt or anhydrite. Hence, given T and fO2, the curve for sulfide melt-saturation of the silicate melt is defined by a particular value of Y=μFeSm that is equal to that of FeS in the saturating sulfide melt. In practice, we apply the formulation of O'Neill (2021) to calculate S2−CSS for the silicate melt, which is conceptually equivalent to holding μFeSm equal to a constant chemical potential of FeS in sulfide melt at the relevant conditions. In our calculations, we assume the sulfide melt is pure FeS (i.e. μFeS=μFeS), although this could be modified. Similarly, the curve for anhydrite-saturation is defined by Y=μCaSO4m. Again, although S6+CAS is calculated based on Chowdhury and Dasgupta (2019), this is conceptually equivalent to holding μCaSO4m equal to a constant chemical potential of pure anhydrite at the relevant conditions. The stable saturation condition (i.e. whether the silicate melt is sulfide melt- or anhydrite-saturated) is the one with the lowest wSTm (see Fig. S4 in the Supplementary material).

If the silm + v assemblage is saturated with either sulfide melt or anhydrite, P and wSTm can no longer be independent variables. Therefore, having chosen values for three variables (i.e. T, fO2 and μFeSm or μCaSO4m, which confines the state of the system to the sulfm or anh curves), the state of the system is fully defined. All other intensive variables in the coexisting silicate melt and vapour phases are fixed (including both P and wSTm), as can be visualized in Figures 35. Moreover, when the silm + v assemblage is saturated with both sulfide melt and anhydrite at a fixed T, three independent variables are specified (T, μFeSm and μCaSO4m), so no other parameters can be chosen independently. This is demonstrated graphically in Figures 3c and 5c by the grey star, which for a given T has fixed values of P, fO2 and wSTm at the unique intersection of the grey continuous and dashed curves. We can thus describe quantitatively with our model a constant-T, increasing-fO2 path for an assemblage containing silm + v + (sulfm and/or anh) by following the grey curves. The variations in the values of the Z variables along this path can be read from the contours crossed by the grey curves in Figures 35, whereas 6a, b show the variation in wSTm and P on such paths.

At low fO2 (ΔFMQ ≲ +0.6), if the independent variables are T, fO2 and μFeS, the silicate melt is S2−-dominated and saturated with vapour and sulfide melt (i.e. the continuous grey sulfm curves in the purple–blue–turquoise regions of Figs 35 and 6a, b). Therefore, wSTm is controlled by S2−CSS from equations (9) and (14a); hence
wSTmwS2m=wS2CSSm.
(17)
As a result, S2−CSS is subparallel to the log10[wSTm]3 contour (Figs 35), although there is a small increase in wS2CSSm with increasing fO2 because: (1) there is always some S present in the silicate melt as S6+ and (2) wS2CSSm depends on P and [Fe3+/FeT]m (O'Neill 2021), both of which increase with increasing fO2. Nevertheless, wSTm is nearly constant and therefore to a good approximation P (and the pi's) behave as for the case when TfO2wSTm are independent variables (Figs 2c, d and 4).
At intermediate fO2 (+0.6 ≲ ΔFMQ ≲ +1.5), the independent variables remain the same (Y=μFeS; the silicate melt is saturated with vapour and sulfide melt and wS2m=wS2CSSm), but the silicate melt contains significant quantities of both S2− and S6+ (i.e. the continuous grey curves in the green region of Figs 35 and 6a, b). As the solubility of S6+ increases, wSTm increases based on equation (10) (e.g. in Fig. 4c where the nearly horizontal grey sulfm curve turns steeply upward as it enters the green region with increasing fO2). With increasing fO2, wSTm increases until the silicate melt reaches anhydrite saturation when wS6+m=wS6+CASm from equation (11). At this point, the silicate melt is saturated with v + sulfm + anh, and the independent variables are T, μFeS and μCaSO4. It is important to emphasize that for this choice of independent variables, at any given T, the vapour-saturated silicate melt defines a unique point in Figures 3 and 5 (the grey star at the intersection of the continuous and dashed grey curves) and fixes the values of all other intensive parameters. Therefore, combining equations (2), (9), and (11)
wSTm=wS2CSSm+wS6+CASm
(18rma)
and substituting equations (9) and (11) into (6c):
wS6+CASmwS2CSSm=CS6+CS2(fO2)2.
(18rmb)
This corresponds to the maximum wSTm at ΔFMQ  +1.5 for our calculations (i.e. the star on the grey curves in Figs 3c, 5c and 6a (labelled ‘SSVSAmax’); this multiply saturated point is not shown in Fig. 4 because it occurs above the wSTm range of the figure). Unlike the SSmin, the [S6+/ST]m of the SSVSAmax depends on melt composition and T (i.e. S2−CSS and S6+CSS). If the sulfur content of the system were above this value, it would go into additional vapour, sulfide melt and/or anhydrite. Hence, the wSTm at the SSVSAmax is the maximum wSTm for a particular composition of the silicate component at a given T. A maximum in P also occurs at this fO2 (Fig. 6a) because a higher P is required to keep the total S content dissolved in the silicate melt as described in the section ‘Regions with mixed speciation in the silicate melt or vapour’.
At higher fO2 levels than those described in the previous paragraph (i.e. ΔFMQ > +1.5), the vapour-saturated silicate melt can no longer be stably saturated with sulfide melt. However, it is still saturated with anhydrite (independent variables are T, fO2 and μCaSO4) and thus follows the dashed grey anh curve. The silicate melt is S6+-dominated and wS6+m is fixed at S6+CAS (which is taken to be a constant); hence, wS2m is given by equation (6c) rather than S2−CSS. Therefore, the amount of S2− in the silicate melt decreases whereas wS6+m remains constant. Hence, wSTm initially decreases with increasing fO2 based on equation (12), resulting in the maximum at the star (this small effect is most visible in Fig. 6b). At sufficiently high fO2 (ΔFMQ > +2.0), wSTm is essentially constant at S6+CAS from equations (11) and (14b):
wSTmwS6+m=wS6+CASm
(19)
and independent of P and fO2 given the parameterization we have used. Hence, wS6+CASm is essentially parallel to the log10[wSTm]4 contour over most of the yellow regions in Figures 3c and 5c, except for a small (but effectively undetectable) deviation because some S2− is always present, which decreases with increasing fO2.

The final multiply saturated assemblage is silm + sulfm + anh (no vapour is present). Because in our treatment the silicate melt is confined to the model ternary system and it is vapour-undersaturated, the silicate melt is the only stable phase. Therefore, φ = 1 and F = 4, meaning that four independent variables are needed to specify the state of the system. If the silicate melt is both sulfide melt- and anhydrite-saturated, only two of the three other variables (T, fO2 and P) can be independent (the other two being dependent variables). Choosing T, P, μFeS and μCaSO4 as the independent variables, fO2 can be calculated from equation (18b) and wSTm from equation (18a). Increasing P from the grey star, wSTm and fO2 decrease (dotted grey curve in Fig. 6a and b).

Variations in wSTm along isothermal and isobaric paths of increasing fO2

For Y values above the grey curves in Figures 35, the stable phase assemblage is no longer silm + v. The contours for vapour-saturated silicate melt shown beyond these curves are metastable in these regions. Above these grey curves, silicate melt is saturated with sulfide melt and/or anhydrite. For given values of T, fO2 and P (and parameterizations of S2−CSS, S6+CAS, CS2 and CS6+), wSTm at sulfide melt-saturation can be calculated using equation (10) and at anhydrite-saturation using equation (12). In this section, we choose independent variables of T, fO2, P and saturation with at least one of vapour, sulfide melt (μFeS) or anhydrite (μCaSO4) and calculate wSTm. The stable saturation condition (i.e. vapour-, sulfide melt- or anhydrite-saturation) is the one with the lowest wSTm at a specified T, fO2 and P (see Fig. S4 in the Supplementary material).

Figure 5f shows in PfO2 space (at T = 1200°C) the phase(s) with which the silicate melt is saturated. The dark and light grey regions are for sulfide melt- and anhydrite-saturation, respectively, and the coloured regions are for vapour-saturation as described in the section ‘Silicate melt + vapour’. We next explore how wSTm changes along isobaric paths of increasing fO2 at various P (the horizontal blue lines in Fig. 5g and shown in Fig. 7a–c). At P = 1543 bar the silicate melt is multiply saturated with vapour, sulfide melt and anhydrite (i.e. the grey star in Fig. 5f and g).

At P > 1543 bar, the silicate melt is vapour-undersaturated but the trend of wSTm with increasing fO2 (the horizontal blue dashed line in Fig. 5g and the blue curve in Fig. 7c) is similar in shape to the case where the silicate melt is saturated with vapour in addition to sulfide melt and/or anhydrite as described in the section ‘Multiply saturated silicate melt with vapour, sulfide melt and/or anhydrite’ (shown in Fig. 6b). The trends are similar because wSTm varies little with P when the silicate melt is saturated with sulfide melt and/or anhydrite. A maximum in wSTm is reached when the vapour-undersaturated silicate melt is multiply saturated with sulfm + anh (i.e. the grey square labelled ‘SSSAmax’ in Fig. 7c). Both μFeS and μCaSO4 are specified in addition to P and T at this point, and therefore fO2 is fixed for this set of four independent variables. The maximum wSTm without vapour but with sulfm + anh saturation (the grey square in Fig. 7c) is lower than with v + sulfm + anh (the star in Fig. 6b), demonstrating that the latter is the maximum wSTm for a particular composition of the silicate component at a fixed T.

At P < 1543 bar, the stable isobaric path crosses into the vapour-saturated field at intermediate fO2 (the horizontal blue dot–dashed and dotted line in Fig. 5g and the blue curves in Fig. 7a and b). At low fO2 (ΔFMQ ≲ +0.8), when the silicate melt is only saturated with sulfide melt, the behaviour is the same as at higher P (i.e. compare the blue continuous curves in Fig. 7a–c: these curves are nearly identical as the dependence of wS2CSSm on P is minor; O'Neill 2021). At P = 200 bar, a local maximum (ΔFMQ  +0.8) occurs when the silicate melt is saturated with both sulfide melt and vapour (i.e. the grey circle labelled ‘SSVSmax’ in Fig. 7a). When fO2 is increased beyond this point (+0.8 ≲ ΔFMQ ≲ +3.6), the silicate melt is only vapour-saturated, and wSTm displays the SSmin as described in the section ‘Regions with mixed speciation in the silicate melt or vapour’ (i.e. labelled ‘SSmin’ on the dot–dash black curve in Fig. 7a). However, at P = 800 bar, there is no local maximum then minimum in wSTm (Fig. 7b). This is because the fO2 of the SSmin is lower than the fO2 of silm + v + sulfm at this P, causing the SSVSmax and SSmin (like those in Fig. 7a) to disappear. For both P = 200 and 800 bar, the silicate melt is saturated with vapour and anhydrite when wS6+m=wS6+CASm (e.g. ΔFMQ  +3.6 or +2.2; the grey diamond in Fig. 7a and b). There is a small maximum at anhydrite and vapour saturation (at the grey diamond labelled ‘SSVAmax’), but it cannot be distinguished in Figure 7a and b because for the isobars chosen it occurs at sufficiently high fO2 that there is essentially only S6+ (i.e. S6+S2) dissolved in the silicate melt. Therefore, although wS2m decreases with increasing fO2 beyond this point, it is already so low that the magnitude of the decrease in wSTm is insignificant. At fO2 values above this point (ΔFMQ ≳ +2.2 or +3.6 at P = 800 and 200 bar, respectively), the silicate melt is only saturated with anhydrite (i.e. the dashed blue curve in Fig. 7a and b) and behaves as described previously under these conditions.

The examples developed in the section ‘Independent variables of T, fO2, fS2, wSTm, P, μFeS and μCaSO4’ and Figures 27 do not apply to the fO2fS2 paths followed by magmas as they degas on decompression (or as a result of crystallization). Modelling of degassing requires calculation of varying phase proportions (i.e. the amounts of vapour and silicate melt) as degassing proceeds. As mentioned in the section ‘The importance of the phase rule in our treatment of S-solubility in silicate melt’, our calculations to this point only constrain the states of the coexisting phases, not their amounts. To extend our modelling to degassing, we use Duhem's theorem, which states that for a closed system (i.e. where the masses of all components remain constant), if the masses of all the components are known, the equilibrium state is completely determined by two independent variables, which can be intensive or extensive (e.g. Prigogine and Defay 1954). Hence, we can simulate closed-system degassing if we have specified the composition of the system (i.e. the fixed amounts of the three components in our model ternary system) and any two independent variables. For the calculations presented in this section, we chose a constant value of T and variable P as the two independent variables (importantly, fO2 is not an independent variable in these calculations). We can also model open-system degassing as a series of increments of closed-system degassing where after each increment, the vapour is removed from the system.

Using this approach, we evaluate whether the SSmin has any impact on closed- and open-system degassing during isothermal decompression using our simple model ternary system. Inclusion of other volatiles (e.g. H and C) would add complexity to the system and is important for getting the details correct for simulation of degassing in nature (e.g. Moretti et al. 2003; Moretti and Papale 2004; Burgisser and Scaillet 2007; Gaillard and Scaillet 2009, 2014; Gaillard et al. 2011, 2015; Wallace and Edmonds 2011; Burgisser et al. 2015; Iacovino 2015; Liggins et al. 2020). However, their exclusion here allows us to isolate the behaviour of sulfur during degassing.

We model closed-system degassing by decreasing P at constant T (1200°C) and constant bulk composition for the two-phase system silm + v; hence, the independent variables are T and P at fixed bulk composition. The bulk composition comprises the proportions of the silicate, S2 and O2 components. The amount of the O2 component in the system is the O2 in excess of the oxygen contained in the model silicate component. A key difference between the degassing calculations and those presented in the section ‘Independent variables of T, fO2, fS2, wSTm, P, μFeS and μCaSO4’ is that at equilibrium mass balance must be satisfied; that is, the weighted sum of the concentrations of each component in the silicate melt and in the vapour must equal the specified bulk composition of the system. The silicate component is only present in the silicate melt phase. The S2 component is mass balanced based on the concentrations of S2− and S6+ in the silicate melt and S2 and SO2 in the vapour. The O2 component is mass balanced based on the O2 and SO2 in the vapor and SO42 plus the excess O2 in the silicate melt (relative to the silicate component end-member), which is partitioned between FeO and Fe2O3 in the silicate melt. For open-system degassing, any vapour present at each P step (1 bar increments) is removed from the system. Hence, the bulk composition of the system changes during depressurization and the composition of the silicate melt phase (including any dissolved S and excess O2) becomes the bulk composition of the system for the subsequent P step. We do not allow sulfide melt or anhydrite to be present; hence, in some regions the silm + v assemblage encountered during degassing is metastable (to the left of the continuous grey curve in Fig. 8f). Details for how these systems are solved are given in the Supplementary material. In Figure 8, the system contains 5000 ppm of the S2 component and the amount of the O2 component varies from 3.5 to 4.0 wt%, shown along the x-axis (the amount of the silicate component covaries such that the system sums to 100 wt%) Results of closed- and open-system degassing calculations are essentially indistinguishable (Fig. 9), hence we only discuss closed-system degassing calculations.

For a particular bulk composition, when P>Psatv (pressure of vapour-saturation; above the continuous black curve in Fig. 8), the system is vapour-undersaturated and only silicate melt (equal in composition to the bulk composition of the system) is present. At P=Psatv (on the continuous black curve), the system is still 100% silicate melt, but the silicate melt is vapour-saturated (i.e. equation (7) is satisfied using equations (3b), (4c), (5c), (8) and the bulk composition of the system). Hence, the 5000 ppm wSTm contour is the locus of Psatv (i.e. the P at which degassing begins; Fig. 8c) for a silicate melt corresponding to the bulk composition of the system. The white diamond indicates SSminfO2 (i.e. ΔFMQ + 1.30 for these parameters, at the maximum of the black curve). When P<Psatv (below the continuous black curve), the silicate melt is supersaturated with respect to vapour, such that vapour exsolves from the silicate melt and the two phases coexist under equilibrium conditions.

If the silicate melt starts off more reduced than SSminfO2, fO2 decreases with decreasing P during closed-system degassing (e.g. the yellow curves in Figs 8d and 9b). This reflects that although the silicate melt is S2−-dominated, the vapour is S2-dominated and/or SO2-dominated (Fig. 8f), so O2 is consumed during degassing to convert S2− dissolved in the silicate melt into the more oxidized S2 and SO2 vapour species (see reaction (4a) and the reaction generated by subtracting (3a) from (4a)). Therefore, reduction of the system (i.e. a decrease in fO2) must occur upon degassing under these conditions, although this effect is partially buffered by changes in [Fe3+/FeT]m, which are included in our calculations. However, if the silicate melt starts off more oxidized than SSminfO2, fO2 increases with decreasing P during closed-system degassing (e.g. the orange curve in Figs 8d and 9b). This reflects that the silicate melt is S6+-dominated but the vapour is SO2-dominated (Fig. 8f), so O2 is generated during degassing to convert SO42 dissolved in the silicate melt into the more reduced SO2 vapour species (see the reaction generated by subtracting (3a) from (5a)). Hence, oxidation of the system (i.e. an increase in fO2) must occur upon degassing under these conditions. Close inspection of Figure 8d demonstrates that for a narrow range of bulk compositions that start degassing close to (but slightly more reducing than, SSminfO2) fO2 initially increases and then decreases at lower P. In contrast to bulk compositions that are displaced in fO2 from SSminfO2 by more than ∼ 0.1 log units, there is minimal change in fO2 during depressurization for bulk compositions that start close to SSminfO2 (e.g. the green paths in Figs 8d and 9b). These closed-system degassing paths are also shown in part (e) of Figures 35.

At constant bulk composition when P<Psatv, wSTm decreases monotonically with decreasing P to 1 bar (Figs 8c and 9c). Hence, the SSmin is not encountered in our simple system during closed-system, depressurization-induced degassing because changes in fO2 do not cause the silicate melt to cross the SSmin. However, as illustrated in Figure 9c, the shapes of the wSTm with decreasing P are different at the SSmin (roughly linear) compared with those to either side (concave down). With decreasing P, fS2 (Figs 8a and 9a) and fSO2 (Fig. 8b) always decrease; the rate of decrease depends on the bulk composition (i.e. the O2 content) of the system. The rate of change is controlled by the trade-off between decreasing wSTm and decreasing or increasing fO2 to the left and right of SSmin.

The SSmin for silm + v and the SSmax for silm ± v + sulfm + anh (SSVSAmax and SSSAmax) occur in a similar fO2 range, as both features depend on the silicate melt having mixed S speciation (see the section ‘Independent variables of T, fO2, fS2, wSTm, P, μFeS and μCaSO4’). The SSmin manifests when T (held constant) and fO2 are independent variables, and fS2, wSTm or P is the third independent variable and is held constant or varied in particular ways (see the section ‘Regions with mixed speciation in the silicate melt or vapour’ and Figs 25). The SSmin reflects the prominent ‘valley’, centred where the silicate melt has mixed S speciation, in the isothermal surfaces showing wSTm as a function of fS2 or P v. fO2 (Figs 2a, c, 3c and 5c). An alternative expression of the SSmin is the prominent ‘ridge’ in the topography of the isothermal P surface plotted as a function of wSTm v. fO2 (Fig. 4c). However, there are simple linear paths that cross the valley (or ridge) with increasing fO2 in Figures 35c for which there is no SSmin. There are also geologically important paths (including closed- or open-system depressurization; see the section ‘Isothermal, decompression-induced degassing’) that nearly parallel the valley floor (or ridge crest) and therefore do no encounter an SSmin. Various SSmax occur along isothermal and isobaric paths when silm + v is additionally saturated with either sulfide melt or anhydrite (SSVSmax and SSVAmax). Summaries of the various SSmin and SSmax are shown in Figure 10.

The two geologically relevant paths have independent variables of T, fO2, wSTm or P. The SSmin for these paths occurs at [S6+/ST]m 0.75, resulting in an SSminfO2ΔFMQ+1.3 for the Hawaiian basaltic composition we have used at 1200°C. For comparison, this is more oxidized than most Hawaiian basalts (ΔFMQ −0.5 to +1.0; Moussallam et al. 2016; Brounce et al. 2017; Lerner et al. 2021) and mid-ocean ridge basalts (MORB), but within the range measured in arc and some other ocean island basalts (OIB) (e.g. Cottrell et al. 2021). From equation (6c), the position of SSminfO2 depends on CS2 and CS6+ and therefore will depend on T, composition of the silicate component and probably on additional volatile components (e.g. H2O, CO2). However, numerous experimental and modelling studies have demonstrated that the SSmin exists within the fO2 range relevant to common terrestrial magmatic environments (e.g. Fincham and Richardson 1954; Katsura and Nagashima 1974; Carroll and Rutherford 1985; Moretti et al. 2003; Clemente et al. 2004; Backnaes and Deubener 2011; Lesne et al. 2015; Matjuschkin et al. 2016; Nash et al. 2019). When these results are combined with our model results, it is likely that this phenomenon could influence terrestrial magmatic and volcanic processes.

In this final section, we explore several possible effects of the SSmin and SSmax in experiments (see the section ‘Solubility experiments’) and in natural processes and describe its possible application as a tool for constraining magmatic fO2 (see the section ‘Using wSTm as an oxybarometer’). These potential impacts on natural magmatic systems are summarized schematically in Figure 1c. We highlight natural processes throughout the magmatic and volcanic system, starting deep with mantle melting (see the section ‘Mantle melting’), up through the crust where magmas undergo mixing (see the section ‘Magma mixing and crustal assimilation’) and degassing (see the section ‘Magma ascent and degassing’), and finally to volcanic emissions into the atmosphere (see the section ‘Volcanic emissions’).

Solubility experiments

The SSmin and SSmax described in the section ‘Independent variables of T, fO2, fS2, wSTm, P, μFeS and μCaSO4’ have been observed in experiments that externally controlled P, T, fO2 and fS2. For example, the SSmin has been observed using 1 atm gas-mixing experiments at constant T and P with varying fO2 (e.g. Fincham and Richardson 1954; Katsura and Nagashima 1974). In these experiments, fS2 was not constant as the volume percentage of SO2 in the gas prior to heating in the gas mixing furnace was controlled instead. SO2 was mixed with other gas species and then heated to attain the required fO2. Despite the resulting complexity of gas speciation of theose experiments, the wSTm v. fO2 variations (including the presence of an SSmin) are reproduced by our modelling (e.g. compare fig. 2a of Fincham and Richardson (1954) with our Fig. 2e). Consistent with our analysis (see the section ‘Regions with mixed speciation in the silicate melt or vapour’), Fincham and Richardson (1954) attributed the first change in slope in their figure 2a to a change in vapour speciation (corresponding to the blue bar in our Fig. 2e) and the second to change in silicate melt speciation (SSmin: corresponding to the green bar in our Fig. 2e). At constant T and P, a different fS2fO2 path can be followed in 1 atm gas-mixing experiments using the bulk composition of the input gas to control fS2 as well as fO2 (e.g. O'Neill and Mavrogenes 2002; Nash et al. 2019). O'Neill and Mavrogenes (2002) did not observe an SSmin because the fO2 was always lower than the expected value of SStminfO2, but Nash et al. (2019) did observe an SSmin.

Piston cylinder experiments by Matjuschkin et al. (2016) produced an SSmin in sulfide melt + vapour-saturated silicate melts at constant P and T, in which fO2 was controlled using solid-state buffers (their fig. 8b). Their experiments were always saturated with sulfide melt (except one set at sufficiently high fO2 that were anhydrite-saturated), whereas the minimum described in the section ‘Variations in wSTm along isothermal and isobaric paths of increasing fO2’ (SSmin in Fig. 6d) is only vapour-saturated (i.e. it is not sulfide melt-saturated). Matjuschkin et al. (2016) attributed their SSmin to the presence of an additional S-bearing melt species with intermediate charge between 2- and 6+ that was not quenchable. Alternatively, their minimum might be influenced by a decrease in S2−CSS reflecting compositional changes in the silicate melt owing to crystallization in their experiments as fO2 varied.

A variety of studies using high-pressure, high-temperature apparatus at constant T and P and varying fO2 have observed an increase in wSTm as the system underwent a transition from sulfide melt- to anhydrite-saturation (e.g. Carroll and Rutherford 1985, 1987; Jugo et al. 2004; Beermann et al. 2011; Botcharnikov et al. 2011). However, these studies suggested a plateau rather than the SSmax seen in our model results (Fig. 6c), probably owing to the narrow fO2 range of the SSmax and/or the small elevation of wSTm above the S6+CAS plateau predicted by our modelling. However, Jugo (2009) modelled the data from Carroll and Rutherford (1985, 1987) and Jugo et al. (2004) and produced an SSmax similar to what we have shown in Figures 6b and 7c (SSVSAmax and SSSAmax although Jugo (2009) did not consider the presence or absence of vapour-saturation as an additional factor). It is encouraging that our calculated trends (Fig. 6c) match figure 1 from Jugo (2009), despite using different sets of equations, solubility mechanisms and thermochemical parameters.

Mantle melting

The maximum S content for silicate melts generated by partial mantle melting of sulfide- and anhydrite-bearing mantle is defined by the SSmax occurring when silicate melt is saturated with sulfm + anh ± v at a given T (SSVSAmax and SSSAmax). After either sulfide melt or anhydrite is exhausted during melting, the S content of the partial melt decreases to wSTCASm or wSTCSSm, respectively, and then decreases further by dilution when both phases are exhausted (e.g. Chowdhury and Dasgupta 2019). Chowdhury and Dasgupta (2019) explored the S content of silicate melts generated by mantle melting when the mantle contains either sulfide or anhydrite. They found that the S content of most arc magmas was equal to or less than wS2CSSm. Hence, the S content could be generated from mantle that initially contained sulfide if S2− was the dominate silicate melt species. For some arc magmas, the S content was higher than wS2CSSm but similar to the S6+ concentration expected after anhydrite is exhausted during melting but they calculated that sulfide should still be present in the mantle source during melting. They attributed these S contents to the presence of significant dissolved sulfate species in the silicate melt. This can be visualized in Figure 7c by partial mantle melting producing silicate melts within the medium-grey bar (i.e. saturated with sulfide melt or anhydrite and with both S2− and S6+ in the silicate melt). For two of their arc magmas, the S concentrations were too high for even anhydrite-saturated melting, which Chowdhury and Dasgupta (2019) suggested requires an additional S source, such as crustal assimilation of sulfate. Alternatively, our modelling suggests that these anomalously S-rich magmas could be generated by melting of mantle sources containing both sulfide and anhydrite (or with only one of these phases but at an fO2 close to SSmax), as this would result in a wSTm of the partial silicate melt higher than both S6+CAS and S2−CSS. A thorough assessment on mantle melting requires exploring how T and melt composition (including H- and C-bearing components) affect the calculated S contents, as these parameters influence CS2, CS6+, S2−CSS and S6+CAS (e.g. Moretti and Ottonello 2005; Chowdhury and Dasgupta 2019; Nash et al. 2019; O'Neill and Mavrogenes 2019; Zajacz and Tsay 2019; O'Neill 2021; Boulliung and Wood 2022). The presence of mixed S speciation in the melt also suggests that partial melting of sulfide- and/or anhydrite-bearing mantle could generate a wide range of S concentrations in the silicate melt (e.g. Jugo 2009).

Magma mixing and crustal assimilation

In nature, approximately constant P paths could be important when mixing reduced and oxidized S-bearing silicate melts. For example, mixing vapour-undersaturated silicate melts from either side of SSminfO2 could generate a mixed melt that is vapour-supersaturated if the S content of the mixed melt were greater than wSTm at the intermediate fO2. This scenario can be visualized using Figure 8. At P = 550 bar (black dashed horizontal line), both a relatively reduced (white circle at 3.26 wt% O2) and oxidized (white square at 4.23 wt% O2) silicate melt containing 5000 ppm wSTm would be vapour-undersaturated (Fig. 8c). Mixing these silicate melts isobarically results in a new bulk composition between the two end-members that depends on the mixing proportions. If the mixed silicate melt has a bulk composition inside the vapour-saturation curve for 5000 ppm ST (i.e. inside the black curve; between +1.0 ≲ ΔFMQ ≲ +1.7 for the example shown in Fig. 8), the mixed silicate melt will be vapour-supersaturated. The mixed magma will then degas until wSTm decreases to the contour at that point (Fig. 8c) and the amount of vapour degassed will depend on the mixing proportions (Fig. 8e).

This simple analysis assumes the oxidized and reduced silicate melts are the same in composition (other than in total O2). If they differed in major element chemistry or T, there would be added complexity because CS2, CS6+ and T will change depending on the proportions of the two silicate melts in the mixture (and the possibility of crystallization of the mixture). These factors would thus influence the solubility of S as a function of the mixing proportions of the two silicate melts. This would be particularly important if a reduced silicate melt assimilated oxidized country rock or an oxidized silicate melt assimilated reduced country rock (e.g. Tomkins et al. 2012; Iacono-Marziano et al. 2017). Nevertheless, encountering SSmin during mixing could be relevant to eruptive dynamics and volcanic SO2 contributions to the atmosphere by rapidly producing a large amount of vapour (e.g. Fig. 8e), potentially driving eruption (e.g. Kress 1997; Di Muro et al. 2008). It could also lead to deeper degassing, wringing out S-rich gas at higher P than would occur if the magmas in the mixture had erupted without mixing (see the section ‘Magma ascent and degassing’).

Magma ascent and degassing

For ascent of vapour-saturated magma toward the surface, P and bulk composition are two key independent variables. Our calculations for the model ternary system show that fO2 paths during isothermal depressurization are unlikely to cross SSminfO2 (Fig. 8d and the green path in Fig. 9b). Moreover, wSTm decreases monotonically with decreasing P for a closed (and open) system once vapour-saturation is reached (Figs 8c, 9c). In this simple system at low P, this is expected because the partial molar volumes of the gaseous S-bearing species are higher than those of the dissolved silicate melt species. Hence, the progressively more degassed (i.e. higher volume) state is the more stable one with decreasing P. We have not evaluated the possible effects of changes in T or silicate melt composition during degassing (e.g. owing to crystallization resulting from heat loss and/or from an increase in the liquidus on degassing). These changes would result in changes in CS2 and/or CS6+ and therefore the value of SSminfO2 at a particular P could vary during degassing, possibly resulting in the SSmin having more or less of an influence on the degassing process.

Degassing of C- and H-bearing species from a silicate melt can cause fO2 to increase (e.g. Sato and Wright 1966; Sato 1978; Mathez 1984; Candela 1986; Holloway 2004; Burgisser and Scaillet 2007; Brounce et al. 2017; Métrich 2021). The magnitude of this fO2 increase is greater when the initial fO2 of the magma is lower and depends on the relative solubilities of oxidized and reduced C- and H-bearing silicate melt species (e.g. Gaillard et al. 2015). Therefore, it may be possible to have a vapour-saturated silicate melt that starts S2−-dominated but oxidizes sufficiently owing to early CO2 and/or H2O degassing to drive the system toward or even across the SSmin during depressurization. If so, this could be manifested by initially decreasing, followed by increasing, wSTm with progressive closed-system, depressurization-induced degassing. It is also possible that a maximum in wSTm could occur during degassing, although this would not be due to changes in fO2 or the SSmin. For example, the loss of large enough quantities of CO2 and H2O prior to significant S degassing could cause wSTm to initially increase (i.e. because the total mass of the silicate melt is decreasing whereas the mass of dissolved S is nearly constant); then, when S begins to degas, wSTm would decrease.

Progressive reduction or oxidation of magmas during degassing is expected based on previous modelling efforts and has been observed in natural samples (e.g. Anderson and Wright 1972; Candela 1986; Carmichael and Ghiorso 1986; Burgisser and Scaillet 2007; Métrich et al. 2009; Gaillard et al. 2011, 2015; Kelley and Cottrell 2012; Moussallam et al. 2016, 2014; Brounce et al. 2017). Our degassing calculations suggest that the fO2 after extensive S-degassing will not represent the initial silicate melt unless the silicate melt begins degassing near SSmin (Fig. 8d). For example, when fO2 differs from that of SSminfO2 (i.e. at the diamond) by more than ∼ 0.1 log unit (either positively or negatively), the fO2 after nearly complete degassing (i.e. P = 1 bar) has increased or decreased by more than 0.5 log units (and up to > 1 log unit) relative to the fO2 of the silicate melt at initial vapour-saturation (Fig. 8d). The magnitude of the change in fO2 will depend on the initial wSTm (i.e. larger changes in fO2 occur for higher initial wSTm). Therefore, as has been highlighted previously, using the fO2 of volcanic gases (or the Fe3+/FeT or S6+/ST of silicate glasses) as a proxy for the fO2 of the initial silicate melt (and potentially of the mantle), should be approached with caution (e.g. Anderson and Wright 1972; Carmichael and Ghiorso 1986; Burgisser and Scaillet 2007; Métrich et al. 2009; Gaillard et al. 2015; Brounce et al. 2017). The key points to emphasize are that (1) the direction of the change in fO2 on degassing differs for more oxidized and reduced magmas relative to SSminfO2 and (2) silicate melts that begin to degas near SSminfO2 can maintain nearly constant fO2 during degassing (Figs 8d and 9b; e.g. potentially many arc and OIB magmas, although the effects of other volatiles, such as H and C, will need to be investigated).

The unlikelihood of crossing the SSmin during closed- or open-system degassing does not mean that SSmin is unimportant during degassing. Decompression paths to the left and right of SSmin begin to degas at lower P than one that passes right through the maximum (Fig. 8d). The shapes of the wSTm contours in Figure 8c are simple expressions of the SSmin for the case Y=wSTm (Figs 2d and 4c). The maximum in P at the SSmin occurs because a higher P is required to keep a given concentration of S dissolved in the silicate melt (5000 ppm in Fig. 8c) at the SSmin than when the S solubility is higher on either side. Based on our modelling, a Hawaiian basaltic melt with 5000 ppm ST reaches vapour-saturation at 666 bar at SSmin (the green curve in Fig. 9). This is 164 bar higher than when the melt is initially 0.5 log units higher in fO2 (the orange curve) and 363 bar higher than when the melt is initially 0.6 log units lower in fO2 (the yellow curve). Thus, SSmin can exert significant control on the P at which S degassing begins as a function of fO2 (Fig. 8). For instance, silicate melts with initial fO2 close to SSminfO2 (e.g. many arc and OIB magmas) would, all other things being equal, begin to degas S deeper than more oxidized or reduced magmas. Therefore, the SSmin could be a contributing mechanism, in addition to S partitioning into the H2O-rich vapour exsolved from H2O-rich magmas (e.g. Wallace and Edmonds 2011; Zajacz et al. 2012; Edmonds and Mather 2017; Edmonds and Woods 2018), for the deep degassing of S in arc and some ocean island settings.

The SSmin also has implications for calculating the Psatv using the volatile concentrations of glasses (e.g. in pillow rims and glassy melt inclusions) (e.g. Anderson et al. 1989; Blundy and Cashman 2008; Hughes et al. 2023). Such calculated pressures are often used to constrain the architecture of magmatic systems and link volcanic products to eruptive vents in submarine systems. These calculations are based on generalizations of equation (7) that include other gaseous species for which partial pressures can be determined (e.g. CO2 and H2O). This exercise is comparable with when T, fO2 and Y=wSTm are taken as independent variables in a silm + v assemblage and the dependent variable Z = P is calculated (e.g. Fig. 4c). The effect of dissolved S is currently not included in calculations of Psatv, which include only H2O and CO2 (e.g. Newman and Lowenstern 2002; Papale et al. 2006; Iacono-Marziano et al. 2012; Ghiorso and Gualda 2015; Allison et al. 2019; Iacovino et al. 2021). Also, the effect fO2 on volatile speciation in the silicate melt and vapour is mostly ignored (but see Scaillet and Pichavant 2004; Wetzel et al. 2015).

MORB magmas have [S6+/ST]m < 0.05 (e.g. Métrich et al. 2009; Jugo et al. 2010; Labidi et al. 2012), which is much lower than the [S6+/ST]m 0.75 of the SSmin. Therefore, including S for such melts would have a negligible effect on calculated Psatv because pS2 and pSO2 would be very low (e.g. Figs 2d and 10a). However, OIB and arc volcanic glasses can have [S6+/ST]m up to ∼ 1 (e.g. Jugo et al. 2010; de Moor et al. 2013; Labidi et al. 2015; Brounce et al. 2017; Muth and Wallace 2021), and therefore the effects of the SSmin on the contribution of partial pressures from S-bearing species to the total P at vapour-saturation (e.g. S2 and SO2 in the S–O system, but potentially H2S and OCS as well in C–O–H–S system) could be non-negligible (e.g. in the green bar of Fig. 2d). Melt inclusions from arc and OIB magmas reach wmST ∼ 5000 ppm (Wallace 2005) and have fO2 levels of ΔFMQ2 to +2 (Cottrell et al. 2021), which could lead to underestimates of Psatv by up to 650 bar for melts at ΔFMQ = +1 to +2 for the example we have considered (Fig. 10a). The exact value of the underestimate would depend on the fO2, CS6+ and CS2 (and hence T and silicate melt composition), and on including the effects of H and C species (e.g. Burgisser et al. 2015; Lesne et al. 2015; Hughes et al. 2023). However, the potential magnitude of this SSmin-related effect implied by our calculations is robust. Therefore, in such cases neglecting S and the SSmin could result in calculated Psatv values that significantly underestimate the Psatv the glass is recording.

We note that the maximum in Psatv at SSmin based on our calculations contrasts with the modelling of Lesne et al. (2015), who predicted a minimum in Psatv when both reduced and oxidized S species are present in silicate melts. The difference could be due to (1) different assumptions regarding the speciation of oxidized S in the silicate melt (i.e. as SO2 rather than SO42; see the section ‘Equilibria between silicate melt, vapour, sulfide melt and anhydrite’ for a listing of potential S-bearing species in the silicate melt) and (2) the effect of H on the SSmin, which adds H2S species to the silicate melt and vapour.

Volcanic emissions

The El Chichón (1982) and Pinatubo (1991) eruptions released the largest quantities of SO2 recorded by satellites during explosive events, and most of this SO2 was sourced from a coexisting vapour already present prior to eruption (e.g. Wallace and Gerlach 1994; Krueger et al. 1995; Gerlach et al. 1996; Keppler 1999; Bluth et al. 2015). Additionally, both magmas contain anhydrite and pyrrhotite/other sulfides in the erupted products (e.g. Luhr et al. 1984; Bernard et al. 1991; Luhr 2008). We infer from this observation that both these magmas were stored at SSVSAmax prior to eruption (i.e. at the condition represented by the grey star in Fig. 5g). Independent fO2 estimates from these magmas are within the range for experimental constraints for the coexistence of sulfide melt and anhydrite (e.g. Luhr et al. 1984; Rye et al. 1984; Carroll and Rutherford 1987; Evans and Scaillet 1997; Luhr 2008). Older eruptions have released even greater quantities of SO2 (e.g. fig. 5 of Vidal et al. 2016). For example, eruptive products from the Samalas 1257 eruption contain sulfides and vapour, and anhydrite has been observed as microcrystals on the walls of fluid inclusions (Vidal et al. 2016; it should be noted that anhydrite is highly soluble in water and is therefore rarely observed in volcanic products; e.g. Luhr et al. 1984). The presence of sulfides + anhydrite in the products of explosive eruptions associated with unusually large releases of SO2 is consistent with these magmas having evolved under conditions near the SSVSAmax. If so, it may suggest a connection between the fO2 range at which the SSVSAmax occurs and these events.

The ‘petrological method’ is often used to estimate the volatile emissions from volcanic eruptions by subtracting the volatile concentration of degassed matrix glass from that in melt inclusions (the latter is assumed to represent undegassed silicate melt; e.g. Devine et al. 1984; Thordarson et al. 1996; Wallace 2001; Sharma et al. 2004). The total SO2 emitted during eruption estimated in this way is often low relative to those measured using other techniques (e.g. satellite-based techniques; e.g. Stoiber and Jepsen 1973; Rose et al. 1982; Andres et al. 1991; Wallace 2001; Shinohara 2008). This ‘excess sulfur’ problem reflects the large amounts of SO2-rich vapour (although still dominated by H2O and CO2) often present during magma storage (i.e. only minor amounts of SO2 are thought to be released from the breakdown of sulfide melt and anhydrite during degassing; e.g. Anderson 1975; Luhr et al. 1984; Andres et al. 1991; Gerlach et al. 1994, 1996; Gerlach and McGee 1994; Wallace and Gerlach 1994; Giggenbach 1996; Keppler 1999; Wallace 2001; Scaillet and Pichavant 2003; Scaillet et al. 2003; Sharma et al. 2004; Shinohara 2008). One possibility is that the amount of S-rich vapour is enhanced by the magma being close to the SSmin (and to some extent to the SSVSAmax). Therefore, using melt inclusions that formed at SSmin or SSVSAmax and trapped only silicate melt will result in large discrepancies with other methods owing to the additional S in other phases. Melt inclusions that co-entrap these additional phases cannot be used because the proportions of the different phases present in the inclusion are unlikely to represent the bulk system. Additionally, as melt inclusions evolve as a closed system to S post-entrapment, SSmin and various SSmax may be encountered during cooling, causing additional phases to form within the inclusion (e.g. vapour bubbles, sulfide blebs). Measuring only the silicate glass within these melt inclusion results in a greater underestimate of SO2 emissions because the S contribution from the silicate melt is additionally underestimated (e.g. Venugopal et al. 2020).

Using wSTm as an oxybarometer

For a silicate melt with a given value of wSTm and T, it is possible to place constraints on its fO2 using the calculations we have presented based on the presence or absence of vapour, sulfide melt and/or anhydrite as saturating phases (e.g. Beermann et al. 2011; Muth and Wallace 2022). Let us suppose it is known that a silicate melt is vapour-saturated. If it is also saturated with sulfide melt it must fall on the continuous grey curve in Figure 10b, so for a known wSTm, the fO2 can be read directly from this curve (e.g. α and β for 5000 and 13 000 ppm wSTm, respectively, in Fig. 10b). Alternatively, if the silicate melt is anhydrite-saturated it must fall on the dashed grey curve in Figure 10b; the implied fO2 for a given value of wSTm can be read from this curve (e.g. γ for 13 000 ppm wSTm in Fig. 10b). If the silicate melt is saturated with both sulfide melt and anhydrite, the fO2 is constrained to the SSmax value at the grey star.

Let us suppose, however, that we do not know whether the vapour-saturated silicate melt is sulfide melt and/or anhydrite saturated (or we know that it is not). Then the fO2 can be constrained to be between the grey curves in Figure 10b for a given value of wSTm (e.g. at 13 000 ppm wSTm, the fO2 must be between β and γ in Fig. 10b). When nearly all of the dissolved S is S2− (i.e. in the purple–blue–turquoise region), the continuous grey curve reaches a plateau at ∼ wS2CSSm, so for values of wSTm<wS2CSSm, no constraint can be placed on a lower bound to fO2 using this approach. Likewise, the dashed grey curve reaches a plateau at ∼ wS6+CASm, and thus no constraint on an upper bound on fO2 can be determined by this approach if wSTm<wS6+CASm. If the silicate melt is not thought to be vapour-saturated, but P is known independently, a figure like Figure 7 at the relevant P can be used instead in the same way as described for the vapour-saturated case. Although the presence of other volatiles (e.g. C and H) will modify the results from the simple S2–O2 system, the principles are the same.

The technique described using Figure 10b is most sensitive when measured S concentrations are high and the S speciation in the silicate melt is mixed (ΔFMQ ∼ +1 to +2); for example, to obtain both lower and upper bounds on fO2 in the example described requires wSTm>wS6+CASm13000ppm (i.e. the black dot–dash line in Fig. 10b). Although this is much higher than the typical total dissolved S contents of common magmas, some studies suggest that S6+CAS decreases significantly with increasing dissolved H2O and decreasing T; for example, S6+CAS is ∼ 3000 ppm at 5 wt% H2O and 1200°C for a basaltic melt (Chowdhury and Dasgupta 2019). Melt inclusions from arcs and ocean islands can be S-rich (up to ∼ 7000 ppm ST), hydrous (up to ∼ 6 wt% H2O for arcs and ∼ 3 wt% H2O for ocean islands) and relatively oxidized (up to ΔFMQ + 3) (e.g. Wallace 2005; Moussallam et al. 2019; Cottrell et al. 2021) and may therefore provide useful fO2 estimates based on this technique. However, this will require accurate knowledge of CS2, CS6+, S2−CSS and S6+CSS at the relevant conditions (especially the effects of T and silicate melt composition, including the influence of H2O), as these parameters strongly influence the fO2 of the transition from S2−- to S6+-dominated silicate melt (e.g. O'Neill and Mavrogenes 2002, 2019; Li and Ripley 2005; Moretti and Ottonello 2005; Baker and Moretti 2011; Chowdhury and Dasgupta 2019; Nash et al. 2019; Zajacz and Tsay 2019; O'Neill 2021; Boulliung and Wood 2022).

Depending on the choice of independent variables, vapour-saturated silicate melts can with increasing fO2 encounter a ‘sulfur solubility minimum’ (SSmin) when both S2− and S6+ are dissolved in the silicate melt in similar concentrations (Figs 25). This minimum occurs because O2 is on different sides of the reactions describing S2− and S6+ dissolution in the silicate melt from S2- and/or SO2-dominant vapour (e.g. reactions (4a) and (5a)). Examples of choices of independent variables and the paths they follow that exhibit a minimum in the dissolved total S content (wSTm) in vapour-saturated silicate melt include isothermal paths of increasing fO2 on which either fS2 or P is held constant (Figs 3c and 5c). For isothermal paths on which wSTm is held constant with increasing fO2, the SSmin is expressed as a maximum in P (Fig. 4c). The SSmin is at [S6+/ST]m = 0.25 for constant fS2, but [S6+/ST]m 0.75 for constant wSTm or P. However, not all choices of independent variables or paths defined by changes in these variables display the SSmin. An important geological example of this is that despite changing fO2, wSTm decreases monotonically (i.e. no minimum in wSTm is encountered) during isothermal, closed-system, decompression-induced degassing in a system in which the vapour contains only S- and O-bearing species (Fig. 8c).

Depending on the isothermal path, there are also maxima in wSTm (SSmax) of silicate melts that are multiply saturated with two of vapour, sulfide melt and anhydrite. An absolute maximum wSTm (for a given T and choice of silicate melt composition) occurs if all four phases are present (SSVSAmax). This is a univariant assemblage, so once the independent conditions of constant T and the coexistence of silm + v + sulfm (which sets the value of μFeS) + anh (which sets μCaSO4) are imposed, P and fO2 (and all other intensive parameters) are fixed (this maximum in wSTm occurs at a single point, the grey star, in Figs 3 and 5). As for the SSmin, the SSVSAmax occurs at an fO2 value at which both S2− and S6+ are dissolved in the silicate melt in similar concentrations (Fig. 6b). This maximum can be explained by the constraints of sulfide melt- and anhydrite-saturation leading to simultaneous maximization of the concentrations of both the S2− and S6+ species in the silicate melt. A maximum in wSTm is also encountered at constant T with increasing fO2 when a vapour-undersaturated silicate melt is both sulfide melt- and/or anhydrite-saturated (SSSAmax in Fig. 7c), also when S2− and S6+ are dissolved in the silicate melt in similar concentrations.

These SSmin and SSmax features can play roles over the entire magmatic and volcanic system, extending from the mantle to eruption (Fig. 1c). However, their influences depend on the independent variables governing the system at each point during the process, the paths followed by these variables, and their effects on various dependent variables. For example, these features can influence the maximum S concentration in mantle melts, volatile release owing to magma mixing and crustal assimilation, the depth at which significant amounts of S begin to degas from silicate melt, the fO2 path followed by melt and volcanic gases during progressive degassing of magmas, the fO2 of erupted magma and emissions of volcanic gases, and the amounts of SO2 released to the atmosphere during explosive eruptions. Of particular interest is the result that melts with initial fO2 values lower than the SSmin will decrease in fO2 with progressive degassing, those with higher initial fO2 values will increase in fO2 with degassing and those at the SSmin will not change in fO2 on degassing. This could have significant implications for understanding measured fO2 values of erupted melts from arc and ocean island magmatic systems, as their fO2 levels extend from lower to higher than the position of the SSmin modelled here. Additionally, SSmin and SSmax may also affect some of the tools used to infer intensive and extensive variables of these systems. Excluding the effects of the SSmin (and SSmax) can lead to significant underestimations in calculations of Psatv, as well as SO2 emissions using the petrological method. SSmax also provides the possibility of constraining fO2 for S-bearing magmas based on limits set by the fO2 dependence of wSTm of silicate melts saturated with vapour, sulfide melt and/or anhydrite.

We emphasize that our approach has been to use a simplified ternary system (silicate–O2–S2) to model thermodynamically the coexistence of silm + (v and/or sulfm and/or anh). This choice allows us to isolate and analyse the interplay of key variables in a system in which all the volatile species in the vapour are on the S2–O2 join. Therefore, we have not included the effects of other volatile components that typically make up most of the gas phase. Although we are confident that the patterns and behaviour expressed in this simple system generalize to natural systems, an important next step will be to include other components and species in the silicate melt and vapour. Our approach can be readily expanded to model such complex natural systems, particularly when H and C are present. Principally, this involves including additional species and homogeneous equilibria to the vapour (e.g. H2, H2O, CO, CO2, CH4, H2S, OCS, etc.) and solubility reactions for the species that dissolve in the silicate melt (e.g. OH, H2O, H2, CO32−, CO2, CO, CH4, H2S, SH, etc.). Full generalization of our results and approach to natural systems will require exploring the effects of variations in composition on solubility and speciation (e.g. the composition of the silicate melt, including H2O and [Fe3+/FeT]m) and the P and T dependence of the equilibrium constants for heterogenous silicate melt–vapour equilibria. This is currently of particular importance for the sulfate capacity (CS6+), which controls the behaviour of dissolved sulfate and is expected to be strongly influenced by T and melt composition (e.g. Moretti and Ottonello 2005; Nash et al. 2019; O'Neill and Mavrogenes 2022; Boulliung and Wood 2022; Moretti 2021). Although such an expanded treatment will be important and useful in detailed modelling of natural systems, the complexity already present in the simple system we have examined (i.e. a single Hawaiian basaltic melt composition in which only S ± O-bearing species are included) highlights in our view the importance of adding complexity (especially in terms of silicate melt and vapour composition) incrementally to such end-member systems for understanding the behaviour of sulfur in magmatic systems.

We would like to thank A. Saal and R. Moretti for their detailed and thoughtful reviews, which greatly improved our work, and G. Shields for editorial handling of the paper.

ECH: conceptualization (equal), formal analysis (equal), investigation (lead), methodology (equal), software (lead), writing – original draft (lead), writing – review & editing (equal); LMS: conceptualization (equal), methodology (supporting), writing – review & editing (supporting); PL: software (supporting), writing – review & editing (supporting); HSCO: formal analysis (supporting), writing – review & editing (supporting); EMS: conceptualization (equal), formal analysis (equal), methodology (equal), writing – original draft (supporting), writing – review & editing (equal)

ECH was funded by a Caltech Geology Option Post-Doctoral Fellowship and a Caltech Centre for Comparative Planetary Evolution (3CPE) research grant, and is supported by the Hazards and Risk Management Programme, which is part of New Zealand Strategic Science Investment Funding (SSIF) from the New Zealand Ministry of Business, Innovation & Employment (MBIE). PL is funded by an Embiricos Trust scholarship from Jesus College, University of Cambridge.

The authors declare no known conflicts of interest associated with this publication.

All data generated or analysed during this study are included in this published article (and its supplementary information files). The code used to generate the data is available at https://github.com/eryhughes/SSminmax.

This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (http://creativecommons.org/licenses/by/4.0/)