Gold Open Access: This paper is published under the terms of the CC-BY license.


Early Cretaceous (145–100 Ma) rocks record a ∼5‰ negative shift in the sulfur isotope composition of marine sulfate, the largest shift observed over the past 130 m.y. Two hypotheses have been proposed to explain this shift: (1) massive evaporite deposition associated with rifting during opening of the South Atlantic and (2) increased inputs of volcanically derived sulfur due to eruption of large igneous provinces. Each process produces a very different impact on marine sulfate concentrations, which in turn affects several biogeochemical phenomena that regulate the global carbon cycle and climate. Here we present sulfur isotope data from Resolution Guyot, Mid-Pacific Mountains (North Pacific Ocean), that track sympathetically with strontium isotope records through the ∼5‰ negative sulfur isotope shift. We employ a linked sulfur-strontium isotope mass-balance model to identify the mechanisms driving the sulfur isotope evolution of the Cretaceous ocean. The model only reproduces the coupled negative sulfur and strontium isotope shifts when both hydrothermal and weathering fluxes increase. Our results indicate that marine sulfate concentrations increased significantly during the negative sulfur isotope shift and that enhanced hydrothermal and weathering input fluxes to the ocean played a dominant role in regulating the marine sulfur cycle and CO2 exchange in the atmosphere-ocean system during this interval of rapid biogeochemical change.


The sulfate content of the ocean plays an important role in regulating the global carbon cycle and the chemical composition of the ocean-atmosphere system. For example, in the present sulfate-rich ocean, microbial sulfate reduction (MSR) remineralizes as much as 50% of the organic matter in coastal marine sediments (Jørgensen, 1982) and affects phosphorus recycling efficiency and therefore marine primary productivity (Ingall et al., 1993; Van Cappellen and Ingall, 1996; Adams et al., 2010). Moreover, MSR increases carbonate alkalinity via bicarbonate production, which in turn affects the partitioning of CO2 between the atmosphere and ocean (Zeebe and Wolf-Gladrow, 2001). However, marine sulfate concentrations have varied considerably over geologic time (Lowenstein et al., 2001), affecting the marine carbon cycle, the evolution of Earth’s climate system, and the long-term redox balance of the ocean-atmosphere system.

The relative importance of marine sulfur inputs and outputs through time can be tracked in part by reconstructing the S isotope composition of sulfate phases (i.e., barite, calcium sulfate, and carbonate-associated sulfate) preserved in marine sedimentary rocks. The S isotope composition of marine sulfate (δ34Ssulfate) represents a balance between the input of 34S-depleted S from riverine (continental weathering) and hydrothermal sources and the removal of 34S-depleted S via MSR and associated pyrite burial (Canfield, 2001).

An ∼5‰ negative excursion extending from ca.126 to 100 Ma represents the most distinguishing feature of the δ34Ssulfate record spanning the past 130 m.y. (Paytan et al., 2004; Fig. 1). While original interpretations of this event implicated a combination of increased hydrothermal and weathering inputs and decreased pyrite burial rates as possible triggers (Paytan et al., 2004), more recent work postulates that increased evaporite deposition (calcium sulfate) associated with the opening of the South Atlantic during this time dramatically lowered marine sulfate concentrations and facilitated a drop in relative pyrite burial rates (Wortmann and Chernyavsky, 2007; Wortmann and Paytan, 2012). Either of these processes operating in isolation could account for the negative δ34Ssulfate excursion. However, each has a very different impact on marine sulfate concentrations: increased volcanism and continental weathering rates increase sulfate delivery to the ocean, whereas evaporite deposition removes it during calcium sulfate precipitation. Although there is geologic evidence for both massive evaporite deposition (Hay et al., 2006; Davison, 2007; Chaboureau et al., 2013) and extensive volcanism due to the formation of several large igneous provinces (LIPs) during the Early Cretaceous (Chandler et al., 2012; Mills et al., 2014), the relative importance of these factors in regulating marine sulfate levels and attendant biogeochemical change remains unconstrained.

The Cretaceous Period (ca. 145–65 Ma) is also punctuated by several ocean anoxic events (OAEs), short-term (<1 m.y.) intervals of carbon cycle disruption during which massive amounts of organic carbon were buried in marine sediments (Schlanger and Jenkyns, 1976). OAE1a occurred at ca. 125 Ma (Ogg et al., 2012), although its temporal relationship with the negative δ34Ssulfate shift has not been documented unambiguously (Gomes et al., 2016). Given the strong coupling between the geochemical cycles of carbon and sulfur, improved understanding of how ocean chemistry and sulfate concentrations evolved over this period could provide insight into the factors that conditioned the oceans for carbon cycle instability (Wortmann and Chernyavsky, 2007; Adams et al., 2010; Gomes et al., 2016).

One way to decipher the importance of increased volcanism and/or weathering rates on the sulfur cycle during this time is through the use of strontium (Sr) isotopes, because the geochemical cycles of S and Sr are linked through shared input fluxes of riverine (weathering) and hydrothermal inputs. The 87Sr/86Sr ratio of seawater represents a balance between hydrothermal (relatively unradiogenic, lower 87Sr/86Sr) and terrestrial weathering (relatively radiogenic, higher 87Sr/86Sr) inputs (Palmer and Edmond, 1989). While the relative ratio of these input fluxes largely controls the isotopic composition of the marine Sr reservoir, the output flux of pyrite burial exerts strong control over δ34Ssulfate due to the large kinetic isotope effect associated with MSR and accompanying pyrite formation (Canfield and Teske, 1996). Thus, changes in the input fluxes (either in magnitude or isotopic composition) should be recorded in both isotope systems, while changes in the output flux (e.g., pyrite burial) should only affect the S cycle. If the Early Cretaceous negative S isotope shift was triggered by massive evaporite deposition and an attendant decrease in global pyrite burial rates, S and Sr isotope records should be decoupled. By contrast, if massive volcanism were responsible for driving the S isotope shift, one would expect a sympathetic shift in the 87Sr/86Sr ratio of seawater (as recorded in marine carbonate).

We present new S isotope data from Resolution Guyot, Mid-Pacific Mountains, that track sympathetically with Sr isotope records through the ∼5‰ negative sulfur isotope shift, and employ a linked sulfur-strontium isotope mass-balance model to identify the mechanisms driving the sulfur isotope evolution of the Cretaceous ocean.


Sulfur and Carbon Isotope Analysis

Carbonate-associated sulfate (CAS) was extracted by standard techniques (Hurtgen et al., 2006) and S isotope results are reported as per mil (‰) deviations from Vienna Canyon Diablo troilite (VCDT), using the conventional (δ34S) notation. Sulfur isotope results were reproducible within ±0.2‰, based on repeat analyses of standards. Carbon isotope results are reported as deviations (‰) from the carbon isotope composition of the Vienna Peedee belemnite (VPDB) standard and were reproducible within ±0.1‰ based on repeat analyses of standards. (See the GSA Data Repository1 information for additional methods information.)

Coupled S-Sr Cycle Box Model

A box model was constructed to track the mass and isotopic composition of marine sulfate and Sr using the following equations: 

where MS is the mass of sulfate in the ocean; MSr is the mass of Sr in the ocean; FHSr, FWSr, and FDSr are the hydrothermal, weathering, and diagenetic input fluxes of Sr, respectively; h is the S/Sr ratio of the hydrothermal flux; w is the S/Sr ratio of the weathering flux; FpyS and FSevap are the burial fluxes of pyrite and evaporite, respectively; and τSr (constant) is the residence time of Sr in the ocean. The S and Sr isotope reservoirs are described by coupled isotope mass-balance equations: 
where RSSW is the S isotope composition of seawater sulfate; RSrSW is the 87Sr/86Sr ratio of seawater Sr; RHSr and RHS are the Sr and S isotope composition of hydrothermal inputs, respectively; RSrW and RSW are the Sr and S isotope composition of weathering inputs, respectively; RSrD is the Sr isotope composition of the diagenetic input; Δpy is the average isotope fractionation factor associated with pyrite deposition; and Fpy is the burial flux of pyrite. The initial steady state for the Sr cycle was determined assuming that the magnitude of the hydrothermal flux and the 87Sr/86Sr ratios of the input fluxes (both hydrothermal and weathering) during the Early Cretaceous were comparable to the modern. The weathering flux was then calculated to achieve the observed pre-excursion RSrSW (see the Data Repository for table detailing initial steady-state conditions). Steady state for the S cycle was determined by modifying the scaling parameters for the hydrothermal and weathering fluxes (h and w) using estimates for the modern S/Sr ratios of hydrothermal and weathering inputs (Palmer and Edmond, 1989; Arthur, 2000; Halevy et al., 2012). Furthermore, the model assumes an initial sulfate concentration of 6 mM (as opposed to modern levels of 28 mM), in accordance with Cretaceous estimates based on the chemical composition of fluid inclusions encased in halite (Lowenstein et al., 2001). The model also assumes that the relative ratios of S/Sr in hydrothermal and terrestrial weathering inputs have remained relatively constant over the past 130 m.y. Additional complexity was added to the model to explore the impact of a sulfate concentration-dependent pyrite burial term (see the Data Repository for discussion).


We extracted CAS from drill core collected at Resolution Guyot (Ocean Drilling Program [ODP] Site 866) to examine the relationship between δ34Ssulfate and δ13Ccarbonate and previously published Sr isotope data (Jenkyns et al., 1995) from Early Cretaceous rocks. These data, plotted versus time (using the updated geological time scale of Ogg et al., 2012) in Figure 2 (see the Data Repository for age assignment details), indicate that while δ34Ssulfate values are variable between ca. 130 and 126 Ma, there is a clear ∼3‰ negative shift from an average of ∼19‰ that stratigraphically correlates with the OAE1a positive δ13C excursion. δ34Ssulfate continues to fall after OAE1a for a total S isotope shift of ∼4.5‰ before rebounding to pre-excursion values. The magnitude and timing of the δ34Ssulfate negative shift are consistent with previously published δ34Sbarite data (Paytan et al., 2004). Importantly, the Sr isotope data (Jenkyns et al., 1995) shift sympathetically with δ34Ssulfate and suggest that the geochemical cycles of S and Sr had roughly similar oceanic residence times and were at least partly coupled through this time period.


The Early Cretaceous record contains evidence for both extensive LIP volcanism and massive evaporite deposition; however, the relative timing of these events is critically important for the interpretation of the S and Sr isotope trends. A recent synthesis of the tectonic events associated with the breakup of South Africa and South America suggests that massive salt deposition in the South Atlantic basin occurred in the late Aptian (ca. 116–113 Ma; Chaboureau et al., 2013), whereas the major LIP events are thought to have occurred earlier. Mills et al. (2014) assigned the emplacement of the Ontong-Java and Manihiki Plateaus to the early Aptian (125 Ma) and the initiation of Kerguelen Plateau activity to 118 Ma.

Within this temporal context, we generated a coupled S-Sr box model to examine whether increases in hydrothermal activity, followed by massive evaporite deposition, are capable of reproducing the sympathetic shifts in δ34Ssulfate and 87Sr/86Sr captured at Resolution Guyot. In the model we were able to reproduce the ∼4.5‰ negative δ34Ssulfate shift by simulating two periods of increased hydrothermal activity (scenario 1, Fig. 3): a large, short time-scale increase in volcanic activity (8× increase for 0.5 m.y.) associated with the emplacement of the Ontong Java Plateau coincident with OAE1a, followed by a period of slightly elevated hydrothermal activity (1.7× hydrothermal flux) for 5.5 m.y.; and a second, smaller pulse of volcanism (2.3× hydrothermal flux for 5 m.y.), associated with the emplacement of the Kerguelen Plateau. The hydrothermal flux is maintained at 1.1× the original flux for the remainder of the simulation. While this elevated hydrothermal flux generates a δ34Ssulfate shift similar in magnitude to that observed in the Resolution Guyot record, it produces unreasonably low 87Sr/86Sr ratios (blue line, Fig. 3A).

We were only able to reproduce the negative S and Sr isotope shifts recorded in Resolution Guyot rocks by increasing both hydrothermal and weathering fluxes (scenario 2, Fig. 3; see caption for flux increases). An important result of this model is that the negative S isotope shift is accompanied by a doubling of the marine sulfate reservoir (red line, Fig. 3C). If a period of increased evaporite deposition is added to the model corresponding to the late Aptian South Atlantic evaporite deposits (scenario 3, 6× increase in evaporite burial from 116 to 113 Ma; green line, Fig. 3), then the associated drawdown in marine sulfate concentration allows the δ34Ssulfate reservoir to rebound on a time scale similar to, or slightly faster than, the Sr isotope reservoir, as observed in the Resolution Guyot record (Fig. 2). This result only holds if pyrite burial rates are assumed to be independent of marine sulfate concentrations.

Although the timing of the δ34Ssulfate negative shift recorded in CAS at Resolution Guyot is consistent with previously published δ34Sbarite data (Paytan et al., 2004), the return of δ34Ssulfate to pre-excursion values occurs ∼10 m.y. prior to the δ34Sbarite record. This discrepancy is not likely the result of temporal reconstruction errors, but rather indicates that this portion of the Resolution Guyot δ34SCAS record was influenced by either local environmental conditions and/or postburial alteration (see the Data Repository for more discussion). Given that the δ34Sbarite record was generated using multiple sections (Fig. 1), we assume it is more representative of global δ34Ssulfate evolution during this time. Wortmann and Chernyavsky (2007) proposed that both the negative δ34Ssulfate shift and the prolonged return to pre-excursion values resulted from evaporite burial; in their model, they prescribed a relationship between sulfate concentration and pyrite burial, where decreased sulfate concentrations (<10 mM) lead to decreased global pyrite burial rates. However, recent work suggests that massive evaporite deposition occurred in the South Atlantic in the late Aptian, ∼10 m.y. after the initiation of the negative S and Sr isotope shifts. Therefore, we propose that a simultaneous increase in the hydrothermal and weathering input fluxes associated with an increase in LIP volcanism best explains the coupled negative shift in the S and Sr isotope records. If we include massive evaporite deposition in our model during the late Aptian (9 m.y. after the initial onset of LIP volcanism) and assign a sulfate concentration-dependent pyrite burial flux (Wortmann and Chernyavsky, 2007; scenario 4, Fig. 3), the time scale of the δ34Ssulfate rebound more closely matches the δ34Sbarite curve (Fig. 1; Paytan et al., 2004) and allows for marine sulfate concentrations to remain low for the interval leading up to OAE2 ca. 95 Ma (Wortmann and Chernyavsky, 2007; Adams et al., 2010; Owens et al., 2013; Gomes et al., 2016).


Our findings, in conjunction with recent temporal constraints for LIP emplacement and evaporite deposition during the Early Cretaceous, demonstrate that coupled increases in hydrothermal and weathering inputs, followed by a period of massive evaporite deposition, best explain the negative S and Sr isotope shifts recorded in rocks drilled at Resolution Guyot. This indicates that marine sulfate concentrations likely increased through most of the Aptian interval (including OAE1a) before dropping to lower levels via late Aptian evaporite deposition. Furthermore, the simultaneous increases in volcanism and weathering required to reconcile the S and Sr isotope records highlight the importance of feedback mechanisms in regulating atmospheric CO2 levels (Berner et al., 1983; Kump et al., 2000). It has long been recognized that enhanced volcanism could increase atmospheric CO2 concentrations and elevate weathering rates. However, in isolation, the marine Sr isotope record cannot be used to infer simultaneous increases (or decreases) in both hydrothermal and riverine fluxes. We propose the use of coupled S and Sr isotope measurements as a new tool for interpreting isotope records and identifying the processes that affect marine sulfate levels, and therefore the C cycle and the evolution of Earth’s climate system, particularly during large-scale plate tectonic reorganization.

This work was supported by National Science Foundation grant EAR-0955969 to Hurtgen. This manuscript benefited from particularly helpful reviews by D. Fike and U. Wortmann.

1 GSA Data Repository item 2017153, detailed method and site description (including age assignments), S-Sr box model development and discussion of fidelity of the sulfur isotope record, is available online at, or on request from