Mitigation of Extreme Ocean Anoxic Event Conditions by Organic Matter Sulfurization

Past occurrences of widespread and severe anoxia in the ocean have frequently been associated with abundant geological evidence for free hydrogen sulfide (H2S) in the water column, so‐called euxinic conditions. Free H2S may react with, and modify, the chemical structure of organic matter settling through the water column and in marine sediments, with hypothesized implications for carbon sequestration. Here, taking the example of Ocean Anoxic Event 2, we explore the potential impact of organic matter sulfurization on marine carbon and oxygen cycling by means of Earth system modeling. Our model experiments demonstrate that rapid sulfurization (ksulf≥ = 105 M−1 year−1) of organic matter in the water column can drive a more than 30% enhancement of organic carbon preservation and burial in marine sediments and hence help accelerate climate cooling and Ocean Anoxic Event 2 recovery. As a consequence of organic matter sulfurization, we also find that H2S can be rapidly scavenged and the euxinic ocean volume reduced by up to 80%—helping reoxygenate the ocean as well as reducing toxic H2S emissions to the atmosphere, with potential implications for the kill mechanism at the end‐Permian. Finally, we find that the addition of organic matter sulfurization induces a series of additional feedbacks, including further atmospheric CO2 drawdown and ocean reoxygenation by the creation of a previously unrecognized net source of alkalinity to the ocean as H2S is scavenged and buried.


Introduction
Oceanic Anoxic Events (OAEs) are severe perturbations of global biogeochemical cycling associated with widespread deposition of laminated organic matter-rich sediments under anoxic/euxinic conditions-so-called black shales (Schlanger & Jenkyns, 1976). Originally envisaged as encompassing a series of three major and five minor events that occurred during the Jurassic and Cretaceous periods (Jenkyns, 2010), evidence for the transient occurrence of extensive ocean anoxia/euxinia has also been found associated with the Late Ordovician, Late Devonian, end-Permian, and Early Jurassic extinction events . Hydrogen sulfide (H 2 S) in the water column is produced from bacterial sulfate reduction once more powerful terminal electron acceptors, like oxygen, are depleted. Thought to be important to the occurrence of euxinic conditions are a warm climate and a nutrient-trapping geography, such as in restricted basins that were predisposed to euxinia . Many OAEs are also linked to episodes of enhanced volcanism and the emplacement of flood basalts, providing a mechanism to drive climate warming (via elevated atmospheric CO 2 concentrations; Naafs et al., 2016) and ocean deoxygenation, via reduced O 2 solubility and increased phosphate weathering and hence ocean productivity (Flögel et al., 2011;Monteiro et al., 2012;Wallmann, 2001). Plausible OAE triggering mechanisms are hence known, even if not necessarily explored yet in any great mechanistic and dynamical detail. This leaves the primary outstanding question: What is (are) the dominant mechanism(s) by which the Earth system recovers from such extreme conditions?
There are two main purported mechanisms to drown atmospheric CO 2 and cool and reoxygenate the oceans: (1) enhanced organic carbon burial in marine sediments and (2) enhanced silicate mineral weathering on the continents with subsequent carbonate burial in marine sediments (e.g., Berner, 1999;Ridgwell & Zeebe, 2005). However, geochemical evidence from OAE2, for example, suggests that the climate recovery was faster than can be explained by enhanced silicate and carbonate weathering alone (Kuypers et al., 1999;Pogge von Strandmann et al., 2013;Wallmann, 2001) and, thus, emphasizes the important role of OM Paleoceanography 10.1029/2018PA003470 burial in terminating anoxic events (Arthur et al., 1988). Enhanced burial of OM removes CO 2 from the ocean/atmosphere system and thus induces global cooling. At the same time the accumulation of O 2 in the atmosphere and surface waters has been shown to act as a negative feedback, eventually leading to the reoxygenation of the deeper ocean and a subsequent shutdown of the positive feedback between anoxia, enhanced P recycling and increased productivity (Handoh & Lenton, 2003;Tsandev & Slomp, 2009;Van Cappellen & Ingall, 1994). While not a true OAE per se, recovery from the warming transient at the Paleocene-Eocene boundary ("PETM") is also inferred to have been amplified by enhanced organic matter burial and not just driven by silicate weathering alone (Bowen & Zachos, 2010;Gutjahr et al., 2017). If so, the greatest outstanding unknown shifts to the specific means by which organic carbon burial is enhanced during OAEs.
Long-established explanations for increased OM burial during OAEs infer increased productivity and/or an increased preservation of OM (Demaison & Moore, 1980;Pedersen & Calvert, 1990). Their relative importance is most likely dynamic, but difficult to assess. Elevated marine productivity and carbon export from the ocean surface very likely played an important role for triggering and maintaining OAEs, potentially driven by an increased source of phosphate from continental weathering in conjunction with reduced removal due to efficient phosphate regeneration under anoxic conditions (Van Cappellen & Ingall, 1994). Indeed, previous Earth system modeling has been unable to explain the observed widespread occurrence of anoxic and euxinic conditions without invoking at least a doubling of ocean phosphate inventory and hence marine productivity (Monteiro et al., 2012;Ruvalcaba Baroni et al., 2014;Tsandev & Slomp, 2009). Yet, unlike preservation, increased productivity is not an efficient recovery mechanism as only a very small fraction of the carbon fixed in OM is sequestered in the sediments for longer time scales (Volk & Hoffert, 1985). The majority of OM is remineralized in the upper ocean, releasing CO 2 and intensifying ocean anoxia due to increased oxidant demand. In contrast, enhanced OM preservation increases the removal of fixed carbon from the ocean/atmosphere system as well as the net oxygen gain by the ocean and atmosphere. Enhanced OM preservation could thus have been critical to the recovery from OAEs. Previously, enhanced OM preservation during OAEs has been explained by reduced degradation rates under anoxic environmental conditions (e.g., Demaison & Moore, 1980;LaRowe & Van Cappellen, 2011), potentially supplemented by physical protection of OM as a consequence of higher clay mineral availability (Kennedy et al., 2002). A further but less widely recognized factor involves the reaction of labile OM compounds, such as functionalized lipids and carbohydrates, with reduced inorganic sulfur species (e.g., H 2 S), resulting in sulfur cross-linked macromolecular networks (Sinninghe Damsté et al., 1988) that are more resistant to bacterial degradation and thus exhibit an enhanced preservation potential .
Intense sulfurization of black shale OM has been previously described for sediments from Tarfaya Basin (e.g., Kolonic et al., 2002Kolonic et al., , 2005, Demerara Rise (e.g., Böttcher et al., 2006;Hetzel et al., 2009), Pont d'Issole (Raven et al., 2018), or the Jurassic Kimmeridge Clay Formation (e.g., Sinninghe van Kaam-Peters et al., 1997van Dongen et al., 2006). At all these sites, intense OM sulfurization is linked to enhanced OM burial. For instance, based on a comprehensive geochemical analysis of samples from the Kimmeridge Clay Formation, van Kaam-Peters et al. (1997 and van Dongen et al. (2006) showed that OM contents can be quantitatively linked to the degree of OM sulfurization. Kolonic et al. (2005) estimated that OAE2 mass accumulation rates at Tarfaya Basin would result in unrealistically high export production if one assumed a purely productivity controlled OM burial flux. They thus conclude that the enhanced preservation of sulfurized OM increased OM burial at the Tarfaya Basin site. Furthermore, using an inverse reaction-transport model approach, Arndt et al. (2006Arndt et al. ( , 2009 demonstrated that intensely sulfurized OM, found in black shale sequences at Demerara Rise, is characterized by extremely low OM degradability and, thus, high preservation efficiencies of more than 80% (sometimes close to 100%) of total deposited black shale OM. These model estimates corroborate the findings of earlier lab experiments with sulfurized sapropel OM revealing that this substrate is quasi unavailable for sulfate reducers, resulting in a close to 100% preservation (Moodley et al., 2005).
Despite observational evidence advocating the importance of OM sulfurization for local OM burial during OAEs, its impact on global biogeochemical cycling and climate and, thus, its significance for Earth system recovery from OAEs is essentially unconstrained. Early observations suggested that OM is sulfurized during early diagenesis over time scales ranging from decades to several thousand years (e.g., Wakeham et al., 1995;Werne et al., 2000). Although such slow OM sulfurization can support an enhanced OM preservation during burial in marine sediments, it cannot serve as an efficient recovery mechanism that rapidly sequesters OM, reduces atmospheric CO 2 , and alleviates water column anoxia/euxinia. However, labora-10.1029/2018PA003470 tory experiments provide evidence for more rapid OM sulfurization (on the order of weeks, e.g., Adam et al., 1998;Kok et al., 2000;van Dongen et al., 2003). In addition, more recent studies in the Cariaco Basin show that reduced sulfur can be rapidly incorporated into sinking OM within days or less (Li et al., 2011;Raven et al., 2016). Small changes in redox stratification, for instance, driven by changes in export production or ocean circulation could thus trigger large changes in sulfurization and OM sequestration (Raven et al., 2018). As a consequence, rapid water column OM sulfurization could be an initially very efficient mechanism for drawing down CO 2 and H 2 S and thus inducing the Earth system's recovery from OAEs. However, although sulfurization of OM is considered by some authors a globally significant process for preserving organic compounds (Arndt et al., 2009;Raven et al., 2018;, the dynamic interplay between primary productivity, ocean redox conditions, OM sulfurization, and its preservation in the sediments has yet to be examined. In this paper, we quantitatively explore the role of OM sulfurization as a global recovery mechanism from ocean anoxic events. More specifically, we test if the occurrence of euxinic conditions and the rates of OM sulfurization in the water column can be sufficient to create an appreciable enhancement of OM sequestration in marine sediments. We use an Earth system model to quantify the impact of this process on key ocean biogeochemical cycles and take as an example event, the archetypal Cenomanian-Turonian boundary OAE (OAE2, ∼93.5 Ma), characterized by high global temperatures caused by elevated atmospheric CO 2 (Friedrich et al., 2012), enhanced primary productivity (Schlanger & Jenkyns, 1976), widespread anoxia (e.g., Brumsack, 2006;Goldberg et al., 2016;Owens et al., 2012), and euxinia with free H 2 S in the water column (Monteiro et al., 2012;Pancost et al., 2004;).

Methods
For our numerical model experiments, we employ the cGENIE model, which is based on a 3-D-ocean circulation and a 2-D energy-moisture balance atmosphere model (Edwards & Marsh, 2005) and coupled to marine biogeochemical cycling of carbon, oxygen, phosphorus, and sulfur (Monteiro et al., 2012;Ridgwell et al., 2007). Here, cGENIE is implemented on a 36 × 36 equal-area horizontal grid with 16 vertical levels in the ocean and uses Cenomanian bathymetry and continental configuration identical to Monteiro et al. (2012) and Zhou et al. (2015, and see supporting information). Monteiro et al. (2012) evaluated the cGENIE setup for the Late Cretaceous and assessed the role of paleogeography, atmospheric CO 2 , and oceanic nutrient content on the distribution of ocean redox conditions. They showed that ocean phosphate content is the main lever for the model to successfully simulate observed patterns of ocean anoxia/euxinia for before and during OAE2. Further model evaluation is provided in Zhou et al. (2015), showing that simulated regional differences in oxygen levels in the proto-North Atlantic Ocean are consistent with the paleoredox proxy I/Ca. We implement a couple of changes compared to the configuration of Monteiro et al. (2012). First, we add a representation of OM sulfurization in the water column. In this, the more labile pool of sinking OM represented in cGENIE (see Ridgwell et al., 2007) is subject to sulfurization in the presence of free H 2 S, which we implement via a bimolecular rate law controlled by rate constant k sulf (see supporting information). Sinking sulfurized OM is hence transferred to the refractory pool and is not subject to further alteration or degradation in the water column and sediments. Note that this is an end-member scenario in that we do not consider here the possibility of partial remineralization of sulfurized OM occurring. However, our assumption is but supported by observations and model results (Arndt et al., 2006(Arndt et al., , 2009Moodley et al., 2005). Second, to simplify the analysis of sulfurization impacts, we omit nitrogen dynamics (and hence consider PO 4 as the only nutrient, as per Meyer et al., 2016). This slightly shifts patterns of the highest primary productivity, however our model still successfully reproduces observed patterns of ocean redox conditions at the seafloor and in the photic zone (see supporting information). Furthermore, we couple the ocean biogeochemistry model to the new 1-D vertically resolved and fully coupled reaction-transport sediment model, "OMEN-SED" (Hülse et al., 2018). Therefore, giving us the ability to avoid having to impose reflective boundary conditions for organic carbon and hence artificially amplifying oxygen consumption at the seafloor (Hülse et al., 2017).
OMEN-SED describes the degradation and burial of two OM pools in the sediments, as well as associated biogeochemical dynamics of the most important terminal electron acceptors (here O 2 , SO 2− 4 ) and methane (CH 4 ), and related reduced substances (H 2 S). For this study, we omit phosphorus dynamics in OMEN-SED and its regeneration under anoxic conditions. Instead, in our initial steady-state experiments, all organic phosphorus raining down on the seafloor gets instantaneously recycled and returned to the deepest ocean grid-cell as phosphate (PO 4 ), thus fixing the global PO 4 inventory. (We will address the role of longer-term  Table 1. Locations of observations (a), simulated TOC wt% at 1-m sediment depth for a pre-OAE2 simulation with 1 × [PO 4 ] (b) and peak-OAE2 with 2 × [PO 4 ] (c), both with observations superimposed and excluding OM sulfurization. Crossplots comparing simulated and observed TOC contents ((d) pre-OAE2; (e) peak-OAE2), with error bars indicating ranges of observed values reported for the specific location. TOC = total organic carbon.
feedbacks, e.g., involving PO 4 in a separate and subsequent paper.) Sediments also receive a nonbiogenic, noncarbonate detrital flux of 0.18 g·cm −2 ·k.y. −1 , a value characteristic of early Cenozoic open-ocean sediments (Zeebe & Zachos, 2007). In OMEN-SED, the degradation rate constant for labile OM follows the tuned spatially uniform value of Hülse et al., 2018Hülse et al., (2018  The entire sulfurized OM pool (together with the associated, scavenged H 2 S) is assumed to be completely preserved in the sediments.
In our initial assessment of the impact of OM sulfurization on marine biogeochemical cycling, we explore the sensitivity of the system to two main controlling parameters (sulfurization rate and ocean nutrient inventory) and create an ensemble of model experiments. These experiments are carried out with atmospheric xCO 2 fixed at 4.0× preindustrial values (1,112 ppmv, as determined by Monteiro et al., 2012) and atmospheric xO 2 at modern (21,000 ppm). We set no weathering flux and hence deliberately do not attempt to account for the role of longer-term (ca. 10-100 kyr) feedbacks between OM burial, CO 2 , ocean temperature, and oxygenation in these particular simulations. While organic carbon is preserved and buried in the sediments, we force a closure of the carbon cycle by means of the fixed atmospheric xCO 2 boundary condition (and hence preventing any climate feedback). All PO 4 in OM reaching the sediments, whether sulfurized or not, is returned to the ocean. While we allow the removal of OM associated H 2 S, this loss is small compared to the initial ocean sulfate inventory of 29,160 mol/kg (e.g., the total ocean dissolved sulfur changes W Shatsky and Hess Rise (DSDP sites 305, 310 Note. For locations reporting a range of TOC wt% the mean value is used for the model-data comparison, for locations with maximum or minimum TOC wt% this value is used and if evidence reports "low TOC" a value of 0.25 TOC wt% is used (as this is the minimum value reported in the database). A dash (−) indicates no observation for this site in the database. Observations adapted from Monteiro et al. (2012). For references and explanation of the evidence see Monteiro et al. (2012). TOC = total organic carbon; OAE2 = Ocean Anoxic Event 2.
less than 0.05% over the course of the experiment applying a moderate sulfurization rate constant of k sulf = 10 5 M −1 year −1 ). (We explore and discuss the role of feedbacks in an open system configuration subsequently.) We test a range of potential pre-OAE2 through peak-OAE2 nutrient inventories, by considering PO 4 inventories varying between ×1 and ×2 modern (in increments ×0.2 modern), following Monteiro et al. (2012). Simultaneously, for each of the different [PO 4 ] assumptions, we vary the OM sulfurization rate constant over a plausible range from k sulf = 10 3 M −1 year −1 , reflecting slow diagenetic sulfurization (Sinninghe Damsté et al., 2007) to k sulf = 10 7 M −1 year −1 representing reaction kinetics similar to those of fast secondary redox reactions (e.g., H 2 S reoxidation by O 2 ). We hence generate a 6 × 5 grid of experiments (30 in total) of differing [PO 4 ] versus k sulf with each experiment run for 20,000 years to achieve steady-state.

Simulated Black Shale Distribution
Inclusion of the coupled OMEN-SED model allows us to make prognostic projections of the preservation of OM in accumulating marine sediments in cGENIE, and hence the global distribution of black shales. Essentially, we are now able to forward proxy model key sedimentological features rather than having to make the assumption that black shale occurrence in the geological record can be equated with modeled occurrence of ocean floor anoxia, as was previously the state-of-the-art in paleo Earth system modeling (Monteiro et al., 2012). We find that even without including the influence of OM sulfurization and assuming a highly simplistic globally uniform, and open-ocean Paleogene relevant detrital accumulation rate in marine sediments, our simulated distributions of black shales are in good agreement with compiled observations (Figure 1 and Table 1). Prior to the onset of OAE2, total organic carbon (TOC) contents were already elevated in the proto-North Atlantic Ocean with contents that can reach up to 10 wt% locally (e.g., Kuypers et al., 2002;. Model results show comparable TOC contents in the proto-North Atlantic regions where seafloor anoxia is simulated in the model (Figure 1a and Table 1). Conversely, both observations and model results reveal low TOC contents (generally <1 wt%) for most of the Tethys Sea and especially shallow waters in the north, the proto-South Atlantic, as well as the Indian and Pacific Oceans (Figure 1a). Under peak-OAE conditions, increased TOC contents and laminated black shales occur globally ( Figure 1b and Table 1), with particularly high TOC contents reported in the proto-North Atlantic Ocean, the proto-South Atlantic, and the Tethys Sea. Simulated TOC contents reach >20 wt% in parts of the proto-Atlantic Ocean and the Tethys Sea, as well as most of the equatorial regions of the Pacific Ocean, qualitatively paralleling these observations. Conversely, no black shale deposition is simulated in the shallower northern Tethys Sea. The degree of general agreement between observations and model results has two reasons. First, the ocean model presumably projects plausible patterns of Cretaceous seafloor oxygenation state, an inference supported by previous work contrasting cGENIE simulations with the I/Ca oxygenation proxy (Zhou et al., 2015). This, in turn, leads to the absence of OM preservation in oxic conditions due to the higher OM degradation rate constant used in OMEN-SED. Second, for anoxic seafloor conditions the OM degradation rate constant in OMEN-SED has been tuned in order to match observed TOC contents with simulated values. Hence, while the magnitude of projected TOC wt% has been deliberately tuned, the spatial pattern largely arises from the distribution of environmental conditions projected by cGENIE for the Cretaceous.

The Role of Organic Matter Sulfurization
When we enable OM sulfurization, we find that reaction rates at the lower end of our considered range (i.e., k sulf ≤ 10 4 M −1 year −1 ) do not exert a significant influence on the redox state of the global ocean and OM burial. However, faster OM sulfurization rates (i.e., k sulf ≥ 10 5 M −1 year −1 ) reduce the availability of labile OM for degradation in the water column (because of its conversion to sulfurized, refractory OM), thus lowering oxidant demand, oxygen consumption, and H 2 S production (Figures 2a and 2b). The development of ocean euxinia is further suppressed as H 2 S is scavenged from the water column as a consequence of OM sulfurization (Figures 3f-3o). Overall, most OM sulfurization takes place in the upper 1,000 m of the water column (not shown) due to the generally high availability of H 2 S and labile OM for sulfurization at these depths.
Increasing the global PO 4 inventory and thus export production, generally enhances ocean anoxic/euxinic conditions as well as OM rain to and burial in the sediment (Figure 2). However, global export production decreases with increasing values of k sulf (Figure 2d) because less OM is degraded in the upper ocean, deepen- ing nutrient regeneration and making the return of recycled PO 4 to the surface ocean less efficient. However, total OM rain to, and preservation within, the sediments always increases with k sulf as the bulk OM pool becomes more refractory (Figures 2e-2g). We find that the effect of increasing [PO 4 ] is to make the system more sensitive to k sulf . The reasons for this increased sensitivity are threefold. First, when PO 4 availability is high, more labile OM is produced and thus available for sulfurization. Second, higher [H2S] in the water column promotes higher sulfurization rates and, third, a nutrient rich ocean can sustain export production even under enhanced OM sulfurization that limits PO 4 recycling. For simulated peak-OAE2 conditions (2 × PO 4 ), a moderate OM sulfurization rate of k sulf = 10 5 M −1 year −1 increases OM burial by ∼0.16 PgC year −1 compared to no sulfurization (k sulf = 10 7 M −1 year −1 increases burial by ∼1.14 PgC year −1 ). These values translate to a 30% increase or more than a tripling of OM burial (Figure 2f). Under these assumptions, the sedimentary OM preservation efficiency reaches ∼30% for k sulf = 10 5 M −1 year −1 and ∼55% for k sulf = 10 7 M −1 year −1 (Figure 2g).
We also find that enhanced OM preservation due to rapid OM sulfurization could lead to a significant reoxygenation of the ocean and a ca. 50% decrease in anoxic volume ([O 2 ] < 5 mol/kg, Figures 3b-3e). The associated reduction in water column euxinia is even more drastic, with the global euxinic volume ([H 2 S] > 10 mol/kg) dropping by about 80% (Figures 3g-3j). Our results show substantial spatial variability in water column integrated OM sulfurization rates, a consequence of the spatial distribution of H 2 S and labile OM in the model. For instance, highest water column integrated sulfurization rates are simulated in the proto-Atlantic Ocean (Figures 3m-3o) enhancing TOC contents by locally more than 15 TOC wt% compared to the no-sulfurization experiment (Figure 4). Observations generally confirm the occurrence of black shales with high TOC contents in these provinces ( Figure 1b and Table 1) and the three sections with observed OM sulfurization are located in areas with elevated sulfurization rates (Figures 3m-3o). Conversely, lowest water column integrated sulfurization rates occur in the high latitudes of the Pacific and Southern Ocean, a consequence of deep-water formation and thus the absence of euxinic conditions in these areas (Monteiro et al., 2012).

Significance of Fast Sulfurization for OM Burial During OAE2
Our Earth system model-based investigation supports the hypothesis that rapid OM sulfurization (i.e., k sulf ≥ 10 5 M −1 year −1 ) is important to the formation of Cretaceous black shales in productive and/or semirestricted, euxinic environments (Arndt et al., 2009;Kolonic et al., 2002;Raven et al., 2018). Assuming a plausible range of oceanic H 2 S concentrations between 10 and 100 M, a threshold sulfurization rate constant of 10 5 M −1 year −1 for the onset of globally important impacts translates to a first-order rate constant of 1.0-10.0 year −1 . This is two to five orders of magnitude higher than empirically determined rate constants for slow, diagenetic OM sulfurization (2·10 −4 year −1 , Werne et al., 2000; 1.3·10 −2 year −1 , Sinninghe Damsté et al., 2007). However, our rate constant is in good agreement with more recent studies in the Cariaco Basin, which show that reduced sulfur can be rapidly incorporated into OM in the water column within days or less (Li et al., 2011;Raven et al., 2016).
The impact of rapid OM sulfurization is to drive up to a threefold increase in global carbon burial rates, making this process potentially key to OAE recovery. We find that even relatively subtle changes in [PO 4 ] and hence export production, in combination with rapid OM sulfurization, have the potential to cause substantial changes in global OM burial (e.g., a 20% increase in [PO 4 ] from ×1 to ×1.2 driving an approximate doubling of OM burial at the highest rates of k sulf , Figure 2f). Furthermore, increased OM burial due to OM sulfurization is not globally uniform but has a strong local response (Figure 4). Simulated OM preservation efficiencies are specifically enhanced in euxinic areas of the ocean (depending on k sulf ranging from 70% to 100%, Figure 5) and compare well with predicted preservation efficiencies for OAE2 black shales from Demerara Rise (79% to 89%, Arndt et al., 2006Arndt et al., , 2009. Considering the implications of our results to the event as a whole-assuming even just a moderate sulfurization rate constant (k sulf = 10 5 M −1 year −1 )-we predict excess peak-OAE2 (×2 PO 4 ) OM burial of 0.630 PgC year −1 above the pre-OAE2 background (×1 PO 4 ) OM burial rate of 0.061 PgC year −1 . This translates to 5.25 × 10 13 mol C year −1 and is thus one order of magnitude higher than the excess burial rate of 0.32 × 10 13 mol C year −1 calculated by Arthur et al. (1988). However, it needs to be considered that our results are more likely representative for the maximum extent of ocean anoxia/euxinia, thus simulating maximum OM burial rates during the peak-OAE2 event. These maximum rates would not be sustained for the complete recovery phase from peak-OAE2 and then instantaneously drop to background rates. It is more likely that repetitive changes in the redox stratification, due to changes in primary productivity and/or ocean circulation, fostered periodically enhanced OM sulfurization and preservation in a highly dynamic ocean system. We expect that areas already characterized by low O 2 quickly tip over and become carbon burial hotspots.
Overall burial rates would gradually decrease from these maximum rates during ocean reoxygenation.

Transient Effects of OM Sulfurization
While our steady-state ensemble analysis allows us, for the first time, to quantify the magnitude and spatial patterns of OM sulfurization, it tells us little regarding how impacts might evolve on geological time scales (e.g., >10,000 years), nor about how important feedbacks induced by the introduction of this new mechanism are to marine carbon cycling. To address this, we present additional model analysis exploring the transient impacts of OM sulfurization. Sulfurization of organic matter effects ocean biogeochemistry and climate in three main ways: First, it has a direct impact via the enhanced burial of organic carbon (C org ) and hence atmospheric CO 2 drawdown and climate cooling. This should act as a negative feedback on sulfurization, as cooling will tend to drive higher oxygen solubility and hence a reduced euxinic area. Second, there is an indirect impact via the ocean alkalinity inventory, caused by the permanent burial of H 2 S. The burial of sulfurized OM prevents the reoxidation of the scavenged H 2 S, which would reverse the gain of alkalinity caused by sulfate reduction (SO 2− 4 reduced to H 2 S). If sulfate reduction is never reversed because the H 2 S is removed through sulfurization and burial, this represents a permanent alkalinity gain by the ocean. In turn, increasing ocean alkalinity drives a shift in carbonate chemistry, away from dissolved carbon dioxide (CO 2(aq) ) and toward carbonate ions (CO 2− 3 ), lowering atmospheric pCO 2 (e.g., Zeebe & Wolf-Gladrow, 2001). This also represents a negative feedback on sulfurization, again via cooling and higher oxygen solubility. Finally, burial of organic phosphorus (P org ) contained in sulfurized OM may decrease the oceanic nutrient content and thus organic matter export production. However, the result of this is more complex, because while lower organic matter export will tend to reduce the euxinic ocean volume, it will also drive higher atmospheric pCO 2 and hence warming-a positive feedback on euxinia and OM sulfurization. One of the important roles of mechanistic models in paleoceanography is clearly then to help elucidate the net impact and dynamics of multiple interacting feedbacks in the Earth system (e.g., see for an analogy Jickells et al., 2005;Ridgwell, 2003).
To explore these feedbacks, both in isolation and combination, we perform a series of "open"-system simulations, that is, where atmospheric CO 2 is not prescribed but can evolve in response to an imbalance between inputs (kerogen weathering) and outputs (OM burial). Atmospheric xO 2 is still restored to modern concentrations (21,000 ppm). Given the large reservoir of oxygen in the atmosphere (ca. 3.7 × 10 19 mol O 2 ) as compared to the maximum burial rate of organic carbon in our model (<1.4 × 10 14 mol C year −1 , Figure 2f), we do not expect this assumption to have any great impact on time scales of a few tens of thousands of years. We show the results of simulations run for 30,000 years and started from the closed-system peak-OAE2 simulation (×2 PO 4 ) without OM sulfurization (the "spin-up," compare, e.g., Figures 3b, 3g, and 3l). We prescribe a constant global kerogen weathering flux of 0.534 PgC year −1 in the open simulations which represents the diagnosed OM burial rate at the end of the spin-up. Any carbon burial simulated beyond this is a net sink of carbon from the system.
We run a total of four different sulfurization experiments (summarized in Table 2), for which we pick a moderate sulfurization rate constant of k sulf = 10 5 M −1 year −1 for each. Each sulfurization experiment includes a different combination of feedbacks ( Table 2). The first sulfurization experiment (#1) includes excess C org removal associated with the burial of sulfurized OM, but with all P org raining down on the seafloor still instantaneously recycled back to the overlying ocean. The second sulfurization experiment (#2) is the same, but now also accounts for P-loss by assuming that all P org in sulfurized OM is permanently buried in the sediments. The third and fourth experiments are the same as the first two, except now we eliminate the permanent ocean alkalinity gain associated with buried sulfurized OM in order to focus on the more direct effect of sulfurization on OM burial alone. We implement these additional two experiments (#3 and #4) by prescribing a negative flux of alkalinity from the sediments back to the ocean that is set equal to twice the loss of sulfur associated with sulfurized OM burial. In other words, we correct ocean alkalinity for the buried H 2 S, that would otherwise have been reoxidized in the water column back to SO 2− 4 and thereby removing (two moles of) alkalinity, but instead is scavenged and permanently lost form the system. In addition, we run a fifth experiment (#5) in which OM sulfurization does not occur and thus serves as a control.
For the sulfurization experiment #3 which only accounts for C org burial (red dashed line), we find that the initial response compared to the control experiment is, as expected, an increase in OM burial (Figure 6d). The enhanced loss of carbon drives a decline in atmospheric pCO 2 ( Figure 6c) and thereby causing global mean temperatures to decrease (Figure 6a + 6b). Cooler temperatures, together with reduced OM degradation rates resulting from a transfer of organic matter from the labile to the sulfurized pool as it settles through the water column, act to increase oceanic O 2 (Figure 6e) and induce further H 2 S decrease in addition to the scavenging associated with OM sulfurization (Figure 6g). The observed decrease in alkalinity (Figure 6i) reflects the decrease in oceanic H 2 S. However, overall, the decrease in oceanic SO 2− 4 reflects the dominance of permanent burial of S within sulfurized OM (Figure 6f). Higher global rates of OM burial as compared to rates of kerogen weathering explain the observed decrease in oceanic DIC ( Figure 6j). As expected, OM burial rates gradually decrease during ocean reoxygenation (Figure 6d + 6e).
In sulfurization experiment #4 which includes both C org and P org burial (green dashed line), oceanic PO 4 declines over time (to 1.36 × [PO 4 ] after 30 ka, Figure 6h) and thus causes a decline in export production. Due to the decreased oxidant demand in the water column, reoxygenation of the ocean occurs faster than in the C org -only experiment #3 (Figure 6e + 6g). In combination, this leads to OM burial rates decreasing faster from the maximum rate, eventually dropping below the control experiment ( Figure 6d) and thus the constant kerogen weathering flux. Therefore, when considering P org burial, atmospheric pCO 2 (and eventually oceanic DIC) increases above the control-experiment concentration, corresponding to increases in mean temperatures (Figure 6a + 6b).
When we now allow the ocean to permanently gain alkalinity (experiments #1 and #2), we see a dramatic and progressive rise in the oceanic alkalinity inventory (dotted lines in Figure 6i). This is a consequence of SO 2− 4 being progressively reduced to H 2 S in combination with a proportion of this H 2 S reacting with labile OM and being permanently buried. As a result of the alkalinity increase, the ocean experiences a higher capacity to sequester and store CO 2 , leading to a very rapid decline in atmospheric pCO 2 (Figure 6c) and mean temperatures (Figure 6a + 6b) which drives a more efficient ocean reoxygenation (Figure 6e + 6g). For a constant oceanic PO 4 inventory (experiment #1) this trend is sustained for the entire simulation (red dotted line). When also allowing for P org burial (experiment #2, green dotted line), as PO 4 is lost from the ocean, the decline in export production eventually dominates over the alkalinity gain, causing an increase in atmospheric pCO 2 and temperatures together with ocean reoxygenation.  Table 2. Our experiments highlight an unexpected result-that is, as H 2 S is permanently lost from the system and prevented from reoxidizing back to sulfate, alkalinity is progressively gained by the ocean. This drives a profound decrease in atmospheric pCO 2 on a time scale of just a few tens of thousands of years. While this effect can be ameliorated if matched by progressive loss of nutrients and declining productivity in the ocean, it should be noted that we have not considered the phosphate regeneration feedback under anoxic conditions here (Van Cappellen & Ingall, 1994), which would tend to reduce the P-loss and hence weaken the amelioration of alkalinity loss. Furthermore, while open in terms of carbon, alkalinity, and nutrients, these experiments lack an sedimentary cycle of calcium carbonate (CaCO 3 ). We are hence missing an additional negative feedback on alkalinity-carbonate compensation (Ridgwell & Zeebe, 2005)-in which as alkalinity and carbonate saturation increase, the preservation and burial of CaCO 3 in marine sediments would increase, removing alkalinity at twice the rate of carbon (as CaCO 3 ) and creating a tendency for atmospheric pCO 2 to increase.

Limitations and Outlook
Although mechanistically based and quantitative, our analysis should be considered only as a first step toward evaluating the potential role and impact of OM sulfurization associated with events such as OAE2. In our closed-system experiments we have taken a deliberately simplified and steady-state numerical model-based approach to simplify and focus the analysis and highlight the key underlying processes. In our open-system experiments we have focused on the intermediate-term (few tens of kyr) implications of sulfurized OM burial for CO 2 sequestration, climate cooling, ocean reoxygenation, and the nutrient and alkalinity inventory of the ocean. There are several important additional feedbacks that ultimately should be considered in the context of the long-term (100-1,000 kyr) dynamical evolution of the system, particularly (1) climate-responsive terrestrial weathering of P-bearing rocks, (2) CaCO 3 weathering and burial, and (3) an anoxia related PO 4 regeneration at the ocean-sediment interface. Furthermore, we have omitted consideration of N cycle dynamics that will serve to further modulate the impact of OM sulfurization. We are also missing the ultimate geological feedback on atmospheric pCO 2 and climate-silicate weathering, which would tend to damp the longer-term trends in pCO 2 and surface temperature although on a time scale longer than the pronounced changes we see occurring in our transient simulations (e.g., Lord et al., 2016). Again, Earth system models are needed to untangle the net results and dynamics of all the interacting feedbacks. Future work is planned considering these subcomponents, in isolation and ultimately all combined, in the context of OAE recovery dynamics.
The wider implications of our novel finding of increased ocean oxygenation and substantially reduced water column euxinia due to OM sulfurization warrants further discussion. Kump et al. (2005) showed that large H 2 S emissions from an euxinic ocean (>2,000 Tg S year −1 ) triggers a threshold effect causing H 2 S to build up rapidly in the atmosphere, leading to an increase in the atmospheric methane lifetime and thus a destruction of the ozone shield. Indeed, in a previous modeling analysis of potential late Permian conditions using an earlier version of the cGENIE Earth system model,  projected the occurrence of both widespread ocean euxinia and sulfur (as H 2 S) fluxes to the atmosphere far exceeding those from volcanism. The importance of this work is helping provide a "common" (land and ocean) kill mechanism for the end-Permian. However, we find that accounting for OM sulfurization in our Late Cretaceous configuration, not only has the potential to reduce ocean euxinia by ∼80% (by volume), but also reduce (by threefold to fourfold) H 2 S emissions to the atmosphere (Figure 2c). This is not to say that we can rule out an H 2 S kill mechanism for the end-Permian, but rather highlight that OM sulfurization, in removing a highly reactive sink for O 2 in the ocean (H 2 S), has more direct implications for upper ocean and atmosphere environmental conditions in addition to the secondary impacts via carbon burial and greenhouse gas warming that have previously been the main focus of interest in this mechanism. Better estimates of the impacts of OM sulfurization specifically for the end-Permian H 2 S kill mechanism would require using an end-Permian configuration of the model.

Conclusions
Our Earth system model simulations show that organic matter sulfurization in the water column can significantly enhance organic carbon burial during OAE2. The analysis demonstrates that rapid sulfurization has the potential to reduce ocean anoxia and euxina, draw down atmospheric CO 2 , and cool down the Cretaceous climate. We thus quantify for the first time by means of 3-D Earth system modeling the role of organic matter sulfurization and enhanced organic carbon burial as an emergency escape mechanism from extreme OAE conditions. Due to the simulated reduction of toxic H 2 S emissions to the atmosphere we propose that OM sulfurization could also have wider implications, for instance for the end-Permian H 2 S kill mechanism. We also highlight an additional mechanism and set of feedbacks induced by permanent H 2 S removal from the ocean-an increase in ocean alkalinity and with it lower pCO 2 and a cooler climate. Overall, the myriad of potential feedbacks that operate in a marine system comprising cycles only of carbon, alkalinity, oxygen, phosphate, and sulfur offer a clear answer to why numerical and particularly Earth system models have considerable utility in paleo studies.

Acknowledgments
We would like to thank two anonymous reviewers for very constructive comments which has led to a significantly improved and more interesting paper. D. H. and A. R. were supported by a Heising-Simons Foundation award. S. A. acknowledges funding from NERC (NE/I021322/1) and the European Union Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement 643052. A. R. also acknowledges support from EU grant ERC-2013-CoG-617313. The code for the cGENIE.muffin model is hosted on GitHub. The specific version used in this paper is assigned a DOI (https://doi.org/10.5281/zenodo.2575070). A corresponding user manual detailing installation, configuration, and tutorials is available from GitHub website (https://github.com/ derpycode/muffindoc). The data used in this paper are available in the supporting information.