Implementation of methane cycling for deep-time global warming simulations with the DCESS Earth system model (version 1.2)

. Geological records reveal a number of ancient, large and rapid negative excursions of the carbon-13 isotope. Such excursions can only be explained by massive injections of depleted carbon to the Earth system over a short duration. These injections may have forced strong global warming events, sometimes accompanied by mass extinctions such as the Triassic-Jurassic and end-Permian extinctions 201 and 252 million years ago, respectively. In many cases, evidence points to methane as the dominant form of injected carbon, whether as thermogenic methane formed by magma intrusions through overlying carbon-rich sediment or from warming-induced dissociation of methane hydrate, a solid compound of methane and water found in ocean sediments.

Abstract.Geological records reveal a number of ancient, large and rapid negative excursions of the carbon-13 isotope.Such excursions can only be explained by massive injections of depleted carbon to the Earth system over a short duration.These injections may have forced strong global warming events, sometimes accompanied by mass extinctions such as the Triassic-Jurassic and end-Permian extinctions 201 and 252 million years ago, respectively.In many cases, evidence points to methane as the dominant form of injected carbon, whether as thermogenic methane formed by magma intrusions through overlying carbon-rich sediment or from warming-induced dissociation of methane hydrate, a solid compound of methane and water found in ocean sediments.As a consequence of the ubiquity and importance of methane in major Earth events, Earth system models for addressing such events should include a comprehensive treatment of methane cycling but such a treatment has often been lacking.Here we implement methane cycling in the Danish Center for Earth System Science (DCESS) model, a simplified but welltested Earth system model of intermediate complexity.We use a generic methane input function that allows variation in input type, size, timescale and ocean-atmosphere partition.To be able to treat such massive inputs more correctly, we extend the model to deal with ocean suboxic/anoxic conditions and with radiative forcing and methane lifetimes appropri-ate for high atmospheric methane concentrations.With this new model version, we carried out an extensive set of simulations for methane inputs of various sizes, timescales and ocean-atmosphere partitions to probe model behavior.We find that larger methane inputs over shorter timescales with more methane dissolving in the ocean lead to ever-increasing ocean anoxia with consequences for ocean life and global carbon cycling.Greater methane input directly to the atmosphere leads to more warming and, for example, greater carbon dioxide release from land soils.Analysis of synthetic sediment cores from the simulations provides guidelines for the interpretation of real sediment cores spanning the warming events.With this improved DCESS model version and paleo-reconstructions, we are now better armed to gauge the amounts, types, timescales and locations of methane injections driving specific, observed deep-time, global warming events.

Introduction
Analyses of carbonate and organic carbon samples have revealed, with ever-finer time resolution, a number of large and rapid negative carbon isotope excursions since the Cambrian radiation of life about 540 million years ago.These isotope G. Shaffer et al.: Implementation of Earth system methane cycling excursions had typical δ 13 C amplitudes of > 2 ‰ and typical timescales of < 10 000 years.Some prominent examples of this have been found at about 56, 120, 183, 201, 252 and 260 million years ago (Hesslebo et al., 2000;Jenkyns, 2003;McElwain et al., 2005;Retallack and Jahren, 2008;Ruhl et al., 2011;Shen et al., 2011;Wignall et al., 2009;Zachos et al., 2001).Such large, rapid carbon isotope excursions can only be explained by large, rapid injections of light carbon (i.e., carbon-13 depleted) to the ocean-atmosphereland biosphere system, injections that forced global warming events for each of these excursions.As a starting point, the source of this carbon could be some combination of volcanism (δ 13 C ∼ −7 ‰), terrestrial biosphere reduction (δ 13 C ∼ −25 ‰), thermogenic methane input (δ 13 C ∼ −40 ‰) and methane hydrate release (δ 13 C ∼ −60 ‰), but the incorporation of further constraints often allows us to rule out the first two options as important contributors to most if not all of these events.
For example, if we assume an initial carbon pool in the ocean-atmosphere-land biosphere of about 40 000 Gt C (1 gigaton of carbon ,Gt C, is 10 15 g C) with a mean δ 13 C of 1 ‰ as at present, mass balance calculations then show that a large carbon isotope excursion of −5 ‰, like for example that of the end-Permian event (Shen et. al., 2011), could be explained by injections of ∼ 44 400, 8900, 5300 and 3400 Gt C of volcanic carbon, biosphere carbon, thermogenic methane or methane hydrate, respectively.The volcanic carbon explanation seems unlikely due to the required size and speed of the input.An explanation in terms of terrestrial biosphere die-off seems unlikely since it would require a biomass at least 4 times greater than the presentday biomass of around 2000 Gt C (Shaffer et al., 2008).Permafrost and peat carbon are other possible sources of terrestrial carbon.Little permafrost could exist for the warm conditions leading up to the warming events considered here (Shaffer et al., 2016).An explanation in terms of global peat conflagration (Kurtz et al., 2003) or degradation seems unlikely since it would require peat carbon reserves almost 20 times greater than present reserves of around 500 Gt C (Brigham et al., 2006).Therefore the most likely explanations for the carbon isotope excursions and associated global warming events involve thermogenic methane and/or methane hydrate inputs (Berner, 2002;Hesselbo et al., 2000;McElwain et al., 2005;Retallack and Jahren, 2008;Ruhl et al., 2011;Svensen, H. et al., 2004;Shaffer et al., 2016).
Thermogenic methane and carbon dioxide are produced when magma intrudes into overlying sediments containing old organic carbon.Many more such intrusions would be expected in association with increased volcanism that creates large igneous provinces and such provinces tend to correlate in time with the ancient carbon isotope excursions (Svensen et al., 2004).Magma contact with shale and coal favor production of CH 4 and CO 2 , respectively, whereas such contact with (isotopically heavy) limestone produces very limited amounts of CH 4 and CO 2 (Arnes et al., 2011).Shale is typically the main host rock through which magma intrudes in large igneous provinces (Svensen et al., 2004(Svensen et al., , 2007)).Methane hydrate is a solid compound of methane and water formed by sufficiently high pressure and low temperature and found mainly in ocean sediments.Within a hydrate stability zone, methane hydrate is formed when methane release from bacterial remineralization of organic matter exceeds that needed to sustain solubility levels in the face of vertical diffusive transport of methane.Warming will lead to methane hydrate destabilization and thereby feedback positively on the warming (Dickens et al., 1995).Methane hydrate dissociation was touted early on in a number of papers as an explanation for some deep-time carbon isotope excursions (e.g., Hesselbo et al., 2000;Berner, 2002;Kemp et al., 2005).
The discussion above underlines the need to consider input and ocean-atmosphere cycling of methane in modeling of such deep-time global warming events, yet in most cases models used to study such events have not included a full prognostic and interactively coupled approach incorporating methane into the climate as well as the biogeochemical cycles and isotope components.In addition to the question of the source of the light carbon inputs that force the events, other important aspects for understanding the workings of these events include the amounts, timings and locations of these inputs.For example, over what timescales is this carbon injected and how are the inputs partitioned between the ocean and the atmosphere?In order to address these questions, paleodata reconstructions should be compared with event simulations using Earth system models that consider methane input and cycling.For such a task, important aspects should be addressed like bacterial oxidation of methane in the ocean, air-sea gas exchange of methane, the dependency of atmospheric methane lifetimes on methane concentration there and the extreme radiative forcing for high atmospheric concentrations of methane.However, to our knowledge, such an Earth system model does not exist at present.Here we upgrade the Danish Center for Earth System Science (DCESS) Earth system model to fill this gap by designing appropriate methane input functions and by implementing methane cycling in the ocean, atmosphere and land biosphere.Since large, rapid inputs of methane to the ocean and oxidation of this methane there can lead to suboxic or anoxic oceanic conditions, we also upgrade our model to deal with ocean denitrification, nitrogen fixation and sulfate reduction (Shaffer, 1989).Finally we present some test simulations of deep-time global warming events using this extended model to illustrate ways forward for future model data analyses.

The DCESS model up to now
In its present configuration, the DCESS Earth system model (Shaffer et al., 2008) features modules for the atmosphere, ocean, ocean sediment, land biosphere and lithosphere (Fig. 1).The model has been designed to simulate global change on timescales of years to millions of years.It has been calibrated and tested against Earth system data and against output from more complex models.Model formulation is described in great detail in Shaffer et al. (2008); below follows a condensed description for present purposes.
The DCESS model consists of one hemisphere, divided by 52 • latitude into a low-mid-latitude (LL) and a high-latitude (HL) zone (Fig. 2).Sea ice and snow cover are diagnosed from estimated meridional profiles of atmospheric temperature (see below).Values for global reservoirs, transports and fluxes are obtained by doubling the hemispheric values.For the atmosphere module we use a simple, zonal mean, energy balance model for the near-surface atmospheric temperature, forced by yearly-mean insolation, outgoing longwave radiation, meridional sensible and latent heat transports and air-sea heat exchange (Fig. 3).In combination with the simple sea ice and snow parameterizations, the model includes the ice/snow albedo feedback and the insulating effect of sea ice.Prognostic equations for mean atmospheric temperature in each zone are obtained by integrating the energy balance over the zones.Meridional variations of atmospheric temperature are represented by a second order Legendre polynomial in sine of latitude, adjusted to fit zonal mean temperatures.Temperatures and gradients at the 52 • dividing latitude from this meridional profile are used in the calculation of meridional transports of sensible heat, latent heat and water vapor.
The atmosphere exchanges heat and gases with the ice-free part of the ocean (Fig. 3).Gases are taken to be well-mixed in the atmosphere and those considered are carbon dioxide and methane (for all three carbon isotopes), nitrous oxide and oxygen.Methane is oxidized to carbon dioxide on timescales on the order of 10 years, increasing as methane concentrations increase (see below).In this original model, air-sea exchange of methane is not considered.The model also deals with oxygen isotopes in atmospheric water vapor.Prognostic equations for atmospheric gases follow from internal sinks, exchanges with ocean surface layers and the land biosphere as well as input/outputs from weathering and lithosphere outgassing.
The model ocean is 270 • wide, extends from the Equator to 70 • latitude and covers 70.5 % of the Earth surface (Fig. 2).Both the LL and HL ocean sectors have 55 vertical layers with 100 m resolution.An ocean sediment sector is assigned to each of the ocean layers with widths representing modern day hypsography (Fig. 3).Ocean circulation and mixing is characterized by four parameters: 1.A transport, q, associated with HL sinking and the deepest LL upwelling.
2. Constant horizontal diffusion between the zones, K h , associated with the wind-driven circulation and deep recirculation.
3. Strong, constant vertical diffusion in the HL zone, K HL v , associated with high-latitude convection.

A weak, depth-dependent vertical diffusion in the LL
zone, K LL v .
Tracers considered are temperature, salinity, oxygen isotope in water, phosphate, dissolved oxygen, dissolved inorganic carbon for all three carbon isotopes, and alkalinity.In this original model, ocean methane is not considered.Biogenic production of particulate organic matter depends on surfacelayer phosphate availability but with lower efficiency in the HL zone.The calcite to organic carbon rain ratio depends on surface-layer temperature.Remineralization and dissolution of particles falling down out of the ocean surface are described by exponential functions with depth.Prognostic equations for each of the tracers in each of the ocean layers follow from ocean transport and mixing, internal sinks and sources, and exchanges with the ocean sediment.Ocean surface layers are also acted upon by exchanges with the atmosphere and river inputs associated with weathering.
The semi-analytical ocean sediment module considers calcium carbonate dissolution and both oxic and anoxic organic matter remineralization.The sediment is composed of calcite, non-calcite mineral and reactive organic matter.Sediment porosity profiles are related to sediment composition and a bioturbated layer 0.1 m thick is applied.Carbonate and organic carbon burial are calculated from sedimentation velocities at the base of the bioturbated layer.Bioturbation rates and oxic and anoxic remineralization rates depend on organic carbon rain rates and dissolved oxygen concentrations.The land biosphere module considers leaves, wood, litter and soil.Net primary production depends on atmospheric carbon dioxide concentration, and remineralization rates in the litter and soil depend on mean atmospheric temperatures.Methane production is a small, fixed fraction of the soil remineralization.The lithosphere module considers outgassing, weathering of carbonate and silicate rocks, and weathering of rocks with old organic carbon and phosphorus.Weathering rates are related to mean atmospheric temperatures.
Up to now the DCESS model has been applied mainly to long-term, future projections of climate and Earth system response to anthropogenic forcing (Shaffer et al., 2009;Shaffer, 2009Shaffer, , 2010)).Model results compare favorably with those of the best other Earth system models of intermediate complexity as seen in recent intercomparison studies of past, present and future carbon cycling and climate (Eby et al., 2013;Joos et al., 2013;Zickfeld et al., 2013).The results from these studies were also included in the recent Fifth Assessment Report of the Intergovernmental Panel for Climate Change (IPCC, 2013).From the 2013 studies onward, the carbon fertilization parameter of the terrestrial biosphere module was decreased somewhat in the model (Eby et al., 2013).
The Paleocene-Eocene Thermal Maximum (PETM) warming event 56 million years ago was addressed in the most recent model application (Shaffer et al., 2016) whereby the very simplified methane cycling sketched above was implemented and the carbon input was taken to occur over a 10 000-year timescale.That is significantly longer than timescales associated with methane oxidation in the atmosphere or ocean overturning, and the system response is rather insensitive to input type and location (ocean or atmosphere).Therefore, for simplicity the carbon input was taken to be in the form of carbon dioxide that was injected into the atmosphere.To deal with this deep-time warming event, several other modifications were made to the basic DCESS model (Shaffer et al., 2016): 1.The solar constant was decreased by 0.5 % from present day, sea level was raised by 100 m while retaining the modern day hypsography, and mean ocean salinity was taken to be 33.8 (the modern value is 34.7).
2. Background albedo was reduced due to less land at subtropical latitudes and more vegetation cover.
3. Initial land biomass was increased to reflect more land at low and mid-latitudes.
4. To account for high-latitude cloud forcing, extra forcing was imposed poleward of the latitude of sea ice extent in the standard model, pre-industrial steady state.
5. Parameters in the simple Budyko-type formulation for outgoing longwave radiation (Budyko, 1969) were adjusted in concert with the other forcing changes to achieve late Paleocene temperatures and climate sensitivity estimates.
6. Ocean calcium and magnesium concentrations were adjusted to late Paleocene values.
7. Model weathering formulation was extended to include the size of the weathering substrate, estimated for the late Paleocene to be less than that for the present.
8. Lithosphere outgassing was increased to conform to late Paleocene estimates.
9. The formulation of the calcite to organic carbon rain ratio was extended to also depend on the calcite saturation state of the ocean surface layer.
10.The dependence of carbon isotope discrimination during ocean surface-layer photosynthesis was extended to depend on concentrations of dissolved carbon dioxide and phosphate there.
3 Methane cycling implementation

Atmosphere
For an oxygenated atmosphere, the dominant sink for atmospheric methane is oxidation to CO 2 by reaction with the OH radical.Since this reaction depletes the concentration of these radicals, methane atmospheric lifetime, τ , grows as methane concentrations increase.In the original model (Shaffer et al., 2008), this effect and that of associated chemical reactions in the troposphere and stratosphere was addressed by fitting results from a complex atmospheric chemistry model (Schmidt and Shindell, 2003) to the equation where τ PI is the pre-industrial lifetime, a and b are fitting constants and M ≡ (pCH 4 − pCH 4,PI )/pCH 4,PI , where pCH 4 is the atmospheric partial pressure of methane and pCH 4,PI is its pre-industrial value taken to be 720 ppb.The Schmidt and Shindell (2003) target results and the model fit to them used in Shaffer et al. (2008) are shown in Fig. 4 (red dots and blue line).For this fit τ PI = 6.9 years, a = 0.96 and b = 6.6.Figure 4 also shows results from several additional chemistry modeling studies for the dependency of τ on pCH 4 (blue dots from Lamarque et al., 2006;black dots from Isaksen et al., 2011).For the range of low to medium methane concentrations considered in Schmidt and Shindell (2003) and Isaksen et al. (2011), these authors found rather low sensitivity of their results to atmospheric water vapor levels citing methane lifetime changes of less than 10 % for tropospheric water vapor increases of 30-40 %.On the other hand, Lamarque et al. (2006) found that atmospheric lifetimes decrease for extreme values of pCH 4 due to much more enhanced water vapor and OH production (and thereby greater methane oxidation) of a very warm climate.The red line in Fig. 4 shows a fit to all the chemistry modeling results of the figure using Eq.(1), a fit for which τ PI = 9.5 years, a = 0.78 and b = 11.This new fit now also captures, in part, enhanced OH production and the associated limitation of τ for very high values of pCH 4 and has been adopted in our enhanced model whereby the model atmospheric methane sink in moles per year is υ a pCH 4 /τ with υ a as the atmospheric mole content.
The standard radiative forcing formulation for methane (including spectral overlap with nitrous oxide) can be used for methane concentrations up to about 5 ppm (Myhre et al., 1998)  here.However, recent work has provided methane radiative forcing formulations, including nitrous oxide overlap, for methane concentrations up to 100 ppm (Byrne and Goldblatt, 2014a) and even up to 10 000 ppm (Byrne and Goldblatt, 2014b).We adopt these results and, for consistency, we also adopt the Byrne and Goldblatt (2014a) expressions for carbon dioxide radiative forcing and for carbon dioxide-nitrous oxide overlap.Figure 5 serves to illustrate these radiative forcings.Note that for very high methane concentrations, radiative forcing levels off and even decreases, an effect caused by enhanced shortwave absorption at such concentrations (Byrne and Goldblatt, 2014b).Taken together with the leveling off of methane atmospheric lifetimes at such concentrations (Fig. 4), this limits the role that can be played by methane in enhancing global warming under such extreme conditions.
To illustrate the effects of these new radiative forcing formulations we consider several specific atmosphere compositions and the warming associated with them.If we consider an atmosphere with pCO 2 , pCH 4 and pN 2 O values of 500, 10 and 0.5 ppm, respectively, this corresponds to radiative forcings relative to pre-industrial conditions of 3.26, 2.32 and 0.73 W m −2 but with a total forcing overlap of −0.16 W m −2 , leading to a total forcing of 6.15 W m −2 .With the original radiative forcing formulations used in the DCESS model (Myhre et al., 1998;Shaffer et al., 2008) the total radiative forcing is similar, 6.25 W m −2 .For a nominal climate sensitivity of 3 • C for a pCO 2 doubling (0.81 such total forcings would lead to global mean temperature increases of 4.98 and 5.06 • C, respectively, above the DCESS pre-industrial mean global temperature of 15 • C. If we consider an atmosphere with pCO 2 , pCH 4 and pN 2 O values of 1000, 100 and 1 ppm, respectively, this corresponds to radiative forcings relative to pre-industrial conditions of 7.45, 6.50 and 1.88 W m −2 but with a total forcing overlap of −0.93 W m −2 , leading to a total forcing of 14.97 W m −2 .With the original radiative forcing formulations, the total radiative forcing is considerably higher, 18.36 W m −2 .For a 3 • C climate sensitivity, such total forcing would lead to global mean temperature increases of 12.13 and 14.87 • C, respectively; for a 5 • C climate sensitivity (Shaffer et al., 2016) the increases would be 20.22 and 24.79 • C, respectively.

Simplified nitrogen and sulfur cycling
To be able to deal with suboxic and anoxic ocean conditions that would arise for significant oxidation of methane in the ocean interior, the model has been generalized to consider nitrogen and sulfur cycling with the introduction of nitrate (NO 3 ), ammonium (NH 4 ) and hydrogen sulfide (H 2 S) as additional ocean tracers (Shaffer, 1989; for simplicity, we often omit ion charges here).Furthermore, for simplicity we take the river input of nitrogen to be in the form of nitrate and to be equal to the (climate-and weathering substrate-dependent) river input of phosphate multiplied by r NP , the N : P Red- field (mole) ratio, taken to be 16.Likewise, organic matter produced, remineralized or buried maintains this ratio of nitrogen to phosphorus.Following Shaffer et al. (2008), we neglect dissolved organic matter such that the rate of export of particulate organic matter down out of the surface layer is equal to new production and we assume an exponential law for the vertical distribution of remineralization of "nutrient" and "carbon" components of this organic matter, each with a distinct e-folding length, λ N and λ C , respectively.
Model new production of organic matter in the surface layer has been generalized to depend on the limiting nutrient, either phosphate or nitrate, such that NP l,h = min(NPP l,h , NPN l,h ), where and l, h refer to LL and HL ocean zones, A l,hni 0 are the ice-free ocean surface areas, z eu is the surface-layer depth (100 m), r NP is a Redfield ratio (taken to be 16), sy is the number of seconds per year, PO l,h 4 and NO l,h 3 are the surface layer, phosphate and nitrate concentrations, and P 1/2 , N 1/2 are half saturation constants (1, 16 µmol m −3 ).L l,h f are efficiency factors, taken to be 1 and 0.36 for the LL and HL zones, respectively (Shaffer et al., 2008).This is how the model accounts for light and iron limitation in the HL zone.
The surface-layer sinks due to new production are −NP l,h for phosphate and −r NP NP l,h for nitrate.Other source/sinks and Redfield ratios are adopted from Shaffer et al. (2008).Since we will deal with water column denitrification, we also consider nitrogen fixation as a surface-layer source of nitrate when NO l,h 3 levels there drop below r NP PO l,h 4 or equivalently when NPP > NPN.In that case, we take the rate of nitrogen fixation, NF l,h , to be www.geosci-model-dev.net/10/4081/2017/Geosci.Model Dev., 10, 4081-4103, 2017 where NF 0 has been chosen to be 1 × 10 6 moles NO 3 per second.
In the water column, remineralization takes place through oxidation of organic matter with dissolved oxygen as long as oxygen levels are above a certain minimum value O 2,min , chosen here to be 3 mmol m −3 .Below that value, remineralization is assumed to occur by way of denitrification as long as nitrate levels exceed a certain minimum level NO 3,min , chosen here to be 0.03 mmol m −3 .The oxidation equation for denitrification is taken to be Here, as for oxidation with dissolved oxygen in the original DCESS model (Shaffer et al., 2008), we consider that organic matter formed in the ocean surface layer consists of proteins and lipids in addition to carbohydrates and adopt the mean composition proposed by Anderson (1995).Such a composition requires more oxygen for complete oxidation than carbohydrate alone would require (Anderson, 1995).This explains the enhanced nitrate sink in denitrification compared to standard Redfield stoichiometry (factors 94.4 and 47.2 compared with 84.8 and 42.4 in Reaction R1).
When oxygen concentrations are below O min and nitrate concentrations are below N min , remineralization is assumed to occur by way of sulfate reduction for which the oxidation equation is taken to be For the same reason as above, the sulfate sink is somewhat enhanced compared to standard Redfield stoichiometry (factor 59 compared with 53 in Reaction R2).Note that the carbon oxidation product in sulfate reduction is bicarbonate that, as opposed to the carbon dioxide produced in denitrification, contributes to alkalinity as well as dissolved inorganic carbon (DIC).Both reactions show a minor contribution to alkalinity by way of ammonia production.In addition, for sulfate reduction with our organic matter model, the alkalinity decrease from the production of hydrogen ions exceeds the alkalinity increase from the production of bisulfide (HS − ; by definition, alkalinity ∝ HCO − 3 +NH 3 +HS − −H + ; Dickson, 1981).In summary, for our organic matter stoichiometry, we find alkalinity : DIC = 0.151 and 1.038 for denitrification and sulfate reduction, respectively.
As shown in Reactions (R1) and (R2), NH 3 and HS − are produced by denitrification and/or sulfate reduction.As we show below, additional bisulfide is produced by microbes for anoxic conditions by way of anoxic methane oxidation (AMO).In the model, these species are oxidized to nitrate and sulfate, respectively, when transported by advection and diffusion to oxygenated ocean layers where O 2 ≥ O 2,min .These oxidation reactions are and Note that each of the reactions leads to an alkalinity decrease from production of hydrogen ions.Oxygenation rates, l,h ox,n (in moles per second), are expressed by where n is the layer number, V l,h n are the layer volumes, [NH 3 ] l,h n , HS − l,h n are the layer species concentrations and τ ox,NS is an ocean lifetime taken to be 200 days.
The simple ocean nitrogen and sulfur chemistry outlined above is intended to capture the most important processes that would occur in the ocean following a massive injection of methane without going into further details of the N and S cycles.Processes that are left by the wayside here include anaerobic ammonium oxidation (anammox), oceanic production of nitrous oxide and iron sulfide formation (Kuypers et al., 2003;Freing et al., 2012;Morse et al., 1987).Anammox has shown to be an important sink of nitrogen nutrients in the ocean but denitrification may still be the most important player here (Ward et al., 2009;Zehr and Kudela, 2011).Treatment of anammox and nitrous oxide would involve in-depth consideration of nitrogen cycling including nitrite production and consumption as well as approaches that transcend our simple, horizontally averaged model geometry.This is beyond the scope of the present paper but is being addressed by our group in other DCESS model developments in progress.Likewise, sulfide formation would require an explicit treatment of iron cycling, also beyond the scope of the model at present.

Methane oxidation and air-sea gas exchange
We now include methane as a tracer in the ocean model acted upon by advection, diffusion, air-sea gas exchange and microbial oxidation in the water column.Here we ignore any weak methane production that may take place in the surface ocean in anoxic microenvironments (Reeburgh, 2007) again with the goal of concentrating on massive methane injections.Oxidation of methane in oxygenated ocean layers (O l,h 2,n ≥ O 2,min ) produces carbon dioxide and consumes oxygen according to Layer oxidation rates of methane under these conditions, l,h ox,CH 4 ,n , are taken to be where [CH 4 ] l,h n are the layer methane concentrations and τ ox,CH 4 is an ocean methane lifetime for oxic conditions.
Geosci.Model Dev., 10, 4081-4103, 2017 www.geosci-model-dev.net/10/4081/2017/From a joint analysis of ocean methane and chlorofluorocarbons data, Rehder et al. (1999) found support for modeling methane oxidation as a first order process (Eq.5) and determined an oceanic methane lifetime of ∼ 50 years.A much shorter lifetime of ∼ 4 months was deduced for strong, fast and localized methane emissions during the 2010 Deepwater Horizon oil spill in the Gulf of Mexico (Kessler et al., 2011).Below we take τ ox,CH 4 to be 50 years since this value is likely more applicable to our global approach but in the following we also check the sensitivity of our results to the much shorter lifetime.
For suboxic/weakly anoxic conditions characterized by O 2 < O 2,min and NO 3 > NO 3,min , we adopt the following reaction for nitrate-dependent microbial AMO (Raghoebarsing et al., 2006;Cui et al., 2015): An alternative and slightly more energetic reaction involves nitrite, a species that we do not consider in our simplified nitrogen-cycling approach (Haroon et al., 2013).
For anoxic conditions for which O 2 < O 2,min and NO 3 < NO 3,min , sulfate-dependent microbial AMO (Treude et al., 2005;Cui et al., 2015) may be expressed as Note that whereas nitrate-dependent AMO (Reaction R6) does not produce alkalinity, the production of bicarbonate and bisulfide in sulfate-dependent AMO (Reaction R7) is a strong alkalinity source leading to an alkalinity : DIC ratio of 2, considerably higher than for sulfate reduction of organic matter (Reaction R2).From carbonate chemistry, an increase in ocean alkalinity tends to depress dissolved carbon dioxide concentrations in the ocean and, by way of air-sea gas exchange, also in the atmosphere.Layer oxidation rates of methane for sulfate-dependent AMO, l,h anox,CH 4 ,n , are taken to be where τ anox,CH 4 is an ocean methane lifetime for anoxic conditions (O 2 < O 2,min and NO 3 < NO 3,min ) that may be considerably longer than τ ox,CH 4 since Reaction (R7) is much less energetic than Reaction (R5) (Lam and Kuypers, 2011).
To take this into consideration and to be consistent with much lower anoxic remineralization rates compared to oxic rates in the DCESS model sediment module (Shaffer et al., 2008), we adopt τ anox,CH 4 = 0.1τ ox,CH 4 or 500 years as our standard case.On the other hand, nitrate-dependent AMO (Reaction R6) is energetically similar to oxic methane oxidation (Reaction R5) so we take τ ox,CH 4 for the methane lifetime associated with nitrate-dependent AMO.Air-sea gas exchange for methane, l,h CH 4 , is calculated as where the gas transfer velocities k l.h w are 0.39u 2 (Sc l,h CH 4 /660) 0.5 whereby u l,h are mean wind speeds at 10 m above the ocean surface, taken to be 8 m s −1 in both zones, and Sc l,h CH 4 are Schmidt numbers for methane that depend on ocean surface-layer temperature (see below).The methane solubility, η l,h CH 4 , was converted from the Bunsen solubility coefficients that depend on surface-layer temperature and salinity (Wiesenburg and Guinasso, 1979) to model units using the ideal gas mole volume.Furthermore, [CH 4 ] l,h is the dissolved methane concentration in the ocean surface layers.
In recent work it was shown that the traditional calculations of Sc CO 2 based on a third order polynomial of surfacelayer temperature (Wanninkhof, 1992) seriously underestimated Sc CO 2 for temperatures (T ) above 30 • C (Gröger and Mikolajewicz, 2011).These authors proposed the form Sc = a exp(−bT ) + c and fit that equation using experimental data of Jähne et al. (1987) to find best fit values for CO 2 of 1956, 0.0663 and 142.6 for a, b and c, respectively.Since we are dealing here with simulations of extreme global warming, we adopt this new Sc CO 2 expression and fit here, thereby superseding the traditional Sc CO 2 calculation previously used in our model (Shaffer et al., 2008).Gröger and Mikolajewicz (2011) did not consider methane but we directly fit their proposed exponential function to the experimental data of Jähne et al. (1987) for Sc CH 4 over 5 • C intervals in the temperature range 0-40 • C. We found best fit values (R 2 = 0.997) for a, b and c, respectively, such that Sc CH 4 = 1771 exp(−0.0650T ) + 128.8.We adopt this relationship here.
4 Sample simulations with methane cycling

Methane input specification
The rate of methane input, F CH 4 , we apply to the atmosphere-ocean system is taken to be where F tot CH 4 is the total input in Gt C, τ in is a timescale in years over which the first 75 % of the input is added and t is the time in years.This mathematical form has been used before within our group to characterize CO 2 input (Shaffer, 2009;Shaffer et al., 2016).As discussed in the Introduction section, the great bulk of both possible major F CH 4 sources to fuel deep-time warming events (methane hydrate and thermogenic methane) are likely to enter via the ocean.Estimates of the timescale for the carbon input driving the deep-time global warming events have varied from years to tens of thousands of years (eg.Zeebe et al., 2009;Cui et al., 2011;Wright and Schaller, 2013).On the lower end of this timescale range, Earth system response to F tot CH 4 will vary significantly depending on whether this methane is dissolved in the ocean or part, or even all, of it escapes as gas to the atmosphere, www.geosci-model-dev.net/10/4081/2017/Geosci.Model Dev., 10, 4081-4103, 2017 in bubble plumes for example (McGinnis et al., 2006).Such response variation may be explained in part by much shorter timescales for methane oxidation to carbon dioxide (in both the atmosphere and ocean) than ocean diffusive and overturning timescales.Escape to the atmosphere will cause immediate global warming since methane is a powerful greenhouse gas, much more so than carbon dioxide (Myhre et al., 1998).Oxidation of dissolved methane in the ocean will consume dissolved oxygen there (Reaction R5).For oxidation of large amounts of methane, suboxic or anoxic conditions may develop with complex nutrient and carbon cycling consequences, upon which we elaborate below.We address these questions here by taking the direct, undissolved methane escape to the atmosphere to be αF CH 4 and the methane dissolution in the ocean to be (1 − α)F CH 4 , whereby α is the fraction of total input that directly enters the atmosphere.Furthermore, we must specify how methane dissolution is distributed between the two ocean sectors and among the 55 layers of each sector.Here we simply assume that 84 % of the dissolution occurs in the LL sector and 16 % in the HL sector, a division that reflects model sector areas.Furthermore, we assume that the amount assigned to each sector is divided equally among the upper 30 model ocean layers, spanning the depth ranges over which methane hydrate and/or thermogenic methane inputs may be expected (e.g., Buffett and Archer, 2004;Svensen et al., 2004).

Pre-event steady state and a long sample simulation
In order to investigate the impacts of these model improvements on simulation results, we chose as an initial case study the late Paleocene (pre-PETM) configuration and associated model steady states constructed in Shaffer et al. ( 2016) but now including some changes like the adaption of the Byrne and Goldblatt (2014a) expressions for carbon dioxide radiative forcing and for carbon dioxide-nitrous oxide overlap.
In this configuration we assume a solar constant decrease of 0.5 % from present day, sea level rise of 100 m, mean ocean salinity of 33.8, slightly reduced background albedo (to represent less land in the subtropics and more vegetation cover), extra high latitude, longwave cloud forcing of 30 W m −2 , pre-industrial atmospheric oxygen level, adjustments in ocean calcium and magnesium concentrations and adjustments in initial lithosphere outgassing and weathering inputs.Furthermore, parameters in the model's simplified treatment of outgoing longwave radiation were adjusted to yield a global mean temperature of 25 • C and a climate sensitivity in the range 4-5 • C for atmospheric pCO 2 levels in the range 800-900 ppm (Shaffer et al., 2016).Open system steady states were then sought (by iteration and/or long time integration) such that external inputs (lithosphere outgassing and weathering) balanced external output (burial down out of the ocean sediment).Table 1 lists some properties of the model steady state we chose among various possibilities given the above constraints.This particular model configuration with a global mean temperature of 25.0 • C and an atmospheric pCO 2 of 860 ppm has a climate sensitivity of 4.8 • C. While this steady state largely reflects late Paleocene conditions to ease comparison with our prior work, we emphasize that a similar approach can be taken to design appropriate initial conditions for any particular deep-time warming event, like the end-Permian event (Shen et. al., 2011), based on paleo-reconstructions for conditions prior to that particular event.
Figure 6 shows time evolutions of some relevant model variables from a 200 000-year simulation that started in equilibrium with the above initial conditions and was forced by methane input of 4000 Gt C over a timescale (τ in ) of 3 kyr, whereby half of the input reaches the atmosphere and half is dissolved in the ocean (α = 0.5).Furthermore, to be specific it is assumed that the δ 13 C of the methane input is −40 ‰.Such forcing leads to maximum event temperature and pCO 2 rises of 5.7 • C and 1040 ppm, a peak atmospheric pCH 4 of 10.3 ppm and a LL calcite compensation depth (CCD) shoaling of 463 m (Fig. 6a, c and h).
The LL ocean warms by up to 5.5 • C at the surface, increasing to 6.0 • C at depth, reflecting modest polar amplification of the warming by about 20 % (Fig. 6b).Mean LL dissolved O 2 concentrations at 1000 m depth approach 0 after less than 2 kyr and remain so for an additional 3 kyr, driven, to a large extent, by oxidation of methane dissolved in the ocean (Fig. 6d).Considerable denitrification occurs during this 3 kyr long anoxic period, reflected in an associated decrease in new production due to nitrogen limitation (Fig. 6e).After this decrease, ocean new production then steadily increases to about 12 % higher than pre-event levels by 24 kyr (due to more ocean phosphate from enhanced weathering) and then slowly decreases again as weathering decreases with temperature.There is only a slight initial decrease of about 100 Gt C in land biosphere carbon inventory as more/less new production from CO 2 fertilization was roughly balanced by more/less soil remineralization from warming/cooling (Fig. 6f).The initial decrease in both ocean new production and land biomass enhance atmospheric pCO 2 and warming while the subsequent ocean new production increase followed by a slow decrease modulate the slow decrease in pCO 2 and slow cooling following peak warming.
Weathering increases over the event are accompanied by an initial drop in carbonate burial as dissolution decreases the amount of sediment CaCO 3 as well as the burial velocity (Fig. 6g and h).There is little initial change in organic carbon burial as the opposing effects tend to balance: decreasing burial velocity and increased organic matter preservation from decreasing O 2 concentrations.During this period there is a net source of carbon to the oceanatmosphere system (in addition to the prescribed methane input) as volcanic/weathering inputs exceed burial outputs.After the methane input event, both CaCO 3 and organic carbon burial increase in response to higher new and CaCO 3 production.Together with decreasing weathering in response to cooling, this results in a weak but long-lived net sink of carbon from the ocean-atmosphere system (Fig. 6g).In combination with decreasing pCO 2 and associated increasing ocean [CO 2− 3 ], the biogenic CaCO 3 production increase drives an overshoot of the LL CCD to depths greater than the pre-event depth, in agreement with theory and sediment core data (Dickens et al., 1997;Leon-Rodriguez and Dickens, 2010;Penman et al., 2016).This overshoot, which is even more pronounced at high latitudes, peaks by 35-45 kyr into the simulation (Fig. 6h).
The LL pelagic δ 18 O excursion in biogenic CaCO 3 is significantly muted compared to the benthic one due to enhanced δ 18 O in the surface layer (δ 18 O w ) from a more vigorous hydrological cycle during the warming event (Fig. 6i) The carbon isotope excursion is slightly muted in the atmosphere and enhanced in the LL ocean surface layer due to warming and temperature-dependent fractionation in airsea gas exchange (Fig. 6j).The LL marine organic matter isotope excursion is greater still, in agreement with paleoreconstructions (McInerney and Wing, 2011), as a consequence of model dependence of carbon isotope discrimination during ocean photosynthesis on concentrations of dissolved carbon dioxide and phosphate in the ocean surface layer (Shaffer et al., 2016).Use of a pCO 2 -dependent carbon isotope fractionation leads to a still greater isotope excursion in terrestrial organic matter, also in agreement with paleoreconstructions (Jahren and Schubert, 2013).Note that the model carbon isotope excursion is more protracted than the initial warming or the δ 18 O excursion related to the warming (Fig. 6a, i and j).Relaxation back to about half of maximum carbon isotope excursion values takes place over about 80 kyr compared to about 40 kyr for comparable temperature or δ 18 O relaxation.This can be explained by enhanced burial of δ 13 C-enriched CaCO 3 after the methane input event.After 200 kyr, most of model Earth system properties have relaxed back toward pre-event values.This long timescale is dictated by external input/output balances of the global carbon cycle (Fig. 6g), in particular by the slow, negative silicate weathering feedback on climate (Shaffer et al., 2008).

Sensitivity to input size, timescale and distribution
Figure 7 shows selected model results for the first 20 000 years of simulations like that of Fig. 6 but with methane input amounts of 2000, 4000 and 6000 Gt C. Both pCO 2 and especially pCH 4 increase more than linearly in response to the linear increase in methane input (Fig. 7a, c  and d).For pCO 2 this is due to positive feedbacks from the land biosphere and ocean.Initial warming from the methane injection leads to greater soil remineralization, land carbon inventory decrease and input of CO 2 to the atmosphere (Fig. 7f).Increasing methane oxidation in the ocean leads to denitrification that draws down ocean new production leading to increased ocean CO 2 outgassing that was already enhanced in response to decreased ocean CO 2 solubility from warming (Fig. 7e).The non-linear increase in pCH 4 is due mainly to the increase in CH 4 atmospheric lifetime with pCH 4 (Fig. 4).On the other hand, warming increases less than linearly (Fig. 7b) as expected from the radiative forcing dependencies shown in Fig. 5 and from our neglect here for simplicity of climate sensitivity increase with warming (Shaffer et al., 2016).Net carbon input from volcanism/weathering minus burial follows process sequences as described for Fig. 6 above (Fig. 7g).The CCD shoals as more methane is oxidized to CO 2 in the atmosphere and ocean leading to more acidic conditions that drive calcite dissolution in the sediment (Fig. 7h).However, CCD shoaling increase is damped for large methane dissolution in the ocean since this leads to anoxic conditions, methane oxidation with sulfate (Reaction R7) and associated alkalinity inputs that oppose dissolution.This also has had the effect of damping the atmospheric pCO 2 increase.In response to warming, oxygen isotopes of the LL ocean surface water are enriched by enhanced evaporation associated with a stronger atmospheric hydrological cycle (Fig. 7i).LL ocean surface-layer salinity (not shown Figure 6.Results from a 200 000-year simulation with a methane input of 4000 Gt C over a timescale of 3000 years and with half the input dissolved in the ocean and half escaping as gas to the atmosphere.(a) Global mean atmospheric temperature, (b) LL ocean temperature, (c) atmospheric pCO 2 and pCH 4 , (d) LL ocean-dissolved O 2 concentration, (e) ocean new and biogenic CaCO 3 production, (f) land biosphere carbon inventories for leaves + wood, litter + soil and their total, and (g) external input/outputs of carbon.The black line is volcanic/weathering inputs minus burial outputs.(h) Depths to the LL 10 % (LL CCD) and HL 50 % CaCO 3 wt%, (i) LL excursions of δ 18 O in biogenic CaCO 3 formed and those excursions corrected for excursions in ambient δ 18 O w .(j) Carbon isotope excursions for the atmosphere, the LL ocean, LL marine organic carbon and terrestrial organic carbon (see text; the methane input has δ 13 C = −40 ‰).same reason and the salinity excursion can be approximated very well by multiplying the results in Fig. 7i by 2.8.Atmospheric carbon isotope excursion values are directly proportional to the size of the methane input and exhibit a long decay timescale as noted in Fig. 6 above (Fig. 7j): Figure 8 shows selected model results for the first 20 000 years of simulations like that of Fig. 6 but with methane injection times τ in of 300, 1000, 3000 and 10 000 years.For the shortest injection times, input rates in terms of carbon equal or exceed those of present-day anthropogenic carbon emissions exceeding 10 Gt C yr −1 (Fig. 8a; IPCC, 2013).Global warming spikes to much higher val-ues for short injection times in response to large instantaneous methane radiative forcing and under the influence of slow ocean exchange times (Fig. 8b and d).Warming tops out at 8.0, 6.7, 5.7 and 4.8 • C for the shortest to the longest injection times.Atmospheric pCO 2 is also enhanced for shorter injection times in part due to the positive feedbacks from the land biosphere and ocean noted above (Fig. 8c, e  and f).Shorter injection times lead to more widespread suboxic/anoxic conditions as evidenced by initial decreases in ocean new production in response to denitrification (Fig. 8e).These conditions are a direct response to increased methane oxidation in the ocean due to enhanced methane input rates .Results for 20 000-year simulations for different methane input amounts over a timescale of 3000 years and with half the input dissolved in the ocean and half escaping as gas to the atmosphere.(a) methane input rate, (b) global mean atmospheric temperature, (c) atmospheric partial pressure of carbon dioxide, (d) atmospheric partial pressure of methane, (e) total ocean new production, (f) total land biosphere biomass, (g) volcanic/weathering inputs minus burial outputs (as black line in Fig. 6g), (h) LL CCD depth (10 % CaCO 3 dry weight in sediment), (i) oxygen isotope excursion of LL ocean surface-layer water, and (j) atmospheric carbon isotope excursion (the methane input has δ 13 C = −40 ‰).
(Fig. 8a).One consequence of this situation is the muted CCD shoaling for shorter injection times following from enhanced methane oxidation with sulfate (Reaction R7) and associated alkalinity inputs that oppose CaCO 3 dissolution (Fig. 8h).Both net carbon input and especially oxygen isotopes from LL ocean surface-layer water follow global warming due in large part to climate-dependent weathering and temperature-dependent atmospheric water vapor transport (Fig. 8a, g and i; Shaffer et al., 2008).The atmospheric carbon isotope excursion is much enhanced for short injection times.This can be understood as follows: the carbon isotope signal from the methane injected to the atmosphere, and transferred there to CO 2 via oxidation of the methane, builds up in the atmosphere and in the ocean surface layer (not shown) due to long exchange times of about 1 kyr with the deep ocean (Fig. 8j).The atmospheric isotope excursion tops out at −5.24, −3.80, −3.21 and −3.03 for the shortest to the longest injection times.
Figure 9 shows selected model results for the first 20 000 years of simulations like that of Fig. 6 but with 0, 50 and 100 % of the methane input reaching the atmosphere (α = 0, 0.5 and 1).Global warming, atmospheric pCH 4 and, to a lesser extent, atmospheric pCO 2 are enhanced for more input directly to the atmosphere (Fig. 9b-d).Warming is enhanced by 1.7 • C for all methane input to the atmosphere compared to all input dissolved in the ocean.Even with all methane input dissolved in the ocean, atmospheric pCH 4 increases initially by more than 2 ppm as methane outgasses .Results for 20 000-year simulations for different methane input timescales for an input of 4000 Gt C with half the input dissolved in the ocean and half escaping as gas to the atmosphere.Properties plotted in (a-j) as in Fig. 7. from the ocean via air-sea gas exchange (Fig. 9d).The positive land biosphere feedback on atmospheric pCO 2 is enhanced for all methane to the atmosphere due to enhanced warming (Fig. 9f).On the other hand, ocean feedbacks on atmospheric pCO 2 are more nuanced.The positive solubility feedback is enhanced for all methane to the atmosphere due to enhanced warming whereas the positive feedback from reduced ocean new production is enhanced for all methane to the ocean (Fig. 9e).The latter effect is due to methane oxidation in the ocean leading to denitrification and nitrogen limitation of new production.Increased CaCO 3 dissolution in the ocean sediment is a further consequence of methane oxidation in the ocean when all methane is dissolved there.This leads to decreased calcite burial and enhanced net carbon input to the ocean-atmosphere system despite less weathering input from less global warming here (Fig. 9g).This also explains the enhanced CCD shoaling for all methane dissolved in the ocean (although modulated somewhat from increased methane oxidation with sulfate and associated alkalinity inputs as discussed above; Fig. 9h).Again, oxygen isotopes from LL ocean surface-layer water follow global warming due to temperature-dependent atmospheric water vapor transport (Fig. 9i).The atmospheric carbon isotope excursion peaks slightly earlier for all methane to the atmosphere but, due to the long timescale of the methane input (3000 kyr) relative to ocean overturning times, the excursion amplitude is essentially the same for all cases (Fig. 9j).The flow diagram presented in Fig. 10 serves to compare and contrast the pathways and feedbacks associated with (1) methane input to the ocean only and (2) methane input to the atmosphere only.These are the cases plotted as blue and maroon lines, respectively, in Fig. 9.
We also carried out a simulation like the case in Fig. 9 for all methane dissolving in the ocean (blue lines) but for a much shorter oceanic methane lifetime of 4 months (compared to the standard lifetime of 50 years; see Sect.3.2.2).
Figure 9. Results for 20 000-year simulations for different fractions of methane escape to the atmosphere for a methane input of 4000 Gt C over a timescale of 3000 years.Properties plotted in (a-j) as in Fig. 7. Also shown are results from a simulation with all methane input to the ocean (0 % fraction to atmosphere) and a much shorter ocean methane lifetime of 3 months (dashed blue lines; see main text).
Results for this new simulation are plotted as dashed blue lines in Fig. 9.The shorter lifetime leads to much reduced methane concentrations in the ocean and less outgassing to the atmosphere (Fig. 9d).Furthermore, faster methane oxidation in the ocean decreases dissolved oxygen there slightly, leading to slightly enhanced anoxia, sulfatedependent anoxic methane oxidation and alkalinity inputs.As a consequence, CCD shoaling is damped somewhat as is the atmospheric pCO 2 increase (Fig. 9c and h).Together with less methane outgassing, this leads to a reduction of maximum global warming of about 7 % compared to the standard case (Fig. 9b).

Modeled ocean distributions
Figure 11 shows selected LL model ocean results for the first 20 000 years of a simulation forced by a methane input of 6000 Gt C with δ 13 C = −40 ‰ over a timescale (τ in ) of 3 kyr, whereby all the input is dissolved in the ocean (α = 0).This forcing configuration was chosen to highlight the workings of the model for suboxic/anoxic conditions and leads to maximum event temperature and pCO 2 rises of 5.7 • C and 1095 ppm, a peak atmospheric pCH 4 of 6.2 ppm and a LL CCD shoaling of 389 m (not shown).
LL ocean warming is greatest at depth, a reflection of polar amplification (Fig. 10a).Maximum warming there (6.2 • C) occurs about 1500 years after maximum surfacelayer (or atmosphere) warming.Subsurface methane concentrations build up to more than 30 mmol m −3 (Fig. 11b).Oxidation of this methane with dissolved oxygen (Reaction R5) forces suboxic conditions (O 2 < 10 mmol m −3 ) for thousands of years at intermediate and mid-depths of the LL ocean (Fig. 11c).Furthermore, a combination of reduced solubility due to warming and enhanced remineralization due to increased ocean new production (Figs.5-7) maintains suboxic conditions at intermediate depths there well beyond the 20 000 simulation years shown here.As dissolved oxygen is forced below O 2,min (3 mmol m −3 ), methane oxidation proceeds by way of nitrate-dependent anoxic methane oxidation (Reaction R6), essentially eliminating nitrate for several thousand years at intermediate and mid-depths of the LL ocean (Fig. 11d).Over a period of several thousand years after most methane is oxidized, nitrate concentrations recover to and even above pre-event values.This can be explained by a combination of nitrogen fixation (Eq. 3) and enhanced new production due to increased ocean phosphate concentrations from warming-enhanced weathering.As nitrate is forced below NO 3,min (0.03 mmol m −3 ), methane oxidation proceeds by way of sulfate-dependent anoxic methane oxidation (Reaction R7).This leads to the production of hydrogen sulfide that reaches concentrations over 10 mmol m −3 at intermediate depths almost 3 kyr into the simulation (Fig. 11e).
Figure 11f shows the model LL time-space distribution of omega ( ) defined as the ratio of carbonate ion concentration to carbonate ion saturation concentration for calcite.Calcite dissolution in the sediment increases as decreases from 1 and model biogenic calcite production in the surface layer is proportional to ( − 1) / {1 + ( − 1)} for ≥ 1 and 0 for < 1 (Shaffer et al., 2016).In response to methane oxidation to CO 2 , carbonate ions and thereby initially decrease in intermediate and mid-depths where this oxidation occurs.However, as dissolved oxygen and nitrate are con-sumed and sulfate-dependent anoxic methane oxidation takes over about 2 kyr into the simulation, increases in these layers in response to alkalinity produced in this reaction.After the period of sulfate-dependent anoxic methane oxidation, decreases again followed by a slow increase driven by decreasing CO 2 levels from carbonate compensation.LL oxygen isotope excursions in biogenic carbonate produced in situ (Fig. 11g) track temperature changes in the deep ocean but are reduced relative to these changes near the surface due to enhanced ambient water isotopes (Fig. 9i).After several thousand years, LL carbon isotope excursions in biogenic carbonate produced in situ (Fig. 11h) are slightly depressed in the deep ocean.This is due to dilution from additional less-depleted carbon inputs associated with sediment CaCO 3 dissolution.
Figure 12 shows selected LL and HL model ocean profiles at 0, 3 and 5 kyr into the simulation of Fig. 11 and serves to illustrate differences between the LL and HL ocean zones.Although methane input is scaled to be comparable per unit area for both zones, much higher methane concentrations build up in the LL zone than in the HL zone (Fig. 12a).Much higher HL vertical exchange spreads methane vertically and promotes some outgassing to the atmosphere there.Phosphate concentrations increase with time in response to warmingenhanced weathering (Fig. 12b), explaining the slow increase in ocean new production shown above .Despite suboxic and eventually anoxic conditions in the LL zone, the HL zone remains oxygenated due to stronger vertical ex-Figure 11.LL ocean property results as functions of water depth and time for a 20 000-year simulation with a methane input of 6000 Gt C over a timescale of 3000 years, all of which dissolves in the ocean.(a) Temperature anomaly, (b) methane concentration, (c) dissolved oxygen concentration (model denitrification for O 2 < 3 mmol m −3 ), (d) nitrate concentration (model sulfate reduction for NO 3 < 0.03 mmol m −3 ), (e) hydrogen sulfide concentration, (f) omega, the ratio of carbonate ion concentration to the carbonate ion saturation concentration for calcite (see Fig. 11g), (g) oxygen isotope excursion in biogenic carbonate produced in situ.(h) Carbon isotope excursion.No corrections for the carbonate effect have been applied in (g) and (h) since these corrections were only derived for pelagic foraminifera (Spero et al., 1997).
change and the downwelling branch of the overturning circulation there (Fig. 12c).As a consequence, HL nitrate concentrations remain elevated and there is no anoxic remineralization nor anoxic methane oxidation there (Fig. 12d).This also explains the essential lack of HL ammonium and hydrogen sulfate, < 0.02 and < 0.1 mmol m −3 , respectively (Fig. 12e  and f).Ammonium or hydrogen sulfate reaching the HL zone via horizontal exchange with the LL zone is quickly oxidized.By way of carbonate chemistry, addition of CO 2 to the water column via methane oxidation depresses carbonate ion concentrations (Fig. 12g).In terms of (see above) this reduction is relatively stronger in the HL surface layer, tending to depress biogenic carbonate production there.However, this reduction is largely offset by warming enhancement.Addition of CO 2 to the water column via methane oxidation increases ocean acidity (Fig. 12h).For the simulation considered here, surface-layer pH is reduced by 0.26 and 0.28 for the LL and HL zones, respectively, from pre-event levels already 0.47 and 0.48 lower than model pre-industrial levels (Shaffer et al., 2008) Figure 12.Low-mid-latitude (LL) and high-latitude (HL) ocean property profiles at selected times for the simulation of Fig. 11.LL and HL results are plotted with red and blue lines, respectively.Line types for the selected times after the simulation start (0, 3 and 5 kyr) are defined in (a).CO sat 3 is the carbonate ion saturation concentration for calcite.
turbated sediment layers, each 10 cm thick, as well as sediment properties exported down out of these layers.In the original model, synthetic sediment cores (SSC) produced in this way only contained information on calcite and organic matter content of the buried sediment (Shaffer et al., 2008).
We have extended this treatment to include carbon and oxygen isotopes and report the first results of this here for bulk carbonate, in the form of calcite in the model.For each sediment layer we consider conservation of calcite as well as conservation of carbon and oxygen isotopes in this calcite.This involves tracking the time-dependent inputs (from the rain of biogenic calcite produced in the ocean surface layer) and outputs (dissolution within the sediment layer and export down out of the layer by burial).Isotope effects for changes in surface-layer CO 2− 3 are applied in the production of bio-genic calcite (Spero et al., 1997).The new model extension also accounts for effects of possible "mining" of buried sediment (i.e., upward directed "burial" velocities) in response to very strong dissolution events (Shaffer et al., 2008).For this, properties buried in each SSC are recorded to provide correct values for properties reentering the active sediment layer from below during any such "mining" event.
Figure 13 shows bulk carbonate results from model LL and HL SSCs at 3000 m depth together with ocean surface-layer properties, all for the simulation of Fig. 10.Burial velocities and sediment calcite content are greater in the HL zone compared to the LL zone due to considerably greater model HL rain rates per unit area (Fig. 13a and b).This can be explained in part by temperature-enhanced HL rain rates in the model: the pre-event HL ocean surface layer is almost 10 • C Figure 13.Comparison of bulk carbonate results from model synthetic ocean sediment cores at 3000 m depth with ocean surface-layer properties for the simulation of Fig. 11.Low-mid-latitude (LL) and high-latitude (HL) results are plotted with red and blue lines, respectively, and the methane input has δ 13 C = −40 ‰.(a) Burial (downward) velocity relative to the base of the bioturbated sediment layer vs. synthetic sediment depth (SSD).SSD is referenced to zero at the start of the simulation and increases downcore by convention, (b) carbonate dry weight fraction vs. SSD, (c) LL bulk carbonate δ 13 C excursion vs. LL SSD (solid line).Also shown are the δ 13 C excursions in biogenic carbonate in the LL ocean surface layer (dashed line) and this δ 13 C including the correction for change in surface-layer CO 2− 3 (dotted line; Spero et al., 1997), both plotted vs. simulation time (right side y axis).The time axis is "stretched/squeezed" relative to the SSD axis, a product of time-varying burial rates, (d) LL bulk carbonate δ 18 O vs. LL SSD (solid line).Also shown are (1) ocean surface-layer temperature anomaly (dashed line) on a temperature scale (top x axis) related to the oxygen isotope scale (bottom x axis) using temperature dependence of biogenic carbonate δ 18 O (Bemis et al., 1998), (2) δ 18 O in biogenic carbonate produced in the LL ocean surface layer (dashed-dot line) that includes the water δ 18 O excursion (see Figs. 7-9i) and (3) this δ 18 O including the correction for change in surface-layer CO 2− 3 (dotted line; Spero et al., 1997), all plotted vs. simulation time, (e) as (c) but for the HL ocean and sediment, (f) as (d) but for the HL ocean and sediment.
warmer than for the pre-industrial simulation (Shaffer et al., 2008).For this reason and due to slightly higher new production (Table 1), pre-event HL biogenic calcite rain is more than 3 times greater than for the pre-industrial case.On the other hand, pre-event, LL biogenic calcite rain is only about 15 % greater than for this case.Over the first 20 kyr of the simulation shown here, 38.7 cm of sediment is laid down at 3000 m depth in the HL zone compared to 14.0 cm in the LL zone.Burial velocities and sediment calcite content decrease initially in response to increased dissolution and, to a lesser extent, decreased rain rate, as a consequence of the methane input and its oxidation in the ocean.Later on, both properties increase to exceed their pre-event levels in response to increasing new production from increasing ocean phosphate levels.
In both model zones the carbon isotope excursion in the SSC bulk carbonate is modestly attenuated relative to the ocean surface-layer excursions (Fig. 13c and e).This is due to reservoir timescales of up to several thousand years for the active sediment layer and, to a lesser extent, to the CO 2− 3 isotope effect (a reservoir timescale is the calcite inventory of a sediment layer divided by rate of calcite inwww.geosci-model-dev.net/10/4081/2017/Geosci.Model Dev., 10, 4081-4103, 2017 put to it).The reservoir effect also leads to somewhat distended bulk carbonate excursion maxima centered 3-4 kyr after surface-layer maxima (note that the time axis of the figure is "stretched/squeezed" relative to the sediment length axis, a product of time-varying burial rates).Similar conclusions and interpretations also apply to the oxygen isotope excursions in the SSC bulk carbonate and the ocean surface layers in both model zones (Fig. 13d and f).Note that if, in lack of further knowledge, a linear timescale would be assigned to SSC length, the bulk carbonate excursions would appear sharper and shorter, particularly in the HL zone.When used as a paleothermometer, SSC bulk carbonate in the LL zone severely underestimates the model temperature excursion in the LL surface layer, primarily due to the positive ambient water δ 18 O excursion demonstrated above, but with a minor contribution from the CO 2− 3 isotope effect (Fig. 13d).On the other hand, the slightly negative ambient water δ 18 O excursion in the HL surface layer and the CO 2− 3 isotope effect there work to improve the estimate of the surface-layer warming from the δ 18 O excursion of SSC bulk carbonate, albeit with the 3-4 kyr time lag due to the reservoir effect (Fig. 13f).

Discussion and conclusions
Here we have extended the DCESS Earth system model to include methane cycling.This is a necessary step for dealing with deep-time global warming events, some of which correspond to major life extinction events, since most of these warming events were probably forced by massive methane inputs (Berner, 2002;Hesselbo et al., 2000;McElwain et al., 2005;Retallack and Jahren, 2008;Ruhl et al., 2011;Svensen et al., 2004;Shaffer et al., 2016).To be able to treat impacts of such inputs more correctly and consistently, we have also extended the model to deal with suboxic/anoxic conditions in the ocean and their consequences.For this we have now included denitrification, nitrogen fixation, sulfate reduction and nitrate/sulfate-dependent anoxic methane oxidation.Furthermore, we have upgraded the treatment of methane for high concentrations in the model atmosphere with the latest radiative forcing relationships and with improved relationships for atmospheric lifetimes.To our knowledge, no other Earth system model of any degree of complexity has yet been formulated that can deal as comprehensively with extreme methane inputs and their Earth system consequences.
After formulating the model extensions we embark on extensive tests of model behavior for methane inputs of various sizes, timescales and locations.We demonstrate model behavior over event timescales exceeding 100 kyr but concentrate on the first 20 kyr of the simulations.The long, > 100 kyr, simulation demonstrates how warming-driven weathering increases lead to a slow buildup of ocean phosphate concentrations and ocean new production.This modulates the evolution of atmospheric pCO 2 and helps to explain model CCD deepening overshoot after the initial CCD shoaling from the methane input.Extensive ocean anoxia develops for larger methane inputs over shorter timescales with more methane dissolving in the ocean.For such anoxia, there is much sulfate-dependent anoxic methane oxidation that produces hydrogen sulfide but also alkalinity that works to oppose calcite dissolution in the sediment and pCO 2 rise in the atmosphere.Furthermore, extensive denitrification also occurs that initially depresses ocean new production, leading to pCO 2 outgassing, until nitrogen fixation steps in to fill the nitrate gap.The global warming excursion is greater when more methane escapes to the atmosphere, leading to higher pCH 4 and more radiative forcing there.Initial methane-driven warming forces increased soil remineralization and CO 2 input to the atmosphere until subsequent pCO 2 rise and associated CO 2 fertilization turns the tables.Carbon and oxygen isotope excursions of bulk biogenic carbonate in model SSCs are attenuated, distended and delayed relative to carbon and oxygen isotope excursions in the ocean surface layer where the carbonate is formed.Oxygen isotope excursions in surface-layer water compromise the use of oxygen isotope excursions in bulk carbonate to gauge surface-layer warming in the low-latitude ocean zone.However, such an oxygen-isotope-based paleo-thermometer works well in the HL ocean zone.
Methane is the most likely main source of carbon-13 depleted carbon inputs forcing deep-time global warming events: not enough organic carbon is available in the land biosphere to supply carbon needed to explain observed carbon isotope excursions and the very large volcanic carbon input needed to explain the observed isotope excursions greatly exceeds carbon inputs consistent with observed CCD shoaling (Shaffer et al., 2016).Such methane input may be in the form of thermogenic methane (δ 13 C ∼ −40 ‰) and methane hydrate (δ 13 C ∼ −60 ‰).There is an ongoing debate as to whether enough methane hydrate for this job could exist for deep-time, warm pre-event conditions (Buffett and Archer, 2004;Gu et al., 2011).To be specific in the present modeling development paper we chose δ 13 C of our prescribed methane inputs to be −40 ‰.However, more extensive and careful work will be required to resolve this issue.As one step in this direction, some of us are working to develop a methane hydrate module for the DCESS model.In combination with the comprehensive methane cycling implemented in the present paper and in the presence of a large methane hydrate reservoir in present ocean sediments, such a module could also be used for improved projections of future, long term warming with the DCESS model, building upon earlier work (Shaffer et al., 2009).It is also possible that different deep-time warming events had different methane sources or the two methane types in some combination (e.g., Freiling et al., 2016).With the present implementation of methane cycling in the DCESS model we are now in the position to systematically zero in on the amounts, types, timescales and locations of methane injections driving specific deep-time global warming events.

Figure 1 .
Figure 1.Modules and interconnections for the Danish Center for Earth System Science (DCESS) Earth system model.

Figure 2 .
Figure 2. Model configuration.The model consists of one hemisphere divided by 52 • latitude into 360 • wide, low-mid-latitude (LL) and high-latitude (HL) sectors.Values for global reservoirs, transports and fluxes are twice the hemispheric values.The model ocean is 270 • wide and extends from the Equator to 70 • latitude.The ocean covers 70.5 % of the model surface.LL and HL ocean sectors have the area ratio 84 : 16.Sea ice and snow coverage are diagnosed from idealized, meridional profiles of atmospheric temperature.

Figure 3 .
Figure3.Atmosphere, ocean and ocean sediment configurations.Temperatures of low-mid (LL) and high (HL) latitude atmosphere sectors derive from energy balance forced by yearly-mean insolation (SW), outgoing longwave radiation (LW), meridional transports of sensible (SH) and latent (LH) heat, and air-sea heat exchange (H).In combination with the simple sea ice and snow (hatched bars) parameterizations, the model includes the ice/snow albedo feedback and the insulating effect of sea ice.Meridional variations in atmospheric temperature are represented by a second order Legendre polynomial in sine of latitude, adjusted to fit sector mean temperatures.Temperatures and gradients at the 52 • dividing latitude from this meridional profile are used to calculate sensible heat, latent heat and water vapor transports.The atmosphere exchanges heat and gases (G) with the ice-free part of the ocean.Gases are well-mixed in the atmosphere and their concentrations follow from internal sinks, exchanges with ocean surface layers and the land biosphere and input/outputs from weathering and lithosphere outgassing.Both LL and HL ocean sectors are continuously stratified with 100 m vertical resolution to maximum depths of 5500 m and each sector has a bottom topography based on real ocean depth distribution.Model ocean circulation and mixing are indicated by the black arrows and associated circulation/mixing type labels (see main text).Model ocean biogeochemical cycling is indicated by the colored arrows (acting in both sectors but for simplicity only shown for LL sector) and associated matter (particulate organic, POC; biogenic carbonate, PIC; and inert mineral, NCM).Each of the 110 model sediment segments has an area set by the real ocean depth distribution and consists of a 10 cm thick, bioturbated layer divided into seven sublayers.Sediment biogeochemical cycling is indicated by colored arrows and associated matter.Exchange of dissolved species at the sediment-water interface is also indicated.The sedimentation rate down out of the bioturbated layer (w s ) is calculated from mass conservation.

Figure 4 .
Figure 4. Atmospheric lifetime of methane as a function of atmospheric methane concentration.See text for a discussion of the model results and the functional fit to them.

Figure 5 .
Figure 5. Radiative forcing of carbon dioxide, methane and nitrous oxide as functions of their atmospheric concentrations.(a) CO 2 radiative forcing used here (new CO 2 ; Byrne and Goldblatt, 2014a) compared to CO 2 forcing in the original DCESS model (old CO 2 ; Shaffer et al., 2008).CO 2 -N 2 O overlap (Byrne and Goldblatt, 2014a) and N 2 O radiative forcing calculated for constant, pre-industrial N 2 O partial pressure of 270 ppb.(b) CH 4 radiative forcing used here (new CH 4 ; Byrne and Goldblatt, 2014a, b) compared to CH 4 forcing in the original DCESS model (old CH 4 ).CH 4 -N 2 O overlap (Byrne and Goldblatt, 2014a) and N 2 O radiative forcing calculated for constant N 2 O partial pressure as in (a).
Figure7shows selected model results for the first 20 000 years of simulations like that of Fig.6but with methane input amounts of 2000, 4000 and 6000 Gt C. Both pCO 2 and especially pCH 4 increase more than linearly in response to the linear increase in methane input (Fig.7a, c and d).For pCO 2 this is due to positive feedbacks from the land biosphere and ocean.Initial warming from the methane injection leads to greater soil remineralization, land carbon inventory decrease and input of CO 2 to the atmosphere (Fig.7f).Increasing methane oxidation in the ocean leads to denitrification that draws down ocean new production leading to increased ocean CO 2 outgassing that was already enhanced in response to decreased ocean CO 2 solubility from warming (Fig.7e).The non-linear increase in pCH 4 is due mainly to the increase in CH 4 atmospheric lifetime with pCH 4 (Fig.4).On the other hand, warming increases less than linearly (Fig.7b) as expected from the radiative forcing dependencies shown in Fig.5and from our neglect here for simplicity of climate sensitivity increase with warming(Shaffer et al., 2016).Net carbon input from volcanism/weathering minus burial follows process sequences as described for Fig.6above (Fig.7g).The CCD shoals as more methane is oxidized to CO 2 in the atmosphere and ocean leading to more acidic conditions that drive calcite dissolution in the sediment (Fig.7h).However, CCD shoaling increase is damped for large methane dissolution in the ocean since this leads to anoxic conditions, methane oxidation with sulfate (Reaction R7) and associated alkalinity inputs that oppose dissolution.This also has had the effect of damping the atmospheric pCO 2 increase.In response to warming, oxygen isotopes of the LL ocean surface water are enriched by enhanced evaporation associated with a stronger atmospheric hydrological cycle (Fig.7i).LL ocean surface-layer salinity (not shown) increases in response to warming for the Figure7.Results for 20 000-year simulations for different methane input amounts over a timescale of 3000 years and with half the input dissolved in the ocean and half escaping as gas to the atmosphere.(a) methane input rate, (b) global mean atmospheric temperature, (c) atmospheric partial pressure of carbon dioxide, (d) atmospheric partial pressure of methane, (e) total ocean new production, (f) total land biosphere biomass, (g) volcanic/weathering inputs minus burial outputs (as black line in Fig.6g), (h) LL CCD depth (10 % CaCO 3 dry weight in sediment), (i) oxygen isotope excursion of LL ocean surface-layer water, and (j) atmospheric carbon isotope excursion (the methane input has δ 13 C = −40 ‰).
Figure8.Results for 20 000-year simulations for different methane input timescales for an input of 4000 Gt C with half the input dissolved in the ocean and half escaping as gas to the atmosphere.Properties plotted in (a-j) as in Fig.7.

Figure 10 .
Figure 10.Comparison of pathways and feedbacks for all methane injection into the ocean or into the atmosphere.

Table 1 .
Pre-event steady-state solution properties compared to pre-industrial properties.CCD is the calcite compensation depth. .