PySulfSat : An open-source Python3 tool for modeling sulﬁde and sulfate saturation

We present PySulfSat , an Open-Source Python3 tool for modeling sulﬁde and anhydrite saturation in magmas. PySulfSat supports a variety of input data types (spreadsheets, Petrolog3 outputs, MELTS tbl ﬁles), and can be directly integrated with alphaMELTS for Python infrastructure to track sulfur solubility during fractional crystallization within a single Jupyter Notebook. PySulfSat allows easy propagation of uncertainty using Monte Carlo methods, and far more customization of calculations than existing tools. For example, the SCSS 2 − could be calculated with one model using the sulﬁde composition from a parameterization released with a different SCSS 2 − model. There are also functions for calculating the proportion of S 6 + /S Tot (allowing modeled SCSS and SCAS values to be converted into total S solubility to compare to natural data), and for modeling mantle melting in the presence of sulﬁdes using a variety of SCSS and K D models. Extensive documentation and worked examples are available at ReadTheDocs ( https://bit.ly/PySulfSatRTD ) along with narrated YouTube videos ( https://bit.ly/PySulfSatYouTube ).


INTRODUCTION
Modeling the solubility of sulfur in a silicate melt provides vital insights into the evolution of sulfur and other S-loving (chalcophile) elements during mantle melting and crustal processes such as fractional crystallization and crustal contamination [Ding and Dasgupta 2018;Reekie et al. 2019;Wieser et al. 2020;Wieser and Jenner 2021;Iacono-Marziano et al. 2022;Muth and Wallace 2022;Virtanen et al. 2022]. Modeling the removal of sulfide and sulfate phases is particularly vital to understand the formation of economical deposits of chalcophile elements, as well as the sulfur and metal flux emitted to the atmosphere during volcanic eruptions [Edmonds et al. 2018;Wieser et al. 2020;Mason et al. 2021]. A number of different models have been proposed over the years to calculate the sulfide content at sulfide saturation (SCSS 2− ), which describes the amount of sulfide (S 2− ) that can dissolve in a silicate melt saturated in a sulfide phase [e.g. Li and Ripley 2009;Fortin et al. 2015;Smythe et al. 2017;O'Neill 2021]. Numerous models also exist to quantify the sulfate content at anhydrite saturation (SCAS), which describes the amount of sulfate (S 6+ ) that dissolves in a silicate melt when saturated in anhydrite [e.g. Li and Ripley 2009;Baker and Moretti 2011;Masotta and Keppler 2015;Chowdhury and Dasgupta 2019;Zajacz and Tsay 2019]. In many magmas with intermediate oxygen fugacity (e.g. in volcanic arcs), S is present as a mixture of S 2− and S 6+ species . O'Neill and Mavrogenes [2022], Nash et al. [2019], and Jugo et al. [2010] present models to quantify the proportion of these two species as a function of melt redox. These speciation models can be used alongside SCSS 2− and SCAS 6+ calculations to obtain the total amount of S that is dissolved in the melt (to compare to measured S contents in volcanic systems). * penny_wieser@berkeley.edu

Previously-available tools
At the moment, SCSS 2− and SCAS 6+ calculations are performed in spreadsheets accompanying each publication [e.g. Fortin et al. 2015;Smythe et al. 2017;O'Neill 2021]. These spreadsheets have a limited number of rows for performing calculations (e.g. = 50 for Smythe et al. [2017], = 194 for O'Neill [2021), making it difficult to apply them to thousands of natural compositions, or outputs of fractional crystallization models with a small temperature step. The prevalence of Excel-based tools also makes it difficult to propagate uncertainty using Monte Carlo methods. Available tools also make it time consuming and difficult to compare different models. Existing spreadsheets require users to paste in their melt compositions with oxides in a specific order, and the order differs between spreadsheets. After reformatting the input structure for each model, users would then have to extract outputs and compile these into a single format and location for plotting. There are also tools for which no published spreadsheets exist [e.g. Blanchard et al. 2021], requiring users to contact the author team to obtain such a tool, or individually interpret the equations (which may contain typographical errors, or ambiguities, particularly regarding which units to use).
The most recent SCSS 2− models have a term accounting for the composition of the sulfide [Smythe et al. 2017;Blanchard et al. 2021;Liu et al. 2021;O'Neill 2021;Li and Zhang 2022], because melts in equilibrium with a sulfide containing Ni and Cu have a substantially lower SCSS compared with melts in equilibrium with pure Fe-S sulfides. However, the spreadsheets for these different models use a variety of approaches to account for the composition of the sulfide, making it hard to directly compare model outputs. The Smythe et al. [2017] Excel workbook has two sheets; one is designed for users to enter a sulfide composition in wt.%, while the other sheet calculates a sulfide composition using partition coeffi-We include numerous narrated worked examples on the PySulfSat YouTube channel to make this package more accessible to non coders * . Some relevant terminology for Python and S modeling is shown in Table 1. * https://bit.ly/PySulfSatYouTube Certain models also require users to input the following parameters ( Figure 1): The import_data function returns a pandas dataframe (see Table 1). The order of the columns in the input spreadsheet does not matter, as columns are identified based on their column heading rather than position. If any column headings are missing in the input spreadsheet, they will be filled with zeros. Any additional columns entered by the user (e.g. temperature, pressure, sulfide composition) are appended onto the end of the outputted dataframe, for easy access for calculations. For example, the O'Neill [2021] and Smythe et al. [2017] models require the Ni and Cu content of the liquid in ppm. These can be stored in a column with any heading the user wishes (e.g. Ni_Liq_ppm, Cu_Liq_ppm), and then obtained from the outputted dataframe (df) using df ['column_name'] to input into the function of interest. For example, to import generic data (perhaps whole-rock, matrix glass or melt inclusion compositions) from a spreadsheet named "Liquids1.xlsx" stored in "Sheet3": This function also supports specific output files from other petrological modelling programs.
For example, users can load in the default spreadsheet-based output from Petrolog3.1.1.3 [Danyushevsky and Plechov 2011]. Here, the Petrolog output is saved to an Excel file named "Petrolog_Model1.xlsx": Similarly, the standard liquid ".tbl" output from MELTS [Ghiorso and Sack 1995;Asimow and Ghiorso 1998;Gualda et al. 2012] can be imported: df_out=ss.import_data (filename='melts-liquid.tbl', MELTS=True) In these examples, the import_data function has identified the appropriate column headings in each default structure, and has changed the column names into the format required by PySulfSat (e.g. converting SiO2_melt from Petrolog3 into SiO2_Liq ).

Units
All temperatures should be entered in Kelvin, all pressures in kbar, and all melt oxides in wt.%, apart from Ni and Cu contents in the liquid which are entered in ppm. All ratios are atomic (e.g. Fe/(Fe+Ni+Cu) in the sulfide).

Available functions
PySulfSat implements the most recent SCSS 2− and SCAS 6+ models ( Figure 1). The open-source nature of PySulfSat means we anticipate continuining to add models as they are published, so users should check the 'Available Functions' tab at ReadTheDocs * . * https://bit.ly/PySulfSatRTD

Calibration datasets
Many SCSS and SCAS models are empirical. Thus, it is not recommended that they are extrapolated too far beyond the compositional range of the calibration dataset. We have compiled available calibration datasets, and incorporated them into PySulfSat (see Figure 1 for available datasets). This means that users can easily plot their melt compositions, and estimates of the pressures and temperatures of their system alongside the dataset used to calibrate each model, to assess its suitability. The function return_cali_dataset returns the calibration dataset for a given model. For example, to obtain the calibration dataset for the Smythe et al. [2017] SCSS model as a pandas.DataFrame: df_S2017=ss.return_cali_dataset(model='S2017_SCSS') Figure 2 shows how these different calibration datasets can be plotted in TAS (total alkali silica) space for visual inspection.

Worked examples
Example Jupyter Notebooks showing a number of workflows are available at ReadTheDocs page † . This list is not exhaustive, and we anticipate that we will continue adding examples in the future: • Notebooks showing how to import different data types (e.g. measured oxide contents, Petrolog3 files, and MELTS tbl outputs).
• Notebooks showing how to calculate the SCSS and SCAS using a variety of models during fractional crystallization from a Petrolog3 output [Danyushevsky and Plechov 2011]. This example also shows how to calculate the trajectory of S if a sulfide phase was not present, and how to calculate the mass fraction of sulfide which has formed during crystallization.
• Notebooks showing how to run a MELTS fractional crystallization paths at a single pressure and at multiple pressures using PyMELTScalc [Gleeson et al. 2023], and then calculate the SCSS and SCAS within the same Jupyter Notebook.
• Notebooks showing how to model the SCSS from a Petrolog3 path, and compare models of S contents and sulfide composition to natural melt inclusion and sulfide data.  Li and Zhang (2022) "calculate_LiZhang2022_SCSS"

Other func�ons
"crystallize_S_incomp" Calculates S le� in the melt for a given Fmelt (assuming S is en�rely incompa�ble "calculate_mass_frac_sulf" Calculates mass frac�on of sulfide removed for a model of changes in SCSS with frac�onal crystalliza�on "convert_d34_to_3432S", "convert_3432S_to_d34"  FeO T Liq (wt%) Figure 2: Plots of SCAS calibration datasets in P-T-X space. An example notebook to produce these plots and overlay user data is available at ReadTheDocs. Similar plots can easily be made for SCSS models.
• Notebooks showing other useful features, including calculating values of K using various models, converting between S isotope ratios and delta notation, and abundances of different S-bearing species.

SCSS 2− MODELS
There are a number of ways to perform SCSS calculations, with various options discussed below (worked examples are available at ReadTheDocs).

Using measured sulfide compositions
The newest SCSS models [e.g. Smythe et al. 2017;Blanchard et al. 2021;O'Neill 2021;Li and Zhang 2022] contain terms for the composition of the sulfide. In some situations, the sulfide composition may have been directly measured in the samples of interest (e.g. using Energy Dispersive Spectroscopy ). If so, the function calculate_sulf_FeFeNiCu can be used to convert measured elemental abundances in wt.% into the atomic Fe/(Fe+Ni+Cu) ratio used by SCSS models. In some systems, the Fe/(Fe+Ni+Cu) may remain approximately constant during fractional crystallization , meaning that a fixed value for this ratio can be used for simplicity. Figure 3 shows a worked example calculating the SCSS 2− using the models of Smythe et al. [2017], O'Neill [2021], and Li and Zhang [2022] for Fe/(Fe+Ni+Cu)=0.65. The expected increase in the S content of the melt with fractional crystallization in the absence of a S-bearing phase is also calculated using the function crystallize_S_incomp for comparison (black dashes), and these different S trajectories are plotted using matplotlib (where they can be compared to natural melt inclusion or quenched submarine glass data).

Calculating sulfide compositions
While using a measured sulfide composition is the simplest and most reliable method to perform SCSS 2− calculations, direct measurements of sulfide compositions do not exist in many systems. PySulfSat allows users to calculate sulfide composition from Ni and Cu contents of the liquid using the approaches implemented in the supporting spreadsheets of O'Neill [2021] and Smythe et al. [2017]. The O'Neill [2021] method is the simplest, calculating the atomic Fe/(Fe+Ni+Cu) ratio using the following empirical expression: where: If the sulfide composition is not known, the spreadsheet of Smythe et al.
[2017] has a sheet which will iteratively calculate the sulfide composition based on the partition coefficients of Cu and Ni in the sulfide from Kiseeva and Wood [2015]. These partition coefficients are sensitive to temperature, liquid FeO content, and the Ni and Cu content of the sulfide. Starting with a first estimate of the sulfide Ni and Cu content, the temperature, and the FeO content of the liquid, a partition coefficient can be calculated. Using this partition coefficient along with the initial estimate of the Ni and Cu content in the sulfide, the amount of Cu and Ni in a melt in equilibrium with this sulfide can be calculated. Smythe et al. [2017] define a residual between this calculated value and the measured Ni and Cu contents of the melt: The Excel solver function varies the Ni and Cu in the sulfide to obtain the values which best minimise this residual. Then, the equation of Kiseeva and Wood [2015] is used to calculate the Fe content of the sulfide for these best fit sulfide Ni and Cu contents, and these 3 parameters are used to calculate the sulfide Fe/(Fe+Ni+Cu) ratio. In PySulfSat, this convergence routine is performed using the scipy optimize minimize function . In Excel, for many compositions, the result obtained can depend slightly on the starting value of the Ni and Cu contents in the sulfide provided by the user. By default, the PySulfSat minimisation starts with initial Ni and Cu contents of 5 wt.%, but these parameters can be overwritten using Cu_Sulf_init=10 and Ni_Sulf_init=5. These parameters are allowed to vary between 0-30 wt.%. In general, we find our python implementation of this solver method is stable and gives identical results to the Excel version for the same starting composition (and the vast majority of samples converge regardless of the starting Ni and Cu contents).
To perform SCSS calculations with modeled sulfide compositions, a string should be entered into the Fe_FeNiCu_Sulf argument. For example, to use the Smythe et al. [2017] SCSS 2− model with the O'Neill [2021] calculated sulfide composition, enter Fe_FeNiCu_sulf='Calc_ONeill'. Users must also specify the Cu and Ni content in the liquid. In the example below, Ni_Liq (ppm) and Cu_Liq (ppm) are columns in the loaded dataframe df_out containing estimated Ni and Cu contents of the melt in ppm:  Blanchard et al. [2021], and Li and Zhang [2022] are sensitive to the amount of H 2 O in the liquid. By default, the SCSS 2− functions for each of these models (Figure 1) use the H 2 O content stored in the data loaded by the user in the column H2O_Liq. However, this can also be overwritten in the function itself, to allow investigation of the sensitivity of calculations to melt water content. For example, to perform all calculation at 3 wt% H 2 O using the Fortin et al. [2015] model: The argument H2O_Liq could also be set to a pandas series (e.g. any other column in the loaded data), which would allow calculations to be performed using several different water contents (e.g. df_out ['Raman_H2O'] for Raman spectroscopy measurements vs. df_out ['SIMS_H2O'] for SIMS measurements in the same samples). Identical inputs to above, only difference is the function name! Scroll Figure 3: Annotated worked example showing how to calculate SCSS 2− for a Petrolog3 fractional crystallization path using a fixed Fe/(Fe+Ni+Cu) ratio in the sulfide. Hypothetical melt inclusion data is overlain. The data initially follows the incompatible fractional crystallization trend, followed by a prominent downturn, indicating the onset of sulfide saturation at ∼ 6-7 wt.% MgO.

Redox sensitivity
A number of SCSS models are also sensitive to the ratio of Fe 3+ , because they contain a term for only Fe 2+ species in the melt (see Figure 1). The input argument Fe3Fet_Liq should be supplied when using these models. If no value is entered, calculations are performed assuming Fe 3+ = 0. Alternatively, users can specify a single value in the function (e.g. Fe3Fet_Liq=0.15), or refer to a column in the input dataframe. Another option is to use the Python package Thermobar [Wieser et al. 2022] to convert a log O 2 value or buffer position into a Fe3Fet_Liq ratio. For models which are not redox-sensitive [e.g. Blanchard et al. 2021;Liu et al. 2021], entering a non-zero value for Fe3Fet_Liq will not affect the SCSS (except through secondary dependencies, e.g. if the model of Smythe et al. [2017] or O'Neill [2021] is used to calculate the sulfide composition).

Calculating sulfide proportions
The difference between the fractional crystallization trajectory and the predicted SCSS 2− can be used to calculate the cumulative mass proportion of sulfide forming over the fractionation interval [after Kiseeva and Wood 2015]: where S init is the initial S content at the start of the fractional crystallization sequence (F melt = 1), F melt is the melt fraction remaining at each step, S model is the modeled solubility of S in the melt, and S sulf is the S content of the sulfide (all concentrations in ppm). In PySulfSat, this is calculated as follows for the example shown in Figure 3: This calculates the mass fraction of sulfide formed for a magma with 800 ppm S initially, a S content in the sulfide of 32 wt.%, and a melt fraction from the Petrolog3 file (column heading Fraction_melt, obtained from the column Melt_%_magma in the Petrolog3 file by the PySulfSat import function). As for SCSS 2− models, these functions return the calculated SCAS 6+ , all intermediate calculations, and the originally loaded compositions. The main simplification relative to SCSS models is the fact that none of the existing SCAS models have a term for the composition of the sulfate-bearing phase, pressure, or the Fe 3+ /Fe T ratio (Figure 1).

MAGMAS WITH A MIX OF S 2− AND S 6+
Silicate melts undergo a relatively abrupt transition in S speciation from sulfide (S 2− ) to sulfate (S 6+ ) dominated with increasing oxygen fugacity [Fincham and Richardson 1954;Wallace and Carmichael 1994;Jugo et al. 2010;Kleinsasser et al. 2022, ; cyan line in Figure 4B]. In systems where both S 2− and S 6+ are present, the calculated SCSS 2− will underestimate the total solubility of S, because this parameter only accounts for the solubility of S 2− species. Similarly, in systems dominated by S 6+ with some S 2− , the total solubility of S will exceed the SCAS 6+ [Jugo 2009].

Demonstrating the importance of S 2− and S 6+ corrections
To demonstrate the importance of accounting for both S 2− and S 6+ species when modeling total S solubility, lets consider a melt with an SCSS 2− of 1000 ppm, and an SCAS 6+ of 5000 ppm. Equation 10 of Jugo et al. [2010] can be used to calculate the proportion of S 6+ /S T as a function of ΔQFM between −1 and +3: This equation can be implemented in PySulfSat for a single ΔQFM value as follows: To produce the cyan line on Figure 4B, we input a linearlyspaced numpy array of 10,001 points between ΔQFM = −1 and ΔQFM = 3 generated using the np.linspace function, and calculate S 6+ /S T for every value in this array (cyan line, Figure 4B). At ΔQFM = −1 (point 1 on Figure 4B), the melt is sufficiently reduced that only S 2− is dissolved in meaningful quantities (S 6+ /S T = 0.00008). Thus, the total solubility of sulfur is well approximated by the SCSS 2− (1000 ppm for this specific example, horizontal magenta line on Figure 4A).
For a moderately oxidized melt at ΔQFM = 1 (Point 2), S 6+ /S T = 0.442, so the presence of S 6+ species substantially increases the total amount of S that is dissolved. Thus, the SCSS 2− must be corrected to obtain the SCSS T using the equation of Jugo et al. [2010]: In PySulfSat this is implemented as follows: The SCSS T comprises with 1000 ppm of S 2− , and 794 ppm of S 6+ (see red and grey lines on Figure 4B). At ΔQFM = 1.39 (Point 3), S 6+ /S T = 0.827. Using Equation 6, the SCSS T is 5786 ppm, with 1000 ppm of S 2− , and 4786 ppm of S 6+ . However, if ΔQFM (and therefore S 6+ /S T ) increases slightly more, Equation 6 becomes invalid, because the amount of predicted S 6+ exceeds the SCAS 6+ (dashed magenta line, Figure 4A). For example, at point 4 (ΔQFM = 2), S 6+ /S T = 0.9876. Equation 6 would predict that the SCSS T is 80,433 ppm, with 1000 ppm of S 2− , and 79,433 ppm of S 6+ . However, this much S 6+ cannot dissolve, because the SCAS 6+ is only 5000 ppm.
Instead of correcting the SCSS for S 6+ , in more oxidising magmas, we can also correct the SCAS 6+ for the presence of S 2− : For example, at point 4: SCAS_Tot=ss.calculate_SCAS_Total(SCAS=5000, S6St_Liq=0.9876) =5063 This total dissolved S comprises 5000 ppm of S 6+ and 63 ppm of S 2− . However, if this equation were applied to point 2, it would predict more dissolved S 2− than the SCSS. These worked examples demonstrate that at certain proportions of S 6+ to S 2− , Equation 7 and Equation 6 are invalid for predicting the total solubility of S. For the specific SCSS 2− and SCAS 6+ values used in this example, ΔQFM = ∼1.4 is the oxygen fugacity where the maximum amount of S dissolves in the system, because at this ΔQFM value, the ratio of S 6+ /S T is such that the amount of S 2− dissolved is equal to the SCSS 2− , and the amount of S 6+ is equal to the SCAS 6+ (yielding the maximum possible sum of these two values). The total amount of dissolved S in ΔQFM space that does not violate the calculated SCSS 2− and SCAS 6+ is defined by the section of the SCSS T curve where S 6+ does not exceed the SCAS 6+ (magenta solid line, Figure 4A), and the section of the SCAS T curve where S 2− does not exceed the SCSS 2− (black solid line, Figure 4A). The combined curve meeting these requirements is shown as a green line in Figure 4B.
In PySulfSat, for any calculated SCSS 2− and SCAS 6+ values, the total amount of S can be calculated using the function calculate_S_Total_SCSS_SCAS. This can be used to produce plots of changing S speciation with fO 2 (e.g. Figure 4).
In  Figure 4: Calculating the total amount of dissolved S by applying corrections for the presence of both S species using the model of Jugo et al. [2010] in the function calculate_S_Total_SCSS_SCAS. These graphs were drawn for SCSS 2− = 1000 ppm and SCAS 6+ = 5000 ppm, although these numbers could be calculated using any SCSS and SCAS model in PySulfSat.

Calculating S 6+ /S T from the Sulfate and Sulfide capacity
In addition to the methods described above where the proportion of S species is estimated from oxygen fugacity or Fe 3+ /Fe T , the ratio of S 6+ /S T can also be calculated using the method of O'Neill and Mavrogenes [2022]. This approach calculates the sulfide capacity (C S 2− ) using the parameterization of O'Neill [2021], and the sulfate capacity (C S 6+ ) using O'Neill and Mavrogenes [2022]. The equilibrium constant for the gasphase equilibrium, lnK, is then calculated using in Kelvin: These values are then used to calculate S 6+ /S 2− , which can be easily converted into a S 6+ /S T ratio: and: The O'Neill and Mavrogenes [2022] supporting spreadsheet also provides an option to input Fe 3+ /Fe T ratio instead of a value for log O 2 . The spreadsheet uses this ratio to calculate ΔQFM using an adapted version of Eq9a of O'Neill et al. [2018] (missing the term for P 2 O 5 , as this oxide is not included in their capacity models): Where X Na , X K , and X Ca are the cation fractions of Na, K, log 10 fO 2 = ∆QFM − 25050/T + 8.58.
Where is in Kelvin.
These equations are all implemented in PySulfSat through the function calculate_OM2022_S6St. For example, to perform calculations using a known logf O 2 value: This allows direct comparison between models.
We also include the function calculate_BW2022_OM2022_S6St to calculate S 6+ /S T using C S 6+ from Boulliung and Wood [2022] and C S 2− from O'Neill [2021].

Calculations for natural samples
When calculating the total solubility of S in a natural system with a non negligible proportion of both S species, using the function calculate_S_Total_SCSS_SCAS ensures that the correction has not exceeded the solubility of either species, unlike functions correcting the SCSS for S 6+ using calculate_SCSS_Total, or SCAS for S 2− using calculate_SCAS_Total.
When comparing measured S contents to total S solubility obtained from SCSS and SCAS models, it is most reliable to use measured S 6+ /S T ratios (e.g. using X-ray absorption near edge structure-XANES [Lerner et al. 2021]). In this ideal scenario, users can enter the measured ratio directly in the calculate_S_Total_SCSS_SCAS function. For example, after calculating the SCSS using Smythe et al. [2017] (saved in df=S2017) and the SCAS using Zajacz and Tsay [2019] (saved in df=Z2019), the total amount of dissolved S can be calculated using a fixed S 6+ /S T ratio of 0.2: Alternatively, it is more common that Fe 3+ /Fe T has been constrained using XANES. Using the Nash et al. [2019] correction, this Fe 3+ /Fe T ratio can be entered directly to calculate the S 6+ /S T ratio, and thus the maximum amount of S that can dissolve: For consistency, in this example, the S2017 dataframe should also have been calculated using the same input value for To use the Jugo et al. [2010] correction, the redox state of the magma must be calculated relative to the QFM buffer position of Frost [1991] (see Section 4.2. If Fe 3+ /Fe T is known, this can be converted into a log O 2 value using Kress and Carmichael [1988] using the Python package Thermobar [Wieser et al. 2022]. Once a log O 2 value is calculated, Thermobar can then be used to calculate the offset from the QFM buffer position (i.e. ΔQFM). Alternatively, the log O 2 may be known independently without having to do a conversion based on Fe 3+ /Fe T first. For example, the Petrolog3 output in Figure 3 has a column for the log of the O 2 value, the temperature, and the pressure. The different buffers stored in the Buffer_calc dataframe can then be input into the PySulfSat function: This function can also take Fe 3+ /Fe T as input, although our code (and the published spreadsheet) convert this into a log( O 2 ) value using Equation 13.

Different buffer positions and melt redox models
It is important to recognize the uncertainty introduced into calculations of S 6+ proportions as a result of different definitions of buffer positions, melt redox models, and XANES data processing strategies. For example, the ΔQFM values for the Jugo et al. [2010] S 6+ correction should be relative to the QFM buffer position of Frost [1991]. Petrolog3 uses the expression of Myers and Eugster [1983] for its QFM buffer. alphaMELTS (including MELTS for MATLAB and Python [Antoshechkina and Ghiorso 2018]) and MELTS for Excel [Gualda and Ghiorso 2015] also use Myers and Eugster [1983], with an additional pressure correction from Frost [1991]. Expressing all these different QFM buffer positions in terms of log( O 2 ) values at QFM yields the following equations: log f 2 at QFM [Frost 1991 where is in bars and is in Kelvin. If users have a ΔQFM value relative to a buffer which is not Frost [1991], they need to convert that into a value relative to the Frost [1991] prior to using Jugo et al. [2010].  Sack et al. [1981], Kilinc et al. [1983], Kress and Carmichael [1988], and Borisov and Shapkin [1990]. For the default Petrolog3 composition at ΔQFM=0, atmospheric pressure and the liquidus position, these 4 models return Fe 3+ proportions between 10 and 14 %, which corresponds to S 6+ proportions using Nash et al. [2019] Brounce et al. [2017] and Muth and Wallace [2021] using the method of Berry et al. [2018] prior to performing calculations of S 6+ proportions. However, we find that the S XANES measurements of Muth and Wallace [2021] are best matched by the O'Neill and Mavrogenes [2022] if measured Fe 3+ /Fe T ratios are input, rather than ratios corrected using Berry et al. [2018]. More comparisons are clearly required to see if this is a one-off occurrence. In many instances, offsets between different Fe 3+ /Fe T XANES reduction methods, and Fe 3+ /Fe Tlog( O 2 ) conversion strategies should perhaps be considered as true error on these methods, given the lack of community consensus [Anenburg and O'Neill 2019].

MONTE CARLO ERROR PROPAGATION
In addition to simplifying calculations and aiding model comparisons, PySulfSat also allows users to propagate uncertainty in input parameters for all calculation types using Monte Carlo methods. There are two main workflows that can be used. First, if errors are known for every input variable, users should load in two dataframes. The first dataframe (df1) should contain the preferred value for each input parameter (e.g. columns MgO_Liq, FeOt_Liq, H2O_Liq). The second dataframe (df2) should have exactly the same column headings with the addition of the suffix _Err. These columns can contain absolute or percentage errors. Additional columns (e.g. temperatures calculated using Thermobar [Wieser et al. 2022]) can be appended onto df1 in the Jupyter Notebook itself, along with an appropriate error in df2. The function add_noise_2_dataframes can then be used to duplicate each input row in the input dataframe df_values N_dups times, adding noise based on the value in the dataframe with errors (df_err). For example, to add normally distributed errors using absolute 1σ values from df2, and create 5000 duplicates for each sample: df_noisy=ss.add_noise_2_dataframes( df_values=df1, df_err=df2, error_type="Abs", error_dist="normal", N_dups=5000) This new dataframe is then entered into any of the functions.
In Figure 5 we use the add_noise_2_dataframes function to generate 5000 synthetic compositions for each melt inclusion, with errors from quoted 1σ values for each variable from Muth and Wallace [2021]. These synthetic compositions were then input into the various functions to calculate S 6+ /S T ratios. Finally, the function av_noise_samples_series is used to calculate statistics for each melt inclusion. Users should input a panda.Series containing the variable of interest into this function as the arguement var (in this case, the calculated S 6+ /S T ratio stored in the dataframe ONeill_S6ST), and a second pandas.Series with the sample names to average over (as sampleID). For example, the first row in the output averages all 5000 simulations for rows with the sample name BBL-5-32.
This function calculates the mean, median, max, min and standard deviation of all 5000 simulations for each melt inclusion (which are used to plot symbols and error bars on Figure 5F-H). Simulations can be conducted for any of the calculations available in PySulfSat (e.g. SCSS, SCAS, K , etc.). A second set of functions can be useful when you want to explore noise in a smaller number of input variables (e.g. just , and H 2 O), or where some errors are absolute, some are percentages, some are normal and some are uniformly distributed. The function duplicate_dataframe takes a dataframe and duplicates the values in each row N_dup times (row1-row2-row3 goes to row1-row1-row1..., row2-row2-row2 ...): Then the function add_noise_series can be used to create a panda.Series of noise for one specific variable with the same length as this larger dataframe. For example, EDS measurements in a suite of lavas may reveal the sulfide composition for sample 1 is Fe/(Fe+Ni+Cu) = 0.65, and sample 2 is 0.8, with an error of ± 0.05 (stored in the column df1['Sulf_X']. Here, we add normally distributed noise, with 5000 duplicates for each input (to match the dataframe above).

Dupdf['Sulf_MC']=sulf_comp_err
As many 'noisy' columns can be added as the user wishes, with different error types and distributions. This dataframe where some columns have noise added and some do not can then be input into any of the PySulfSat functions.

INTEGRATION WITH MELTS
While PySulfSat can load the results from a MELTS calculation as a .tbl file, recent advances in the MELTS computing infrastructure means that MELTS fractional crystallization calculations can be performed directly in Python in the same Jupyter Notebook as PySulfSat calculations. There are currently two options for performing MELTS calculations in Python; Thermoengine [Johnson et al. 2022] and alphaMELTS for Python [Antoshechkina and Ghiorso 2018]. We make use of the PyMELTScalc Python package [Gleeson et al. 2023], which provides neatly wrapped functions for fractional crystallization using alphaMELTS for Python, and returns output structures consistent with the required inputs for PySulfSat.
After installing PyMELTScalc (see example on ReadThe-Docs), this package must be imported into the notebook: import pyMELTScalc as M After loading data using the ss.import_data function as df_out, a specific melt composition can be selected as a starting composition (here, we select the first row):

sample=df_out.iloc[0]
Then, a MELTS fractional crystallization model can be initiated at a single pressure using the multi_path function: This runs a fractional crystallization model at 1000 bars (P_bar), starting at the wet liquidus (find_liquidus=True), and runs until 750 • C (T_end_C). If the MELTS calculation does not converge after 100 quadratic minimisation attempts, the simulation may end at a higher temperature. The temperature step is 5 • C (dt_C), the initial Fe3Fet_Liq ratio is set at 0.1, and both fluids and solids are fractionated.
This multi_path function outputs a dictionary containing a series of dataframes. There is a dataframe for each phase, but most relevant for this work is the dataframe named 'All'. This contains all the relevant outputs stitched together, and can be obtained from the overall output as follows:

MELTS=MELTS_FC['All']
This dataframe named MELTS contains system properties ( , , enthalpy, entropy, volume) and the composition of each phase with the phase name as an underscore (e.g. SiO2_Liq, SiO2_Plag etc.). It can be fed directly into the PySulfSat code. For example, lets use the model of Li and Zhang [2022] for a specified sulfide composition: PyMELTScalc can also be used to investigate a wide range of different fractional crystallization paths using parallel processing for computational efficiency, with hundreds to thousands of different fractional paths initiated with a single function call. For example, coupling of PyMELTScalc and PySulfSat would allow users to investigate S behavior during fractional crystallization for a single melt or range of melt compositions over a wide variety of different starting pressures, oxygen fugacities, and melt water contents. Figure 6 shows the SCSS 2− calculated for fractional crystallization models run at 4 different pressures from a single call to the PyMELTScalc multi_path function. PyMELTScalc can run calculations at a redox buffer or unbuffered, so calculations can be implemented with the various options for the treatment of S 6+ to investigate changes in S speciation during fractional crystallization.

MANTLE MELTING CALCULATIONS
Modeling the concentrations of S, Cu and other chalcophile elements during mantle melting is complicated by the fact that these elements are held in silicate minerals and mantle sulfides. Because mantle melts contain higher S contents than the mantle residue, the mantle becomes more and more depleted in sulfide during progressive melting until the sulfide phase is eventually exhausted [Lee et al. 2012;Ding and Dasgupta 2018;Wieser et al. 2020]. Exhaustion of sulfide in the mantle residue drives a large change in the bulk partition coefficient of chalcophile elements during the melting interval. Lee et al. [2012] provide an Excel spreadsheet for calculating the concentration of Cu during near-fractional melting. This model removes small batch melts, updating the composition of the remaining mantle residue before the next melting step proceeds. The equation for batch melting is as follows: Where melt is the concentration in the melt, source is the concentration in the mantle source, D 0 is the bulk partition coefficient (sulfide+silicate) at the start of that melting step, F is the degree of melt produced in that melt step, and is the bulk partition coefficient weighted for the proportion that each component enters the melt. For simplicity, Lee et al. [2012] assume that 0 = (e.g. sulfide and silicate minerals melt at the same rate). Wieser et al. [2020] update this model to account for non-modal melting behavior, because the sulfide preferentially melts, so contributes more to the partition coefficient of highly chalcophile elements such as Cu than the silicates. It should be noted that at a small enough step size (i.e. small enough Δ ), the results from these two approaches converge. However, using the limited number of columns supplied in the spreadsheet of Lee et al. [2012], the divergence can be several 10s of ppm at a given extent of melting ( ).
We implement the non-modal melting version of Wieser et al. [2020] in PySulfSat with the function Lee_Wieser_sulfide_melting. This function can be used to model the concentration of any element during near fractional batch melting, and allows the contrasting behavior of chalcophile and lithophile elements to be modeled (e.g. Ba vs. Cu ). The user must supply a dataframe with partition coefficients for silicate and sulfide phases, and the mass proportion of each phase. In Figure 7A For simplicity in this example, we assume that the silicate modes stay fixed throughout the melting interval. This assumption makes very little difference for Cu, as the partition coefficient is substantially higher for sulfides than any silicate phases. Even for Ba, this is a reasonable first-order assumption because it is extremely incompatible in all high abundance silicate phases. The other required inputs are: 1. The number of iterative steps (N=3000) 2. The S content of the mantle source in ppm (S_Mantle=200) 3. The concentration of S in mantle sulfides in ppm (S_Sulf=360000) 4. The initial concentration of the element of interest in the mantle prior to melting in ppm (elem_Per=30) 5. The S 2− concentration of the melt in ppm (S_Melt_SCSS_2=1000).
6. The proportion of S 6+ (here Prop_S6=0), which will be used alongside the S 2− concentration to calculate the total amount of S in the melt using Equation 6: These inputs are then used as follows for Cu: These calculations were run at S_Mantle contents of 100 ppm, 200 ppm, and 300 ppm to produce Figure 7A-B).
In addition to the ease of the above calculations vs. existing tools, the other substantial advantage of PySulfSat is that it allows integration of melting models with models for partition coefficients in sulfides, and models of the SCSS within a single calculation environment. This enables a more sophisticated modeling approach than existing studies, which assumed a fixed SCSS throughout the melting interval [e.g. Lee et al. 2012;Wieser et al. 2020;Muth and Wallace 2022]. In reality, the major element composition of instantaneous melts will change as melting proceeds, particularly for incompatible elements such as Na 2 O and K 2 O. Consequently, the SCSS will change during melting, rather than being set at a fixed value. PySulfSat can be used to calculate the SCSS for instantaneous melt compositions from melting models. For example, the cyan line in Figure 7C-E shows calculations using instantaneous melt compositions estimated from a Thermocalc melting model [Jennings and Holland 2015]. This model using a calculated SCSS has a higher S content in the initial melts than the model assuming S = 1000 ppm throughout, resulting in a lower sulfide mode, a lower bulk K , and thus a higher Cu concentration in mantle melts at low F values (cyan vs. dashed magenta line, Figure 7). Sulfide is also exhausted at a lower (black star, Figure 7C). Changing silicate melt modes can also be used instead of a fixed modal abundance, which will create more realistic trajectories for elements with an affinity for both sulfide and silicate phases.
While both the cyan and magenta models on Figure 7 assume K for sulfide-melt is fixed at 800, PySulfSat can also be used to calculate K as a function of temperature, liquid FeO content, and the Ni and Cu content of the sulfide using the model of Kiseeva and Wood [2015]. This more rigorous K approach results in a substantially lower K , and thus higher Cu contents in the melt. Additional information on how to perform these more advanced melting calculations can be found at ReadTheDocs. Overall, PySulfSat gives substantially more flexibility to explore concentrations in instantaneous and aggregated melts for all elements during melting in the presence of sulfide phases.

OTHER USEFUL FUNCTIONS
We also include a number of functions for other common workflows associated with S. For example, the functions convert_d34_to_3432S and convert_3432S_to_d34 can be used to convert between δ 34 S values and 34 S/ 32 S ratios. By default, these functions use the the Vienna-CDT value of 1/22.6436 from Ding et al. [2001], although this can be overwritten with any value of interest (using the input st_ratio). For example, if a dataframe is loaded in with a column for d34S the isotope ratio can be calculated as follows: We also include a function which allows users to enter the amount of S present in the melt as S in wt.%, ppm, or as SO 2 , SO 3 , or SO 4 . It then converts this concentration into an equivalent concentration expressed as different species (useful when converting EPMA data measured as SO 2 into S in ppm for example):

df=ss.convert_S_types(S_ppm=df['S_ppm'])
Additionally, the studies of Kiseeva and Wood [2015] and Brenan [2015] parameterize K s as a function of melt composition, and sulfide composition for Kiseeva and Wood [2015]. The function calculate_sulfide_kds can be used to calculate these partition coefficients.

FUTURE WORK AND CITATION
The open-source nature of PySulfSat, along with recent increase in interest in the behavior of S in magmas, means that   et al. [2020] where the K in the sulfide, the modal proportion of silicate minerals and S in the melt is kept constant throughout the melting interval. Variation in elemental concentrations correlate with the initial S content of the mantle source.
[C]-[E] More complex models combining melting models with K D and SCSS functions within PySulfSat. For 200 ppm S in the mantle source, substantially different trajectories can be generated by varying the model for the amount of S in the melt, or the partition coefficient of Cu. The cyan and blue lines use a mantle melting model from Thermocalc to obtain the major element contents and temperature of instantaneous melts [Jennings and Holland 2015]. This allows the S content of these melts to be determined using the SCSS model of O'Neill [2021], assuming mantle sulfides contain 20 wt.% Ni and 5 wt.% Cu (after Ding and Dasgupta [2018]). The cyan line uses a fixed K D for Cu (800, after Lee et al. [2012]). The blue line uses K D calculated from the instantaneous silicate melt composition and an estimated mantle sulfide composition from Kiseeva and Wood [2015]. All models assume there is 30 ppm Cu in the mantle source.
this tool will continuously evolve. The current author team will endeavor to add new models as they are released, and anyone can submit new code using a pull request on GitHub (or by contacting the authors). Thus, users should check the ReadTheDocs page, where examples demonstrating new functionality beyond that described in this manuscript will be added in the future. New versions of PySulfSat can be obtained by running the following code in a Jupyter environment: !pip install PySulfSat --upgrade When citing calculations performed in PySulfSat in papers, users should be sure to specify which version they used, which can be obtained using:

ss.__version__
For example, the text may read "SCSS calculations were performed using the model of Smythe et al. [2017] implemented in PySulfSat v.1.0.3 (Wieser and Gleeson, 2023)." It is important to cite all the original papers used to perform calculations (e.g. the SCSS model, the model for S 6+ ), as well as citing PySulfSat.
At present, there is no open-source code that can model sulfide and sulfate saturation with all the most recent models, and the behavior of S during degassing from a a silicate melt. We hope that in future, the PySulfSat source code can be integrated with the wide variety of S degassing tools becoming available [e.g. Ding et al. 2023] to produce a single, coherent model engine for modeling S behavior in silicate melts.

REPORTING BUGS AND REQUESTING FEATURES
No software is free of bugs, particularly when new features are being constantly added. We have extensively benchmarked PySulfSat to existing spreadsheets, and before the package is published on PyPI, automatic unit tests are run through GitHub in the attempt to catch problems introduced by changing Python dependencies/updates. However, if users spot any bugs, or wish to request features, they should submit an 'issue' on the GitHub page. Alternatively, they can email the author team.

CONCLUSIONS
PySulfSat is a open-source Python3 tool motivated by the FAIR (Findable, Accessible, Interoperable, and Reusable) research framework [Bloemers and Montesanti 2020]. It will greatly speed up calculations, allow more intercomparison between models, and through its ease of implementation with Python, allow more detailed and robust investigations of the behavior of sulfur in magmatic systems (with a rigorous consideration of errors).

AUTHOR CONTRIBUTIONS
PW conceived the project, wrote the S-based code and the manuscript. MG built the fractional crystallization MELTS functions allowing integrating of pyMELTScalc with PySulfSat, and contributed to manuscript editing and code testing.