Understanding the long-term variations in seawater sulfate concentrations ([SO42−]sw) is crucial to our understanding of the dynamic relationships between the sulfur, carbon, calcium and oxygen cycles, and their influence on the habitability of the Earth. Here, we explore how [SO42−]sw has changed throughout the Phanerozoic and its impact on other elemental cycles. We do this by utilizing the biogeochemical box model GEOCARBSULFOR. The model suggests that [SO42−]sw increased throughout the Paleozoic, decreased during the Mesozoic and then increased once more in the Cenozoic, generally matching geochemical proxies. Atmospheric oxygen mirrors [SO42−]sw changes during the Paleozoic and Mesozoic, but, intriguingly, decouples during the Cenozoic. We further explored the controls on [SO42−]sw by modifying the modelled gypsum fluxes via the incorporation of evaporite data from the geological record. We found that forcing gypsum burial with the observed evaporite deposition data causes the model to better match proxy records at some times, but worsens predictions at others. We also investigated the reliance of the model on a prescribed record of marine calcium concentrations, finding that it is a dominant control on modelled Phanerozoic [SO42−]sw and that removing this control seriously degrades the model predictions. We conclude that no model can yet simulate a reasonable evolution of both the calcium and sulfur cycles.

Supplementary material: Figures S1–S5 are available at https://doi.org/10.6084/m9.figshare.c.7164928

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

Alongside the carbon cycle, the biogeochemical sulfur cycle has exerted a major influence on the Earth system through geological time (Holland 1973; Garrels and Perry 1974; Garrels and Lerman 1981; Berner and Raiswell 1983; Holser et al. 1988; Halevy et al. 2012; Wortmann and Paytan 2012; Fike et al. 2015). Sulfur is a fundamental component for many different forms of life, but is used dominantly by prokaryotes, which metabolize sulfur in its various oxidation states, often also oxidizing organic carbon when reducing sulfate to sulfide (e.g. Jørgensen 1982; Postgate 1984; Widdel 1988; Shen and Buick 2004; Newton and Bottrell 2007; Leavitt et al. 2013; Fike et al. 2015). Ultimately, in the presence of iron, the terminal product of sulfate reduction is pyrite and the subsequent long-term burial of this reduced mineral represents a major net source of oxygen to the atmosphere, whereas pyrite weathering represents a major sink (Garrels and Perry 1974; Berner 1982). Reduced sulfur compounds can also act as electron donors for photosynthesis by some prokaryotes (Bryant and Frigaard 2006).

In modern day settings, because sulfate reduction in marine sediments produces an excess of sulfide over the total available reactive iron, c. 80–95% of the sulfide initially produced by microbial sulfate reduction is reoxidized back to sulfate (with greater preservation occurring at greater depths), either abiotically or by microbial processes, rather than forming pyrite (Jørgensen 1982). Pyrite formation and preservation rates (and thus oxygen production and consumption) and the microbial influence on these reactions have waxed and waned through time (Canfield 2013), but nevertheless microbial processes have intimately coupled the sulfur, carbon and oxygen cycles on Earth on geological timescales (Garrels and Perry 1974; Garrels and Lerman 1981; Newton and Bottrell 2007; Fike et al. 2015).

The largest pool of biologically available sulfur in the oceans is in the form of dissolved sulfate, with a modern day concentration of c. 28–29 mM, equivalent to a total of 3.8 × 1019–3.9 × 1019 moles and a residence time of c. 12.5 Myr (Canfield 2004; Kah et al. 2004; Weldeghebriel et al. 2022). Over geological time, the major inputs of sulfate to the oceans have come from the following sources: the uplift and subsequent oxidation of sulfides (e.g. pyrite); the weathering of sulfates, such as those associated with evaporite deposits (e.g. gypsum), carbonate-associated sulfate (CAS) or barite; the thermal decomposition at depth of these same species, leading to volcanic degassing of sulfur; and hydrothermal inputs at mid-ocean ridge spreading centres, either from mantle degassing or the dissolution of previously precipitated sulfur minerals (Fig. 1) (Garrels and Lerman 1981; Berner 1987; Fike et al. 2015).

Combined, these sources have been calculated to provide c. 3.12 × 1012 moles of sulfur per year to the ocean, with an average isotopic composition of c. 3‰ (Canfield 2004; Kah et al. 2004). More recent studies suggest that c. 2.80 × 1012 moles of this comes from the weathering of all forms of terrestrial sulfur (Burke et al. 2018), with sulfide and sulfate weathering contributing 46 and 54% to the total riverine flux, respectively. However, there is much uncertainty regarding flux sizes as previous models have used sulfide to sulfate weathering ratios of 1 : 2 (Kump and Garrels 1986; Wu et al. 2010) or closer to 1 : 4 (Lenton et al. 2018; Mills et al. 2021) based on data from observations and/or mass balance calculations (assuming a modern day steady state). If 2.80 × 1012 moles of the sulfur input budget does come from weathering, then this would imply a more minor role for sulfur degassing at the present day (0.32 × 1012 mol a−1), but some models suggest a greater contribution from degassing to the total sulfur input (0.75 × 1012 mol a−1) (Berner 2006; Lenton et al. 2018) and this was likely exacerbated in the past when tectonic activity, and thus degassing rates, were sometimes higher (e.g. Brune et al. 2017; Merdith et al. 2019; Wong et al. 2019; Mills et al. 2021).

The major sinks of sulfate from the ocean are via the burial of sulfides and sulfates and, although both fluxes are comparable in size and sulfate burial imparts little to no isotopic fractionation from seawater, the formation and burial of sulfides can incur a significant isotopic fractionation, leading to an ocean that is isotopically heavier (c. 21‰) (Rees et al. 1978) than the sulfur sources (e.g. Wortmann et al. 2001; Canfield 2004; Kah et al. 2004; Brunner and Bernasconi 2005; Fike et al. 2015; Sim et al. 2023). Disproportionation processes may complicate the picture here (e.g. Canfield and Teske 1996; Bottrell and Newton 2006; Johnston 2011; Sim et al. 2011 and references cited therein); however, a recent study has suggested that disproportionation in the natural environment has less of an effect on sulfur isotope fractionation than has been inferred from earlier culture experiments (Tsang et al. 2023).

Dissolved marine sulfate ([SO42−]sw) concentrations are believed to have varied significantly over geological time (Kah et al. 2004, 2016; Turchyn and DePaolo 2019; Weldeghebriel et al. 2022), although the range over which they varied remains debatable (e.g. Algeo et al. 2015 and references cited therein). The interest in understanding these variations in [SO42−]sw through time, from an Earth system perspective, is multifaceted, but principally lies with its influence on atmospheric and marine oxygen levels (e.g. Krause et al. 2018; Shields et al. 2019; Shi et al. 2022), although there is a growing interest in how sulfur can affect atmospheric carbon dioxide and climate (Charlson et al. 1987; Shields and Mills 2021) – for example, the enhanced release of carbon dioxide as a result of carbonate weathering via sulfuric acid generation from pyrite oxidation has been debated (Torres et al. 2014; Maffre et al. 2021).

A number of methods have been used to reconstruct past [SO42−]sw, which we discuss here. These include:

  1. Determining sulfate concentrations from authigenic minerals such as halite and francolite (carbonate-bearing fluorapatite).

    •    a. Concentrations from fluid inclusions in marine halite can be back-calculated by taking into account processes such as: the degree of evaporation of the initial seawater; the precipitation of mineral phases (e.g. gypsum/anhydrite) during halite formation; and the possible alteration of the initial seawater via contributions from non-marine aqueous solutions or from seawater reactions with sediments on the way to the evaporitic basin, depending on the type of basin forming (e.g. lagoonal or salina) (Horita et al. 2002; Brennan et al. 2004; Lowenstein et al. 2005; Bąbel and Schreiber 2014; Warren 2021; Weldeghebriel et al. 2022).

    •    b. Sulfate can be structurally substituted into authigenic minerals, such as francolite, and precipitated in less restricted marine settings than halite (e.g. Shields et al. 2004; Hough et al. 2006). Concentrations of sulfate in francolite are typically determined by analysing powdered samples of bulk rock, either by X-ray fluorescence (Goldberg et al. 2011) or by electron probe microanalysis (Broom-Fendley et al. 2021).

  2. Applying ‘forwards’ biogeochemical modelling to predict [SO42−]sw (Arvidson et al. 2013; Krause et al. 2018; Lenton et al. 2018; Mills et al. 2021, 2023), where the sulfur fluxes are calculated based on external forcings (e.g. changes to degassing rates through time) and internal model parameters (e.g. the extent of ocean anoxia). More generally, the modelling of carbonate chemistry has been used to explain the relative scarcity of Archean and Paleoproterozoic gypsum deposits by way of exhaustion of [Ca2+]sw due to extensive evaporitic calcium carbonate formation when oceans were bicarbonate-rich (Grotzinger and Kasting 1993; Warren 2021), although substantial gypsum deposits occasionally formed (i.e. during the Lomagundi Event; Schröder et al. 2008; Blättler et al. 2018).

  3. Inferring concentrations, or relative changes in concentrations, from the modelling of sulfur isotopes (e.g. Habicht et al. 2002; Berner 2004; Kah et al. 2004; Wortmann and Paytan 2012; Crowe et al. 2014; Algeo et al. 2015; Zhu et al. 2021) and, very rarely, calcium isotopes (Blättler and Higgins 2014). For example, Algeo et al. (2015) calculated a maximum [SO42−]sw using the ‘rate method’ based on modelling of the isotopic fractionation observed between coeval sedimentary sulfate and sulfides (Δ34Ssulfate-sulfide) and the maximum rate of variation in CAS sulfur isotopes (which are used to infer the isotopic composition of seawater through time). Kah et al. (2004) and Gill et al. (2007) were the first to develop this kind of approach in using sulfur isotopes to reconstruct [SO42−]sw.

  4. An alternative method of reconstructing [SO42−]sw using sulfur isotopes was developed by Algeo et al. (2015), whereby they derived a relationship between the Δ34Ssulfate-sulfide associated with microbial sulfate reduction (MSR) and [SO42−] based on analyses of a wide range (in terms of salinity) of modern depositional environments. Hypersaline environments aside, they found a strong relationship between Δ34Ssulfate-sulfide and [SO42−] in aqueous settings and applied this MSR trend method to the record of Δ34Ssulfate-sulfide from CAS and pyrite samples to derive [SO42−]sw for the Phanerozoic.

  5. Using a numerical diffusion–advection–reaction model, originally used to calculate the δ34S of syndepositional pyrite formation in sediment pore waters (Lang et al. 2020), but modified such that the input of iron speciation, δ34Spy and δ34Ssulf (from CAS) data into the model allows for the back-calculation of [SO42−]sw (Zhu et al. 2021). This method has so far only been applied to the Cambrian and is subject to relatively high uncertainty, as it is assumed that the reaction rate constant and isotopic fractionation for dissimilatory sulfate reduction are equivalent to those observed in modern day marine sediment pore waters. However, for example, the isotopic fractionation for dissimilatory sulfate reduction can encompass a very broad range (Sim et al. 2011) and iron speciation data can only inform local redox conditions (Sperling et al. 2015; Poulton 2021). Similarly, iron speciation data combined with a one-dimensional water column reactive transport model has been used to estimate a low [SO42−]sw of <600 μM during Oceanic Anoxic Event 1a (Bauer et al. 2022).

We investigated changes to past [SO42−]sw using the biogeochemical box model GEOCARBSULFOR (Krause et al. 2018; Mills et al. 2023). The model is more traditionally used to compute variations in atmospheric oxygen and carbon dioxide levels, but can also estimate variations in the concentration of oceanic sulfate through Phanerozoic time. We compare our findings with other reconstructions from both geochemical data and other modelling studies and assess some of the potential impacts of variable oceanic sulfate on atmospheric oxygen.

Model background

GEOCARBSULFOR is a modified version of GEOCARBSULF (Berner 2006). It is a box model that uses a combined ocean–atmosphere box and four crustal reservoir boxes: organic and carbonate carbon and reduced (represented by pyrite) and oxidized (represented by gypsum) sulfur, with these boxes split into ‘young’ (i.e. freshly deposited sediments) and ‘ancient’ boxes (Fig. 2). Because it uses a combined ocean–atmosphere box, the model does not fully reflect the nuances of ocean carbonate chemistry and the impact this has on the operation of other elemental cycles and, consequently, model outputs such as climate predictions. As such, the model is useful for examining secular trends over geological time, but less so for individual Earth system perturbations on timescales <1 Myr. In addition, we do not include any impact that pyrite oxidation may have on the long-term climate because this is the subject of an ongoing debate (Torres et al. 2014; Maffre et al. 2021).

Although the first GEOCARBSULFOR model (Krause et al. 2018) used singular model runs, the model has recently been upgraded (Mills et al. 2023), such that it is now capable of running a Monte Carlo style ensemble simulation, with 5000 model runs conducted using parallel computing and encompassing the uncertainties of some key parameters, such as the δ13Ccarb record, which affect the operation of the carbon, sulfur and oxygen cycles, to produce an envelope of possible results.

The model is written in the MATLAB® language and uses an implicit variable-order, variable-step-size, ordinary differential equation solver (Shampine and Reichelt 1997), meaning that the model inputs and outputs are not constrained to the 10 Myr time intervals that the original GEOCARBSULF model (Berner 2006; Royer et al. 2014) was subject to. Although some of the geological data (e.g. the size of the total land area) that the model uses to help derive the sizes of the different carbon and sulfur fluxes remains coarse (varying over 10 Myr intervals), most of the data (e.g. the δ13C record, uplift rates and tectonic degassing) have been updated in recent years (Krause et al. 2018; Mills et al. 2023) to be of a higher resolution (0.1–1 Myr intervals).

Burial fluxes

As per the original GEOCARBSULF model (Berner 2006), the carbon cycle uses the geological δ13C record of marine carbonates (Cramer and Jarvis 2020) through time to estimate the organic carbon burial flux. However, instead of being determined by the geological δ34Ssw record, as in GEOCARBSULF, the sulfur cycle in GEOCARBSULFOR takes a forwards approach, where pyrite burial is dependent on the concentration of sulfate, the degree of anoxia in the oceans and the availability of organic carbon, with sulfate concentrations and the degree of anoxia being generated by internal model processes and organic carbon availability assumed to be related to the δ13C record. Krause et al. (2018) detail the necessity of this change to the sulfur cycle, but, briefly, the original GEOCARBSULF model used an oxygen-dependent equation to derive δ34S fractionation, whereby low oxygen concentrations resulted in unrealistically small δ34S fractionation values and consequently very large pyrite burial fluxes. In effect, this constrained minimum modelled atmospheric oxygen concentrations to be at approximately modern day levels for the entirety of the Phanerozoic, inconsistent with geochemical proxy data inferring widespread ocean anoxia in the early Paleozoic (Sperling et al. 2015).

As a result of this change, we can derive the isotopic signatures of the sedimentary sulfur reservoirs and all input and output fluxes to/from the ocean to use mass balance calculations to generate a synthetic δ34Ssw signal, which can be compared with proxy data from the geological record. A more elegant method of calculating pyrite burial and the δ34Ssw record through time may be to use reaction transport modelling (e.g. Wortmann and Chernyavsky 2007), but this would require separate atmosphere and ocean boxes in the model. In our baseline model runs, and similar to other models (Berner 2004; Lenton et al. 2018), gypsum burial is dependent upon the sulfate and calcium concentrations in the oceans, where sulfate abundance is calculated explicitly by the model and calcium concentrations are imposed from geological data (see later discussion).

Utility of the model and a major update

A key advantage of the GEOCARBSULFOR model, compared with its predecessors (Garrels and Lerman 1984; Berner 2001, 2006), is that it can make predictions of [SO42−]sw through time, whereas, for example, the oceanic reservoir in the GEOCARBSULF model (Berner 2006) was fixed at its present day value for the entire Phanerozoic because of isotope mass balance constraints. In addition, this change from a fixed to variable [SO42−]sw reservoir allows the decoupling of the δ13Ccarb and δ34Ssulf records and can take into account the different residence times of carbon and sulfur.

Here, we have made an update to the model and altered the normalized seawater calcium concentration ([Ca]sw) forcing, which affects the gypsum burial rate. Instead of using the record of Phanerozoic marine calcium concentrations ([Ca]sw) inferred from fluid inclusion data (Horita et al. 2002), GEOCARBSULFOR now uses an updated [Ca]sw record based on an expanded database of fluid inclusion data (Weldeghebriel et al. 2022). The model incorporates into the Monte Carlo ensemble the estimated range from fluid inclusion analyses (Fig. S1), with the assumption that the concentration product of [Ca]sw × [SO4]sw has the bounds of 150–450 mmol2 with an average of 319 mmol2 (see Weldeghebriel et al. 2022, their table S2). In order to be used by the model, the [Ca]sw data were first normalized to the present day value (11 mM) and a smoothing spline was then fitted to the data. To avoid confusion between the different model iterations, we rename the earlier versions as G18 (Krause et al. 2018) and G23 (Mills et al. 2023), whereas the version in this study remains named as GEOCARBSULFOR.

GEOCARBSULFOR predicts a general trend of increasing [SO42−]sw across the Phanerozoic (Fig. 3). The average of our modelled concentrations ranges between 9 and 16 mM from the terminal Ediacaran to the end of the Devonian and then increases to a peak of c. 25 mM at roughly the Guadalupian–Lopingian boundary (c. 259 Ma). The model estimates that concentrations then decrease such that, for the Jurassic and Cretaceous, [SO42−]sw is between 14 and 18 mM, before increasing to the present day value (28–29 mM).

Our GEOCARBSULFOR model results compare favourably with sulfate concentrations estimated from the MSR trend method (Algeo et al. 2015) for the whole of the Phanerozoic, although this method suggests variations on shorter timescales (Stebbins et al. 2019). Our results also compare reasonably with the fluid inclusion data, although they trend towards the higher end of the range (Weldeghebriel et al. 2022), but are often significantly higher than those predicted by CAS isotopes (see Fig. 3 legend for references). A reasonable fit to the sulfate data from fluid inclusions is expected because we use calcium concentrations from the same inclusions (Weldeghebriel et al. 2022) to partially drive gypsum burial in the model. However, other inputs and outputs to/from the ocean may still hold as great, or greater, control on [SO42−]sw, especially with the higher sulfur degassing rates in the past; we use degassing rates from the SCION model (Mills et al. 2021).

We note that while the calcium estimates (Weldeghebriel et al. 2022) reasonably assume that seawater total alkalinity has not significantly changed throughout the Phanerozoic (Turchyn and DePaolo 2019), some restricted bodies of water, such as lakes, have high alkalinity (Warren 2021), thus there is potential for inaccurate back-calculations of calcium concentrations if there is mixing of seawater with non-marine solutions on the way to the evaporitic basin.

There is a difference between the model and proxy reconstructions during the terminal Ediacaran (Fig. 3a). It could be that the model is underestimating the amount of gypsum burial occurring at this time (see Fig. 4b) and that [SO42−]sw should be closer to the minimum end of the uncertainty band, closer to levels previously suggested (Shi et al. 2022). Currently the model estimates [SO42−]sw to be between that predicted by the MSR trend and the fluid inclusion data used here. Indeed, when comparing previous modelling efforts (Fig. 3b), there is no consensus on [SO42−]sw for the late Ediacaran, with the COPSE and Berner models (Berner 2004; Lenton et al. 2018) agreeing with the MSR trend estimate, whereas the SCION (Mills et al. 2021) and G23 models suggest that [SO42−]sw was closer to 20 mM. Slightly later, during Cambrian Stage 2 (but possibly Stage 3, as there are some chronostratigraphic uncertainties), numerical modelling of the pyrite content and pyrite-sulfur isotopes of shale samples from the Yurtus Formation in the Tarim Block of northwestern China suggest a [SO42−]sw of 8.9–14 mM (Zhu et al. 2021), which broadly falls within our modelled range during this Epoch. However, CAS isotopes (He et al. 2019) suggest that [SO42−]sw levels during Cambrian Stages 2–4 (c. 524–512 Ma) were much lower (0.4–6.6 mM).

Generally, the new GEOCARBSULFOR model predicts a Phanerozoic [SO42−]sw record that exhibits the same overarching trend as other models: that of a long-term increase in oceanic sulfate with a ‘hump’ in the Carboniferous–Permian. The model predicts the same variation in [SO42−]sw as the MSR trend method (Fig. 3a), albeit the pre-Cenozoic peak occurs slightly before that predicted by Algeo et al. (2015), with the modelled average lying within the ±1SD MSR trend band or matching the fluid inclusion estimates (and sometimes both) for most of the last 540 Myr. A nuanced difference between the model results and the proxy data occurs during the Carboniferous, when GEOCARBSULFOR predicts higher [SO42−]sw than is inferred by either the MSR trend or fluid inclusion proxy reconstructions. There is some overlap between the GEOCARBSULFOR-calculated minimum and the fluid inclusion maximum value during the Carboniferous and, for the Early Mississippian (c. 359 Ma), the modelled uncertainty band overlies some of the range previously modelled by Gill et al. (2007) using CAS isotopes, albeit they use a version of the ‘rate method’ that can only provide estimates for the maximum [SO42−]sw level at that time.

Effects of evaporite deposition and weathering

In the current GEOCARBSULFOR model, the removal of marine sulfate into sediments is assumed to occur continuously and at a rate scaled to both the marine sulfate and calcium concentrations, while the weathering inputs through sulfate dissolution are largely dependent on the overall amount of sulfate in the continental crust. While large-scale formation may be episodic, mass reconstructions show evaporites to be nearly ubiquitous across the Phanerozoic (Warren 2021) and calcium sulfate in its various mineralogical phases can form in a number of deep marine settings via a variety of processes (see review in Van Driessche et al. 2019), thus continuous sulfate burial across multi-million-year timescales is a reasonable assumption.

Nevertheless, we build upon the current approach here in two ways. In Scenario 1, we use information regarding evaporite exposure and depositional area from palaeogeographical maps to use as scaled forcings for gypsum weathering and burial (Fig. 4a; Bluth and Kump 1991) and assess the effect they have on [SO42−]sw. In Scenario 2, we take published estimates of evaporite volumes deposited at various intervals throughout the Phanerozoic (Fig. S2) (Bąbel and Schreiber 2014), presume that c. 20% of the evaporites consisted of calcium sulfates, then convert this mass to a flux (in moles per million years; purple line, Fig. 4e) and add this additional flux to the background gypsum burial derived by the model from marine sulfate and calcium concentrations (resulting in the green line, Fig. 4e). Of course, with Scenario 1, this is data on all evaporites, not just gypsum/anhydrite and sulfate-bearing bitterns such as kainite, and halite is a significant evaporitic mineral (Hay et al. 2006; Warren 2021). However, by following other work (Berner 2004) and normalizing the evaporite areas to their present day values, we can use total evaporites as a proxy for sulfate evaporites through the Phanerozoic without making an assumption about the absolute amount of sulfate present in evaporites, only that this has not changed significantly over time.

Under Scenario 1, [SO42−]sw is lower than the baseline model for most of the Phanerozoic (Fig. 4b), with sulfate not exceeding 10 mM until the mid-Eocene, apart from a significant ‘hump’ during the Triassic–Jurassic. The modelled [SO42−]sw generally plots towards the low end of the range suggested by the MSR trend proxy, but is often a reasonable fit to the sulfate levels proposed by CAS isotopes. However, in the early Triassic, our modelled [SO42−]sw is either high (baseline) or begins to increase rapidly (Scenario 1), in stark contrast with CAS isotope estimations of a few mM (Luo et al. 2010; Song et al. 2014). Furthermore, at the Triassic–Jurassic boundary, [SO42−]sw is higher than alluded to by any of the three proxies. It is likely that, because some geological forcings used in the model are represented on (multi)-million-year timescales, transient perturbations to the sulfur cycle, such as at the Triassic–Jurassic boundary, are missed by the model and warrant further investigation with a different model that incorporates, for example, a multi-box ocean, carbonate chemistry and an iron cycle.

In addition, if we compare the gypsum burial forcing of Scenario 1 (Fig. 4a) with the gypsum burial flux of Scenario 2 (Fig. 4e), it is possible that there was more gypsum burial taking place at this time than can be inferred from the palaeogeographical maps of Bluth and Kump (1991). However, Scenario 1 does show that, if gypsum burial in the model is not solely constrained by the prescribed [Ca]sw and model-calculated [SO42−]sw, then low [SO42−]sw, as might have occurred during short-term (c. <2 Myr) events in the Phanerozoic (Wortmann and Chernyavsky 2007; Algeo et al. 2015), is achievable and thus the simple rate dependence on [SO42−]sw in our baseline version of the model may need further consideration in future work.

The introduction of scalings to the gypsum fluxes generally lowers the model computed δ34Ssw record (Fig. 4c) by a few per mil. Both Scenario 1 and the baseline model capture the long-term trend of the δ34Ssw record (Crockford et al. 2019), but miss some of the variability seen in the data, particularly in the Paleozoic, but also during the Cenozoic. This may be partly due to changes in sea-level and, consequently, the flooded continental area (see the compilation of van der Meer et al. 2017). This would influence the extent of shallow shelf settings where local processes, such as rapid sedimentation, increased the delivery of iron and organic carbon. The resuspension of marine sediments can dominantly affect the sulfur cycle and, consequently, the local δ34Spy (and thus presumably δ34Ssulf) recorded, particularly from the late Paleozoic to the present day (e.g. Aller et al. 2010; Leavitt et al. 2013; Pasquier et al. 2017, 2021; Lang et al. 2020). Such changes in sea-level and shelf area may occur over shorter timeframes than accounted for in the data used by the model – for example, the total land area available for subaerial weathering is represented as 10 Myr time slices and changes to the ocean topography and the effect on sediment burial processes is poorly represented or absent (Berner 2006). Again, this can be explored in the future with a more complex model.

Alternatively, our treatment of sulfur isotope fractionation in the model could potentially explain some of the differences between the model-generated δ34Ssw record and the proxy data. We investigated this by running the baseline version of the model several times, with the only change being the method by which sulfur isotope fractionation is calculated (see Fig. S3a and SI for further details). No singular method for deriving sulfur isotope fractionation can generate a modelled δ34Ssw record that is a good fit to the proxy data for the entire Phanerozoic (Fig. S3b), but some methods (e.g. using combined δ34Spy and δ34Ssulf records; Wu et al. 2010) reproduce δ34Ssw for certain timeframes (e.g. the Cenozoic) very well. None of the methods results in as low a δ34Ssw record as evidenced by the proxy data during the Pennsylvanian through to the end-Permian. However, under Scenario 1, where the minimum and maximum sulfur isotope fractionation from the last 570 Myr (Wu et al. 2010) is incorporated as a variable in the Monte Carlo ensemble, there is a closer match between the model and the data, indicating that changes to gypsum deposition, which are evidenced in the geological record (e.g. Bluth and Kump 1991; Šušnjara et al. 1992; Andeskie and Benison 2021; Johnson 2021), may be a key control on the δ34Ssw record (Han et al. 2023).

If [SO42−]sw was much lower in the Paleozoic, as suggested by Scenario 1 and the CAS isotopes (Gill et al. 2007, 2011; He et al. 2019), then the residence time of sulfur in the oceans might have been shorter than the 12.5 Myr at present, meaning that both the [SO42−]sw and δ34Ssw records were susceptible to short-term perturbations in sulfur weathering and/or burial. Indeed, our modelling suggests that the residence time of [SO42−]sw was less than 6 Myr for most, if not all, of the early to mid-Paleozoic (end-Devonian) and only briefly (across the Permo-Triassic) approached the 12.5 Myr of the present day (Fig. S4). The shorter residence time may explain why different proxy methods produce incongruent [SO42−]sw values in the Cambrian (He et al. 2019; Zhu et al. 2021), allowing either rapid fluctuations in the oceanic sulfate reservoir or spatial heterogeneity. With regards to atmospheric oxygen (Fig. 4d), our changes to the gypsum side of the sulfur cycle serve to depress oxygen levels by a few per cent throughout because pyrite burial is partially dependent on the amount of sulfate in the ocean. However, the predicted high levels of [SO42−]sw in the Triassic–Jurassic extend the Phanerozoic atmospheric oxygen peak to almost three times its original timeframe as pyrite burial is greatly enhanced, as evidenced by the increase in modelled δ34Ssw (Fig. 4c).

In Scenario 2, increases in gypsum burial are more episodic in nature, thus [SO42−]sw is more similar to our baseline results (Fig. 4f) than in Scenario 1. Although the baseline results produce [SO42−]sw levels higher than the MSR trend and/or fluid inclusion proxies during the Carboniferous, Scenario 2 is higher during the Devonian as well as the Carboniferous. Comparing the model outputs with the proxy data for the δ34Ssw record (Fig. 4g), the model prediction is a few per mil too low during the Devonian, but too high during the Carboniferous–Permian, implying that for the baseline and Scenario 2 (there is a negligible difference between the two outputs), the model could be better resolved with the data by burying more pyrite in the Devonian and more gypsum, as discussed previously, in the Carboniferous–Permian. There is little difference in the atmospheric oxygen predictions (Fig. 4h), except in the Cenozoic, where Scenario 2 generates slightly higher oxygen as a result of enhanced pyrite burial due to greater [SO42−]sw levels.

Comparing the overall trend of modelled [SO42−]sw with atmospheric oxygen during the Phanerozoic provides some interesting results. [SO42−]sw and oxygen increase in tandem across the Paleozoic, resulting in a peak in the Permian, and then decrease during the Mesozoic. However, at the end of the Cretaceous for our baseline model run (Fig. 3) and the first attempt at modifying the gypsum fluxes (Fig. 4b), this relationship decouples and [SO42−]sw begins to increase again while oxygen generally decreases. A possible explanation for this is that various reconstructions of long-term uplift rates, using either strontium isotopes or back-calculating from terrigenous sediment deposits, indicate that uplift and consequently erosion rates began to markedly increase during the Cenozoic (Fig. S5a) (Berner and Kothavala 2001; Hay et al. 2006; Mills et al. 2019). Previous work (e.g. Calmels et al. 2007) has shown that there is a strong relationship between pyrite oxidation and erosion rates. As pyrite oxidation is a source for sulfur in the oceans, but is a sink for atmospheric oxygen, it seems possible that this process exerted a considerable control on the Earth system during the Cenozoic.

At the same time, organic carbon availability, which affects not only organic carbon burial but also pyrite burial (via sulfate reduction; see Fig. 2), likely decreased across the Cenozoic according to our model (based on lower predicted rates of organic carbon burial; Fig. S5b), while gypsum burial is also moderate, apart from a potential pulse in the Neogene (Fig. 4a, e) and gypsum weathering is sustained (Fig. 4a). Thus, more sulfate is entering the oceans, but less is leaving, while more oxygen is being consumed by oxidative processes and the sources of oxygen are shrinking. However, increased pyrite oxidation in tandem with decreased pyrite burial would result in a decrease in δ34Ssw values across the Cenozoic and, although a large proportion of samples exhibit very low δ34Ssw during this Era (Crockford et al. 2019), tabulated data instead indicate an increase of c. 3‰ across the early to mid-Eocene, with relatively stable values thereafter (Wu et al. 2010).

The stark rise in sulfate concentrations across the Cenozoic suggested by both our model and fluid inclusion data (Fig. 3a) implies not only an increase in pyrite weathering, but also sulfate weathering, especially considering the present day uncertainty about the contributions of pyrite v. gypsum weathering to the total riverine sulfur flux. A large increase in sulfate weathering at this time may have been due to the widescale dissolution of Late Neoproterozoic to Early Cambrian evaporites during the initial collision of India and Eurasia under a hot-house climate (Wortmann and Paytan 2012; Rae et al. 2021; Shields and Mills 2021). Thus it is possible that any variance in δ34Ssw due to changes in the reduced sulfur fluxes was masked by an increase in oxidized sulfur fluxes with very positive δ34S.

The two different updates to the gypsum fluxes implemented in this study have a negligible influence on atmospheric carbon dioxide concentrations, likely because the model does not include carbonate chemistry and calcium carbonate saturation, unlike other models (Shields and Mills 2021). In GEOCARBSULFOR, the effects from changing the sulfur fluxes are largely self-contained to the sulfur cycle, thus the variations in atmospheric oxygen (Fig. 4d, h) are mainly due to alterations in the amount of pyrite being buried.

The calcium cycle

Interestingly, from biogeochemical box model perspectives, although the COPSE (Lenton et al. 2018), MAGic (Arvidson et al. 2013), Berner's Ca–Mg–SO4 model (2004) and the three iterations of the GEOCARBSULFOR model (Krause et al. 2018; Mills et al. 2023; this study) all predict an increase in [SO42−]sw starting at roughly the Devonian–Carboniferous boundary, leading to a peak in concentrations during the Permo-Triassic, the SCION model (Mills et al. 2021) suggests a minor decrease over this interval (Fig. 3b). The difference in results from the SCION model (for the Phanerozoic in general, but especially during this period of time) may be due to the fact that, unlike the other models, gypsum burial in the SCION model is not partially constrained by the [Ca]sw record. Instead, gypsum burial in the SCION model follows other models, such as COPSE, in being dependent on [SO42−]sw, but is also reliant on a newly introduced normalized forcing: the extent of palaeo-shorelines, which acts as a proxy for basin restriction (Mills et al. 2021). The COPSE model is also unable to reproduce the general trend of Phanerozoic sulfate levels when its calcium forcing is removed (see Lenton et al. 2018, their fig. 7). We therefore investigated the influence of [Ca]sw on various model outcomes by running the baseline version of GEOCARBSULFOR in the Monte Carlo mode once more, but this time excluding normalized [Ca]sw as a forcing for gypsum burial rates.

We found that by excluding calcium from GEOCARBSULFOR, the model now predicts a Phanerozoic [SO42−]sw record that is very similar to that generated by the SCION model (Fig. 5a), not only in terms of absolute concentrations, but also with regard to some of the temporal variance, including the minor decrease across the Carboniferous to Triassic. [Ca]sw levels, as inferred by both fluid inclusions (Weldeghebriel et al. 2022) and modelling using calcium isotopes (Farkaš et al. 2007), were very high in the Ordovician and, although they experienced a long-term decrease across the rest of the Paleozoic, this was still to relatively high levels (about 1.5 times the present day levels). With such elevated [Ca]sw, gypsum burial rates would also likely be promoted, which can be seen in the areal extent of evaporite deposition in the geological record (Bluth and Kump 1991), thus [SO42−]sw would be expected to be low. Removing the dependence of the normalized [Ca]sw record in the model reduces sulfate formation as an exit channel and, because pyrite formation is constrained by both organic carbon availability and ocean oxygenation, as well as [SO42−]sw, pyrite burial is limited in its ability to compensate for a decrease in gypsum burial rates and therefore [SO42−]sw builds up during the Paleozoic in this model run. High [Ca]sw in the Cretaceous operates in the same manner, leading to [SO42−]sw that plots at the high end, or above, that predicted by the MSR trend.

The inclusion of a dynamically calculated calcium cycle in the model may help to improve [SO42−]sw predictions while also freeing up the calcium fluid inclusion data to be used as a further proxy for bench testing the model. The original version of COPSE (Bergman et al. 2004) included a partial calcium cycle, but the modelled [Ca]sw estimates missed a peak in values in the Mesozoic. Other models with the inclusion of a dynamic calcium cycle have successfully captured the Mesozoic hump in [Ca]sw, but to the somewhat detriment of present day estimates of [Mg]sw (Berner 2004; Farkaš et al. 2007) and/or [SO42−]sw (Berner 2004; Arvidson et al. 2013). The modelling work of Hansen and Wallmann (2003) does finish, uniquely, with the correct values for all [Mg]sw, [Ca]sw and [SO42−]sw reservoirs for the present day (as estimated at the time), but the model estimates only cover the last 150 Myr. Suffice to say, the evidence suggests that models without a consideration of calcium concentrations are not capturing the observed trend of [SO42−]sw for at least the Paleozoic. [Ca]sw across the Phanerozoic appears to exhibit a significant control on sulfate concentrations, but modelling this is not without its issues.

In terms of its effect on other modelled outputs, excluding normalized [Ca]sw when calculating gypsum burial rates results in atmospheric oxygen that is generally higher by a few per cent (Fig. 5c) due to a small shift in the fraction of sulfur leaving the oceans in the form of pyrite, as can be seen by the generally elevated marine δ34S record (Fig. 5b) generated by the model. In GEOCARBSULFOR, calcium carbonate burial is not dependent on the [Ca]sw record, therefore there is a negligible effect on the atmospheric carbon dioxide record (Fig. 5d), but other work has shown that the sulfur cycle may have a significant role in modulating the seawater carbonate system and climate via evaporite deposition/weathering (Shields and Mills 2021). The inclusion of normalized [Ca]sw as a forcing for calcium carbonate burial may have a significant effect on all modelled outcomes and should therefore be investigated further in future studies.

To investigate how marine sulfate concentrations have changed across the Phanerozoic, we ran the biogeochemical model GEOCARBSULFOR, finding that [SO42−]sw steadily increases through the Paleozoic, reaching a near-modern-day level peak in the Permian, decreases significantly in the Mesozoic and then increases again to the present day level of 28–29 mM, in good agreement with geochemical proxies and other modelling studies. The increase in [SO42−]sw is mirrored by an increase in atmospheric oxygen levels until the end of the Cretaceous, where this relationship starts to break down and atmospheric oxygen decreases to the present day level of c. 21%.

The model predicts a Carboniferous [SO42−]sw that is a little higher than suggested by proxies. We updated the model's gypsum fluxes by incorporating geological evidence from the evaporite record and analysed the effect this had on [SO42−]sw and other modelled outputs. We show that incorporating the observed high rates of gypsum burial during the Carboniferous could cause the model [SO42−]sw record to better match the MSR trend proxy and, similarly during the Cretaceous, this could cause the modelled [SO42−]sw values to be closer to those suggested by CAS isotopes. However, when forcing the model with whole-Phanerozoic records of gypsum burial, the results outside of these timeframes are generally less consistent with the geological record. This could be partially due to a variable abundance of gypsum within individual evaporites.

We also investigated the impact that the calcium concentration of the oceans may have on the Phanerozoic [SO42−]sw record by excluding it as a forcing in the calculation for gypsum burial. Based on the model results, it seems that calcium has a major role in controlling [SO42−]sw over time, with high calcium levels leading to much lower modelled sulfate concentrations during the Paleozoic, in agreement with various lines of geological evidence, consequently affecting atmospheric oxygen predictions. Overall, our study highlights the considerable influence of the sulfur cycle on other elemental cycles, and vice versa, and we conclude that understanding and modelling the coupled calcium and sulfur cycles (e.g. without prescribing the calcium concentration) remains a significant challenge that has not yet been achieved for the entire Phanerozoic.

We would like to thank both the editor and the referees for their constructive comments and advice which helped to improve this work.

AJK: conceptualization (lead), data curation (lead), investigation (lead), methodology (lead), project administration (lead), writing – original draft (lead), writing – review & editing (lead); GAS: data curation (supporting), project administration (supporting), writing – original draft (supporting), writing – review & editing (supporting); RJN: data curation (supporting), project administration (supporting), writing – original draft (supporting), writing – review & editing (supporting); BJWM: conceptualization (supporting), data curation (supporting), project administration (supporting), writing – original draft (supporting), writing – review & editing (supporting).

We gratefully acknowledge funding support from a Leverhulme Research Fellowship (UK) (RF-2019–435) and Natural Environment Research Council (NERC) grants NE/P013643/1 (BETR programme) and NE/R010129/1 to GAS and BJWM and a NERC large grant NE/N018559/1 to RJN.

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

The model code and the various datasets required to run the model are freely available to download from https://doi.org/10.5281/zenodo.10814777 and all published versions of the GEOCARBSULFOR model (including this one) are available to download from https://github.com/Alexjkrause/GEOCARBSULFOR

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/)