Oceanic and atmospheric methane cycling in the cGENIE Earth system model – release v0.9.14

The methane (CH4) cycle is a key component of the Earth system that links planetary climate, biological metabolism, and the global biogeochemical cycles of carbon, oxygen, sulfur, and hydrogen. However, currently lacking is a numerical model capable of simulating a diversity of environments in the ocean, where CH4 can be produced and destroyed, and with the flexibility to be able to explore not only relatively recent perturbations to Earth’s CH4 cycle but also to probe CH4 cycling and associated climate impacts under the very low-O2 conditions characteristic of most of Earth’s history and likely widespread on other Earth-like planets. Here, we present a refinement and expansion of the ocean–atmosphere CH4 cycle in the intermediate-complexity Earth system model cGENIE, including parameterized atmospheric O2–O3–CH4 photochemistry and schemes for microbial methanogenesis, aerobic methanotrophy, and anaerobic oxidation of methane (AOM). We describe the model framework, compare model parameterizations against modern observations, and illustrate the flexibility of the model through a series of example simulations. Though we make no attempt to rigorously tune default model parameters, we find that simulated atmospheric CH4 levels and marine dissolved CH4 distributions are generally in good agreement with empirical constraints for the modern and recent Earth. Finally, we illustrate the model’s utility in understanding the timedependent behavior of the CH4 cycle resulting from transient carbon injection into the atmosphere, and we present model ensembles that examine the effects of atmospheric pO2, oceanic dissolved SO2− 4 , and the thermodynamics of microbial metabolism on steady-state atmospheric CH4 abundance. Future model developments will address the sources and sinks of CH4 associated with the terrestrial biosphere and marine CH4 gas hydrates, both of which will be essential for comprehensive treatment of Earth’s CH4 cycle during geologically recent time periods.


Introduction
The global biogeochemical cycle of methane (CH 4 ) is central to the evolution and climatic stability of the Earth system. Methane provides an important substrate for microbial metabolism, particularly in energy-limited microbial ecosystems in the deep subsurface Chapelle et al., 1995) and in anoxic marine and lacustrine sediments (Lovley et al., 1982;Hoehler et al., 2001). Indeed, the microbial production and consumption of CH 4 are amongst the oldest metabolisms on Earth, with an isotopic record of bacterial methane cycling stretching back nearly 3.5 billion years (Ueno et al., 2006;Hinrichs, 2002;Hayes, 1994).
Published by Copernicus Publications on behalf of the European Geosciences Union. 5688 C. T. Reinhard et al.: Oceanic and atmospheric methane cycling in the cGENIE Earth system model As the most abundant hydrocarbon in Earth's atmosphere CH 4 also has a significant influence on atmospheric photochemistry (Thompson and Cicerone, 1986), and because it absorbs in a window region of Earth's longwave emission spectrum it is an important greenhouse gas. This has important implications over the coming centuries, with atmospheric CH 4 classified as a critical near-term climate forcing (Myhre et al., 2013), but has also resulted in dramatic impacts during certain periods of Earth's history. For example, high steady-state atmospheric CH 4 has been invoked as an important component of Earth's early energy budget, potentially helping to offset a dim early Sun (Sagan and Mullen, 1972;Pavlov et al., 2000;Haqq-Misra et al., 2008), while time-dependent changes to the atmospheric CH 4 inventory have been invoked as drivers of extreme climatic perturbations throughout Earth's history (Dickens et al., 1997;Dickens, 2003;Bjerrum and Canfield, 2011;Zeebe, 2013;Schrag et al., 2002). Because it is cycled largely through biological processes on the modern (and ancient) Earth and is spectrally active, atmospheric CH 4 has also been suggested as a remotely detectable biosignature that could be applied to planets beyond our solar system (Hitchcock and Lovelock, 1967;Sagan et al., 1993;Krissansen-Totton et al., 2018).
A number of low-order Earth system models incorporating a basic CH 4 cycle have been developed, particularly with a view to addressing relatively "deep time" geological questions. These include explorations of long-term changes to the chemistry of Earth's atmosphere (Claire et al., 2006;Catling et al., 2007;Bartdorff et al., 2008;Beerling et al., 2009), potential climate impacts at steady state (Kasting et al., 2001;Ozaki et al., 2018), and transient impacts of CH 4 degassing on climate (Schrag et al., 2002;Bjerrum and Canfield, 2011). In some cases these models explicitly couple surface fluxes to a model of atmospheric photochemistry (Lamarque et al., 2006;Ozaki et al., 2018;Kasting et al., 2001), but in general atmospheric chemistry is parameterized based on offline 1-D or 2-D photochemical models, while surface fluxes are specified arbitrarily or are based on a simple one-box oceanbiosphere model. A range of slightly more complex "box" model approaches have been applied to simulate transient perturbations to Earth's CH 4 cycle and attendant climate impacts on timescales ranging from ∼ 10 5 years (Dickens et al., 1997;Dickens, 2003) to ∼ 10 8 years (Daines and Lenton, 2016). In addition, offline and/or highly parameterized approaches toward simulating the impact of transient CH 4 degassing from gas hydrate reservoirs have been developed and applied to relatively recent periods of Earth's history (Archer and Buffett, 2005;Lunt et al., 2011) or projected future changes Hunter et al., 2013). However, the most sophisticated and mechanistic models of global CH 4 cycling currently available tend to focus on terrestrial (soil or wetland) sources and sinks (Ridgwell et al., 1999;Walter and Heimann, 2000;Wania et al., 2010;Konijnendijk et al., 2011;Melton et al., 2013) or on explicitly modeling atmospheric photochemistry .
Much less work has been done to develop ocean biogeochemistry models that are both equipped to deal with the wide range of boundary conditions characteristic of Earth's history and are computationally tractable when running large model ensembles and/or on long (approaching ∼ 10 6 year) timescales, as well as being able to simulate the (3-D) redox structure of the ocean, allowing for localized zones of production and oxidation (which provides more accurate estimates of emission to the atmosphere). For instance, Elliot et al. (2011) advanced modeling of marine CH 4 cycling by developing and employing a 3-D ocean circulation and climate model (CCSM-3) to simulate the impact of injecting clathrate-derived CH 4 into the Arctic Ocean. However, microbial consumption of CH 4 in the ocean interior was parameterized via an empirical log-linear function that implicitly neglects anaerobic oxidation of methane (AOM) via dissolved sulfate (SO 2− 4 ), which on the modern Earth is an enormously important internal CH 4 sink within Earth's oceans (Egger et al., 2018). Their simulations did not explore atmospheric chemistry. Similarly, Daines and Lenton (2016) also innovated over traditional box modeling approaches by applying an ocean general circulation model (GCM) to examine the role of aerobic methanotrophy in modulating ocean-atmosphere fluxes of CH 4 during the Archean (prior to ∼ 2.5 billion years ago, Ga). However, this analysis likewise did not include AOM, and the GCM results were not coupled to atmospheric chemistry. In contrast, Olson et al. (2016) included AOM in a 3-D ocean biogeochemistry model coupled to an atmospheric chemistry routine and found that AOM represents a critical internal CH 4 sink in the oceans even at relatively low dissolved SO 2− 4 levels. Though this represented an important further step forward in understanding marine CH 4 cycling on the early Earth, Olson et al. (2016) employed a simplified parameterization of aerobic CH 4 consumption, neglected the thermodynamics of CH 4 -consuming metabolisms under energylimited conditions, and employed a parameterization of atmospheric O 2 -O 3 -CH 4 photochemistry that is most readily applicable to only a subset of the atmospheric pO 2 values characteristic of Earth's history (Daines and Lenton, 2016;Olson et al., 2016). While all of these studies provided new modeling innovations and advances in understanding, important facets of global CH 4 cycling, particularly as relevant to the evolution of the early Earth, were lacking.
Here, we present a new framework for modeling the ocean-atmosphere biogeochemical CH 4 cycle in the "muffin" release of the cGENIE Earth system model. Our goal is to make further progress in the development of a flexible intermediate-complexity model suitable for simulating the global biogeochemical CH 4 cycle on ocean-bearing planets, with an initial focus on periods of Earth's history (or other habitable planets) that lack a robust terrestrial biosphere. We also aim to provide a numerical modeling foundation from which to further develop a more complete CH 4 cycle within the cGENIE framework, including, for example, dynamic C. T. Reinhard et al.: Oceanic and atmospheric methane cycling in the cGENIE Earth system model 5689 CH 4 hydrate cycling and the production and consumption of CH 4 by terrestrial ecosystems.
The outline of the paper is as follows. In Sect. 2 we briefly describe the GENIE-cGENIE Earth system model, with a particular eye toward the features that are most relevant for the biological carbon pump and the oceanic CH 4 cycle. In Sect. 3 we describe the major microbial metabolisms involved in the oceanic CH 4 cycle and compare our parameterizations to data from modern marine environments. In Sect. 4 we describe two alternative parameterizations of atmospheric O 2 -O 3 -CH 4 photochemistry incorporated into the model and compare these to modern and more recent observations. In Sect. 5 we present results from a series of idealized simulations meant to illustrate the flexibility of the model and some potential applications. The availability of the model code, plus configuration files for all experiments described in the paper, is provided in Sect. 7, following a brief summary in Sect. 6.
2 The GENIE-cGENIE Earth system model

Ocean physics and climate model -C-GOLDSTEIN
The ocean physics and climate model in cGENIE are comprised of a reduced physics (frictional geostrophic) 3-D ocean circulation model coupled to both a 2-D energy-moisture balance model (EMBM) and a dynamicthermodynamic sea ice model (Edwards and Marsh, 2005;Marsh et al., 2011). The ocean model transports heat, salinity, and biogeochemical tracers using a scheme of parameterized isoneutral diffusion and eddy-induced advection (Griffies, 1998;Edwards and Marsh, 2005;Marsh et al., 2011), exchanges heat and moisture with the atmosphere, sea ice, and land, and is forced at the ocean surface by the input of zonal and meridional wind stress according to a specified wind field. The 2-D atmospheric energy-moisture balance model (EMBM) considers the heat and moisture balance for the atmospheric boundary layer using air temperature and specific humidity as prognostic tracers. Heat and moisture are mixed horizontally throughout the atmosphere and exchange heat and moisture with the ocean and land surfaces with precipitation occurring above a given relative humidity threshold. The sea ice model tracks the horizontal transport of sea ice and the exchange of heat and fresh water with the ocean and atmosphere using ice thickness, areal fraction, and concentration as prognostic variables. Full descriptions of the model and coupling procedure can be found in Edwards and Marsh (2005) and, more recently, in Marsh et al. (2011). As implemented here, the ocean model is configured as a 36 × 36 equal-area grid (uniform in longitude and uniform in the sine of latitude) with 16 logarithmically spaced depth levels and seasonal surface forcing from the EMBM.

Ocean biological pump -BIOGEM
The biogeochemical model component -BIOGEM -regulates air-sea gas exchange as well as the transformation and partitioning of biogeochemical tracers within the ocean, as described in Ridgwell et al. (2007). By default, the biological pump is driven by parameterized uptake of nutrients in the surface ocean, with this flux converted stoichiometrically to biomass that is then partitioned into either dissolved or particulate organic matter for downstream transport, sinking, and remineralization. Dissolved organic matter is transported by the ocean model and decays with a specified time constant, while particulate organic matter is immediately exported out of the surface ocean and partitioned into two fractions of differing lability. In the ocean interior, particulate organic matter is remineralized instantaneously throughout the water column following an exponential decay function with a specified remineralization length scale.
In the simulations discussed below, photosynthetic nutrient uptake in surface ocean grid cells is controlled by a single limiting nutrient, dissolved phosphate (PO 4 ): where DOP represents dissolved organic phosphorus, υ represents the proportion of photosynthetic production that is initially partitioned into a dissolved organic phase, λ represents a decay constant (time −1 ) for dissolved organic matter, and represents photosynthetic nutrient uptake following Doney et al. (2006): Rates of photosynthesis are regulated by terms describing the impact of available light (F I ), nutrient abundance (F N ), temperature (F T ), and fractional sea ice coverage in each grid cell (f ice ). Rates of photosynthetic nutrient uptake are further scaled to ambient dissolved PO 4 ([PO 4 ]) according to an optimal uptake timescale (τ bio ). Note that this parameterization differs from that in Ridgwell et al. (2007). Specifically, the impacts of light and nutrient availability are both described via Michaelis-Menten terms: where shortwave irradiance I is averaged over the entire mixed layer and is assumed to decay exponentially from the sea surface with a length scale of 20 m. It is assumed that nutrient uptake and photosynthetic production only occur in surface grid cells of cGENIE (e.g., the upper 80 m), which is similar to the "compensation depth" z c in Doney et al. (2006) of 75 m. The terms κ I and κ P represent half-saturation constants for light and dissolved phosphate, respectively. In addition, the effect of temperature on nutrient uptake is parameterized according to where k T0 and k eT denote pre-exponential and exponential scaling constants and T represents absolute in situ temperature. The scaling constants are chosen to give approximately a factor of 2 change in rate with a temperature change of 10 • C (e.g., a Q 10 response of ∼ 2.0). Lastly, the final term in Eq.
(3), not present in the default parameterization of Ridgwell et al. (2007), allows biological productivity to scale more directly with available PO 4 when dissolved PO 4 concentrations are elevated relative to those of the modern oceans.
Particulate organic matter (POM) is immediately exported out of the surface ocean without lateral advection and is instantaneously remineralized throughout the water column according to an exponential function of depth: where F POM z is the particulate organic matter flux at a given depth (and z h is the base of the photic zone), z is depth, and r POM i and l POM i refer to the relative partitioning into each organic matter lability fraction i and the e-folding depth of that fraction, respectively. The simulations presented here employ two organic matter fractions, a "labile" fraction (94.5 %) with an e-folding depth of ∼ 590 m and an effectively inaccessible fraction (5.5 %) with an e-folding depth of 10 6 m ( Table 1).
We employ a revised scheme for organic matter remineralization in the ocean interior, following that commonly used in models of organic matter remineralization within marine and lacustrine sediments (Rabouille and Gaillard, 1991;Van Cappellen et al., 1993;Boudreau, 1996a, b). Respiratory electron acceptors (O 2 , NO − 3 , and SO 2− 4 ) are consumed according to decreasing free energy yield (Froelich et al., 1979), with consumption rates (R i ) scaled to both electron acceptor abundance and the inhibitory impact of electron acceptors with higher intrinsic free energy yield: with the exception that in the biogeochemical configuration used here we do not consider nitrate (NO − 3 ). The total consumption of settling POM within each ocean layer is governed by the predetermined remineralization profiles (Eq. 7). The R i terms denote the relative fraction of this organic matter consumption that is performed by each respiratory process. We specify a closed system with no net organic matter burial in marine sediments (see below), and hence the POM flux to the sediment surface is assumed to be completely degraded, with the same partitioning amongst electron acceptors carried out according to local bottom water chemistry. For DOM, the assumed lifetime (λ) determines the total fraction of DOM degraded (and Eqs. 8-10 again determine how the consumption of electron acceptors is partitioned). The κ i terms represent half-saturation constants for each metabolism, the κ i i terms give inhibition constants acting on less energetic downstream respiratory processes, and brackets denote concentration. Default parameter values used here are shown in Table 1.

Microbial methanogenesis
Methanogenesis represents the terminal step in our remineralization scheme and follows the overall stoichiometry: This can be taken to implicitly include fermentation of organic matter to acetate followed by acetoclastic methanogenesis: or the fermentation of organic matter to acetate followed by anaerobic acetate oxidation and hydrogenotrophic methanogenesis: both of which have the same overall net stoichiometry provided that H 2 is assumed to be quantitatively converted to CH 4 by hydrogenotrophic methanogens. We thus ignore the scenario in which some fraction of H 2 is converted directly to biomass by hydrogenotrophic methanogens acting as primary producers (Ozaki et al., 2018).
Because we specify a closed system with no net organic matter burial in marine sediments, all organic matter not remineralized by more energetic respiratory metabolisms is converted into CH 4 (e.g.,  where the κ i and κ i i terms are as described above (Table 1). We disable nitrate (NO 3 ) as a tracer in the simulations presented here such that anaerobic remineralization of organic matter is partitioned entirely between sulfate reduction and methanogenesis ( Fig. 1). Using our default parameter values (Table 1), aerobic respiration dominates organic matter remineralization at [O 2 ] values significantly above 1 µmol kg −1 (Fig. 1a), while anaerobic remineralization is dominated by methanogenesis at [SO 2− 4 ] values significantly below 1 mmol kg −1 (Fig. 1b). An important outcome of the revised "inhibition" scheme is that metabolic pathways with differing intrinsic free energy yields can coexist, which more accurately reflects field observations from a range of natural settings (Curtis, 2003;Bethke et al., 2008;Kuivila et al., 1989;Jakobsen and Postma, 1999). In particular, it allows us to roughly capture the impact of oxidant gradients within sinking marine aggregates (Bianchi et al., 2018), which can facilitate nontrivial anaerobic carbon remineralization within sinking particles even in the presence of relatively high [O 2 ] in the ocean water column (Fig. 1c).
While the model tracks the carbon isotope composition of oceanic and atmospheric CH 4 (δ 13 C, reported in per mill notation relative to the Pee Dee Belemnite, PDB), the only significant isotope effect we include here is that attendant to acetoclastic methanogenesis. We specify a constant isotope fractionation between organic carbon and CH 4 during methanogenesis of −35 ‰ by default (Table 2), which will tend to produce microbial CH 4 with a δ 13 C composition of roughly −60 ‰ when combined with the default isotope fractionation associated with photosynthetic carbon fixation in the surface ocean (e.g., Kirtland Turner and Ridgwell, 2016). The model does not currently include any potential isotope effects associated with aerobic and anaerobic methanotrophy, air-sea gas exchange of CH 4 , or photochemical breakdown of CH 4 in the atmosphere. It does, however, include a comprehensive 13 C scheme associated with ocean-atmosphere cycling of CO 2 (Kirtland Turner and Ridgwell, 2016;Ridgwell, 2001).

Aerobic methanotrophy
Microbial aerobic methanotrophy proceeds according to This reaction is highly favorable energetically, with a free energy yield under standard conditions of ∼ 850 kJ per mole of methane consumed (Table 2). We represent rates of aerobic methanotrophy (R AER ) with a mixed kinetic-thermodynamic formulation Bethke, 2005, 2007;Regnier et al., 2011), in which CH 4 oxidation kinetics are controlled by substrate availability, thermodynamic energy yield, and temperature: A rate constant for aerobic methanotrophy (yr −1 ) is defined as k AER , while F i terms denote kinetic (k) and thermodynamic (t) factors as defined below and a temperature (T ) factor as given in Eq. (6) above. . Shown in (c) are our estimated anaerobic remineralization fractions (grey curve) compared to estimates from a particle biogeochemical model applied to oxygen minimum zones (OMZs) in the eastern tropical South Pacific (ETSP) and Mauritanian upwelling (Bianchi et al., 2018). The kinetic factor (F k ) for aerobic methanotrophy is controlled by substrate availability according to where brackets denote concentration and the κ term denotes a half-saturation constant with respect to O 2 . We employ a hybrid parameterization in which kinetics are first-order with respect to CH 4 but also scaled by a Michaelis-Menten-type term for O 2 . This formulation is based on the rationale that half-saturation constants for CH 4 are typically similar to (or greater than) the dissolved CH 4 levels attained in anoxic water column environments (Regnier et al., 2011) but is also meant to allow for rapid CH 4 consumption under "bloom" conditions with an appropriately scaled rate constant (see below). The effect of thermodynamic energy yield on aerobic methanotrophy is given by where G r denotes the Gibbs free energy of reaction under in situ conditions, G BQ represents the minimum energy required to sustain ATP synthesis (Hoehler et al., 2001;Hoehler, 2004;Jin and Bethke, 2007), χ is the stoichiometric number of the reaction (e.g., the number of times the ratedetermining step occurs in the overall process), and R and T represent the gas constant and absolute in situ temperature, respectively. The available free energy is estimated according to where, in addition to the terms defined above, G 0 r represents the Gibbs free energy of the reaction under standard conditions, and γ i values represent activity coefficients. Note that we assume an H 2 O activity of unity.

Anaerobic oxidation of methane (AOM)
The oxidation of methane can also be coupled to electron acceptors other than O 2 , including nitrate (NO − 3 ), sulfate (SO 2− 4 ), and oxide phases of iron (Fe) and manganese (Mn) (Reeburgh, 1976;Martens and Berner, 1977;Hoehler et al., 1994;Hinrichs et al., 1999;Orphan et al., 2001;Sivan et al., 2011;Haroon et al., 2013;Egger et al., 2015). Because it is by far the most abundant of these oxidants on the modern Earth and has likely been the most abundant throughout Earth's history, we focus on anaerobic oxidation of methane (AOM) at the expense of SO 2− 4 : This process is currently thought to be performed most often through a syntrophic association between Archaea and sulfate-reducing bacteria (Boetius et al., 2000), though the mechanics controlling the exchange of reducing equivalents within the syntrophy remain to be fully elucidated (Milucka et al., 2012;McGlynn et al., 2015). In any case, consumption of CH 4 at the sulfate-methane transition zone (SMTZ) represents an extremely large sink flux of CH 4 in modern marine sediments (Regnier et al., 2011;Egger et al., 2018). Anaerobic methanotrophy is much less energetically favorable under standard conditions, with a free energy yield of ∼ 30 kJ per mole of CH 4 ( Table 2). As a result, the influence of thermodynamics on rates of AOM is potentially much stronger than it will tend to be in the case of aerobic methanotrophy. As above, rates of AOM are controlled by the combined where k AOM is a rate constant for anaerobic methane oxidation (yr −1 ), while F i terms denote kinetic (k) and thermodynamic (t) factors as defined below and a temperature (T ) factor as given in Eq. (6) above. The kinetics of anaerobic methane oxidation are specified according to where brackets denote concentration and the κ term denotes a half-saturation constant with respect to SO 2− 4 . We employ a hybrid parameterization in which kinetics are first-order with respect to CH 4 but are also scaled by a Michaelis-Mententype term for SO 2− 4 for reasons discussed above. The effect of thermodynamic energy yield on anaerobic methane oxidation is specified as follows: As above, G r denotes the Gibbs free energy of reaction under in situ conditions, G BQ is the minimum energy required to sustain ATP synthesis (the "biological quantum"), χ is the stoichiometric number of the reaction, and R and T represent the gas constant and absolute in situ temperature, respectively. The available free energy for AOM under in situ conditions is estimated according to where G 0 r again represents the Gibbs free energy of the net AOM reaction given above under standard conditions, and γ i values represent activity coefficients. Again, we assume an H 2 O activity of unity.

Default parameters for aerobic and anaerobic methanotrophy
We choose default rate constants according to a dataset of compiled rates of aerobic and anaerobic methanotrophy in oxygenated and anoxic marine water column environments (see the Supplement), after correction to in situ temperature (Fig. 2a, b). Our default values for both rate constants are on the low end of the observational dataset but are very roughly tuned to yield steady-state diffusive CH 4 fluxes from the ocean that are consistent with recent observational constraints (Fig. 2c). It is important to note, however, that these  Table 2. Shown in (c) are globally integrated diffusive fluxes of CH 4 from the ocean for a range of rate constants for aerobic methanotrophy, including our default simulation. The bar to the right of (c) shows the median (black bar) and 90 % credible interval (grey shading) for estimates of the modern oceanic diffusive flux from Weber et al.  (Kessler et al., 2011;Crespo-Medina et al., 2014). Default values for other kinetic parameters (Table 2) are chosen to be broadly consistent with field measurements and pure or mixed culture experiments with aerobic methanotrophs Conrad, 1992, 1993;Hanson and Hanson, 1996;Dunfield and Conrad, 2000;van Bodegom et al., 2001), as well as to remain roughly consistent with previous work for comparative purposes (e.g., Olson et al., 2016), though the parameters have not been formally tuned and we explore model sensitivity below. Thermodynamic energy yields of each reaction under standard conditions are calculated based on the standard molal thermodynamic properties given in Regnier et al. (2011). Stoichiometric numbers are assumed identical for both metabolisms, with default values of 1.0 (Jin and Bethke, 2005;Dale et al., 2006). We assume a default biological quantum ( G BQ ) of 15 kJ mol −1 for both aerobic and anaerobic methanotrophy, though these can be expected to vary somewhat as a function of metabolism and environmental conditions (Schink, 1997;Hoehler, 2004;Dale et al., 2008). These can be varied independently for aerobic and anaerobic methanotrophy in the model, and we explore model sensitivity to this parameter below. Lastly, for simplicity and to minimize computational expense we assume constant activity coefficients for each species throughout all ocean grid cells (Table 2). For some applications it may ultimately be important to add a scheme for estimating activity coefficients according to ambient salinity and ion chemistry, for example estimating methane fluxes in planetary scenarios with very different major ion chemistry or much higher or lower salinity than those characteristic of Earth's modern oceans.
4 Atmospheric methane cycling 4.1 Air-sea gas exchange Ocean-atmosphere fluxes of CH 4 (J gas ) are governed by temperature-and salinity-dependent solubility and surface wind speed above a given grid cell: where A denotes the area available for gas exchange (e.g., the area of ice-free surface ocean), [CH 4 ] cell denotes the ambient dissolved CH 4 concentration in a given surface ocean grid cell, [CH 4 ] sat represents the dissolved CH 4 concentration at saturation with a given atmospheric pCH 4 , temperature, and salinity, and k gas represents a gas transfer velocity. Solubility is based on a Bunsen solubility coefficient (β) corrected for ambient temperature (T ) and salinity (S) according to ln β =a 1 + a 2 (100/T ) + a 3 ln(T /300) Note that the Henry's law constant K 0 is related to the Bunsen solubility coefficient by K 0 = β/ρV + , where ρ is density and V + is the molar volume of the gas at standard temperature and pressure (STP). Gas transfer velocity (k gas ) is calculated based on the surface wind speed (u) and a Schmidt number (Sc) corrected for temperature assuming a constant salinity of 35 ‰: where k is a dimensionless gas transfer coefficient, u is surface wind speed, and Sc is the temperature-corrected Schmidt number according to All default constants and coefficients for the gas exchange scheme are given in Table 3. Overall, the scheme for air-sea gas exchange of CH 4 follows by default that for other gases accounted for in BIOGEM, such as O 2 and CO 2 , as described in Ridgwell et al. (2007)

Parameterized O 2 -O 3 -CH 4 photochemistry
Once degassed to the atmosphere, CH 4 becomes involved in a complex series of photochemical reactions initiated by hydroxyl radical (OH) attack on CH 4 (Kasting et al., 1983;Prather, 1996;Pavlov et al., 2000;Schmidt and Shindell, 2003). Following Claire et al. (2006) and Goldblatt et where M i terms represent the atmospheric inventories of O 2 and CH 4 , respectively, and k eff denotes an effective rate constant (Tmol −1 yr −1 ) that is itself a complicated function of atmospheric O 2 , CH 4 , and CO 2 (Claire et al., 2006). We note that in this parameterization, O 3 abundance is not calculated explicitly, but rather the photochemical destruction rate of CH 4 in the atmosphere is controlled by the combined atmospheric chemistry implicitly embedded within k eff (Goldblatt et al., 2006;Claire et al., 2006). At each time step, the distribution of chemical species (e.g., other than temperature and humidity) in the atmosphere is homogenized (Ridgwell et al., 2007) and k eff is estimated based on the resulting instantaneous mean partial pressures of O 2 and CH 4 according to a bivariate fit to a large suite of 1-D atmospheric photochemical models. These photochemical model results (Mark Claire, personal communication, 2016) are derived following Claire et al. (2006). Briefly, values for k eff are computed by a 1-D model of atmospheric photochemistry assuming a range of fixed surface mixing ratios of O 2 and CH 4 and a constant atmospheric CO 2 of 10 −2 bar. We then fit a fifth-order polynomial surface to these k eff values as a function of atmospheric pO 2 and pCH 4 (Fig. 3).
Our default parameterization of O 2 -O 3 -CH 4 chemistry (C06) is fit over a pO 2 range of 10 −14 to 10 −1 bar, a pCH 4 range of 10 −6 to 2 × 10 −3 bar, and a constant high background pCO 2 of 10 −2 bar (Claire et al., 2006). We thus truncate the atmospheric lifetime of CH 4 at a lower bound of 7.6 years in our default parametrization and provide an alternative parameterization of photochemical CH 4 destruction at roughly modern pO 2 and pCO 2 (SS03) derived from the results of Schmidt and Shindell (2003) for use in more geologically recent, high-O 2 atmospheres (Reinhard et al., 2017) (Fig. 4a). Although this parameterized photochemistry scheme should represent an improvement in accuracy relative to that implemented in Olson et al. (2016) (see Daines and Lenton, 2016), it is important to point out that a range of factors that might be expected to impact the photochemical destruction rates of CH 4 in the atmosphere, including atmospheric pCO 2 , the atmospheric profile of H 2 O, and spectral energy distribution (SED), have not yet been rigorously assessed. Ongoing model developments in ATCHEM are aimed at implementing a more flexible and inclusive photochemical parameterization that will allow for robust use across a wider range of atmospheric compositions and photochemical environments.
As a basic test of our photochemical parameterization, we impose a terrestrial (wetland) flux of CH 4 to the atmosphere (balanced by stoichiometric consumption of CO 2 and release of O 2 ) and allow the oceanic and atmospheric CH 4 cycle to spin up for 20 kyr. We then compare steady-state atmospheric pCH 4 as a function of terrestrial CH 4 flux to estimates for the last glacial, preindustrial, and modern periods. Our default parameterization is relatively simple and spans a very wide range in atmospheric O 2 and CH 4 inventories. Nevertheless, both the default scheme and the alternative parameterization for recent geologic history (and analogous planetary environments), with high-pO 2 ,low-pCO 2 atmospheres, accurately reproduce atmospheric pCH 4 values given estimated glacial, preindustrial, and modern terrestrial CH 4 fluxes (Fig. 4c), and both display the predicted saturation of CH 4 sinks at elevated atmospheric CH 4 observed in more complex photochemical models. We note, however, that the alternative parameterization tends to yield slightly higher atmospheric pCH 4 at surface fluxes greater than ∼ 50 Tmol yr −1 (Fig. 4c). (In the remainder of the paper, we employ the default (C06) parameterization for atmospheric O 2 -O 3 -CH 4 chemistry and do not discuss the simple high-pO 2 −low-pCO 2 alternative further.) Figure 4. Comparison of steady-state atmospheric pCH 4 as a function of terrestrial CH 4 flux with modern and more recent estimates. Shown in (a) is an exponential fit to the 2-D photochemistry model of Schmidt and Shindell (2003) (SS03), with individual model runs shown as black crosses. Shown in (b) is a plane through the bivariate fit shown in Fig. 3 (grey curve) compared with the ensemble of 1-D atmospheric photochemical models at pO 2 = 0.1 atm (black crosses; see text). Shown in (c) are steady-state atmospheric CH 4 values as a function of imposed terrestrial CH 4 flux in our "modern" configuration (circles) compared to estimates for the glacial, preindustrial, and modern CH 4 cycles (Kirschke et al., 2013;Bock et al., 2017;Paudel et al., 2016).

High-pO 2 ("modern") steady state
We explore a roughly modern steady state with appropriate continental geography and simulated overturning circulation (as in Cao et al., 2009) and initialize the atmosphere with pO 2 , pCO 2 , and pCH 4 of [v/v] 20.95 %, 278 ppm, and 700 ppb, respectively, and with globally uniform oceanic concentrations of SO 2− 4 (28 mmol kg −1 ) and CH 4 (1 nmol kg −1 ). We fix globally averaged solar insolation at the modern value (1368 W m −2 ) with seasonally variable forcing as a function of latitude and set radiative forcing for CO 2 and CH 4 equivalent to preindustrial values in order to isolate the effects of biogeochemistry on steady-state tracer distributions. The model is then spun up for 20 kyr with atmospheric pO 2 and pCO 2 (and δ 13 C of atmospheric CO 2 ) restored to preindustrial values at every time step and with an imposed wetland flux of CH 4 to the atmosphere of 20 Tmol yr −1 that has a δ 13 C value of −60 ‰. Atmospheric pCH 4 and all oceanic tracers are allowed to evolve freely.
Surface, benthic, and ocean interior distributions of dissolved oxygen (O 2 ), sulfate (SO 2− 4 ), and methane (CH 4 ) are shown in Fig. 5 for our roughly modern simulation. Dissolved O 2 ([O 2 ]) approaches air saturation throughout the surface ocean, with a distribution that is largely uniform zonally and with concentrations that increase with latitude as a result of increased solubility at lower temperature near the poles (Fig. 5a). Benthic [O 2 ] shows patterns similar to those expected for the modern Earth, with relatively high values in the well-ventilated deep North Atlantic, low values in the deep North Pacific and Indian oceans, and a gradient roughly between air saturation near regions of deep convection in the high-latitude Atlantic and much lower values in the tropical and northern Pacific (Fig. 5d). Distributions of [O 2 ] in the ocean interior are similar to those of the modern Earth ( Fig. 5g) with oxygen minimum zones (OMZs) at intermediate depths underlying highly productive surface waters, particularly in association with coastal upwelling at low latitudes.
Concentrations of dissolved SO 2− 4 ([SO 2− 4 ]) are largely invariant throughout the ocean, consistent with its expected conservative behavior in the modern ocean as one of the most abundant negative ions in seawater (Fig. 5). Slightly higher concentrations in both surface and benthic fields are seen in association with outflow from the Mediterranean and are driven by evaporative concentration (Fig. 5b) (Fig. 5e), though again the differences are very small relative to the overall prescribed initial tracer inventory of 28 mmol kg −1 and disappear almost entirely when salinitynormalized (not shown). In the ocean interior, [SO 2− 4 ] is largely spatially invariant with a value of approximately 28 mmol kg −1 (Fig. 5h).
Dissolved CH 4 concentrations ([CH 4 ]) in the surface and shallow subsurface ocean are much more variable but typically on the order of ∼ 1-2 nmol kg −1 , with slightly elevated concentrations just below the surface, both of which are consistent with observations from the modern ocean (Reeburgh, 2007;Scranton and Brewer, 1978). The benthic [CH 4 ] distribution shows locally elevated values up to ∼ 300-400 nmol kg −1 in shallow regions of the tropical and northern Pacific and the Indian oceans (Fig. 5f), which is also broadly consistent with observations from shallow marine environments with active benthic CH 4 cycling (Jayakumar et al., 2001). Within the ocean interior, dissolved CH 4 can accumulate in the water column in excess of ∼ 100 nmol kg −1 in association with relatively low-O 2 conditions at intermediate depths, with zonally averaged values as high as ∼ 70 nmol kg −1 but more typically in the range of ∼ 20- 40 nmol kg −1 (Fig. 5i). These concentrations are comparable to those observed locally in low-O 2 regions of the modern ocean (Sansone et al., 2001;Chronopoulou et al., 2017;Thamdrup et al., 2019). The major metabolic fluxes within the ocean's microbial CH 4 cycle for our "modern" configuration are shown in Fig. 6. Methanogenesis is focused in regions characterized by relatively low [O 2 ] and is particularly vigorous in the eastern tropical Pacific, the North Pacific, and the Indian Ocean (Fig. 6a). The highest zonally averaged rates of methanogenesis are observed in northern tropical and subtropical latitudes and are focused at a depth of ∼ 1 km (Fig. 6d). Rates of microbial CH 4 consumption are generally spatially coupled to rates of methanogenesis, both in a column-integrated sense (Fig. 6b, c) and in the zonal average (Fig. 6e, f). This is particularly true for AOM, the rates of which are highest within the core of elevated methanogenesis rates observed in the northern subtropics. Zonally averaged AOM rates of ∼ 10-15 nmol kg −1 d −1 compare well with field measurements of AOM within oceanic OMZs (Thamdrup et al., 2019). In general, the bulk of CH 4 produced via microbial methanogenesis is consumed via AOM, either near the seafloor or within the ocean interior.

Low-pO 2 ("ancient") steady state
Next, we explore a low-pO 2 steady state, similar to the Proterozoic Earth (Reinhard et al., 2017) but played out in a modern continental configuration and overturning circulation, by initializing the atmosphere with pO 2 , pCO 2 , and pCH 4 of [v/v] 2.1 × 10 −4 atm (equivalent to a value 10 −3 times the present atmospheric level, PAL), 278 ppm, and 500 ppm, respectively, as well as globally uniform oceanic concentrations of SO 2− 4 (280 µmol kg −1 ) and CH 4 (50 µmol kg −1 ). We again fix globally averaged solar insolation at the modern value (1368 W m −2 ), with seasonally variable forcing as a function of latitude, and set radiative forcing for CO 2 and CH 4 equivalent to the modern preindustrial state in order to isolate the effects of biogeochemistry on steady-state tracer distributions. The model is then spun up for 20 kyr with atmospheric pO 2 and pCO 2 (and δ 13 C of atmospheric CO 2 ) restored to the initial values specified above at every time step, with an imposed "geologic" flux of CH 4 to the atmosphere of 3 Tmol yr −1 at a δ 13 C value of −60 ‰. Atmospheric pCH 4 and all oceanic tracers are allowed to evolve freely.
Surface, benthic, and ocean interior distributions of [O 2 ], [SO 2− 4 ], and [CH 4 ] are shown in Fig. 7 for our low-pO 2 simulation. Dissolved O 2 concentrations are now extremely heterogeneous throughout the surface ocean, ranging over an order of magnitude from less than 1 to over 10 µmol kg −1 , with concentrations that are regionally well in excess of air saturation at the prescribed pO 2 of 2.1×10 −4 atm (Fig. 7a). Previous studies have shown that these features are not unexpected at very low atmospheric pO 2 (Olson et al., 2013;. We note, however, that the distribution and maximum [O 2 ] in our low-pO 2 simulation are both somewhat different from those presented in Olson et al. (2013) and Reinhard et al. (2016). We attribute this primarily to the differ- Figure 6. Major biological fluxes in the marine methane cycle for our "modern" configuration. Panels show column integrated (a-c) and zonally averaged (d-f) rates of methanogenesis, aerobic methanotrophy, and anaerobic methane oxidation (AOM) in the ocean interior. ent parameterizations of primary production in the surface ocean. In the biogeochemical configuration of cGENIE we adopt here, we allow rates of photosynthesis to scale more directly with available PO 3− 4 than is the case in these previous studies (Eq. 3), which allows for higher rates of oxygen production in regions of deep mixing and relatively intense organic matter recycling below the photic zone (Fig. 7a). In any case, as in previous examinations of surface [O 2 ] dynamics at low atmospheric pO 2 (Olson et al., 2013;Reinhard et al., 2016), our regional [O 2 ] patterns still generally track the localized balance between photosynthetic O 2 release and consumption through respiration and reaction with inorganic reductants, rather than temperature-dependent solubility patterns (Fig. 5a). Within the ocean interior, O 2 is consumed within the upper few hundred meters and is completely absent in benthic settings (Fig. 7d, g).
In our low-pO 2 simulations we initialize the ocean with a globally uniform [SO 2− 4 ] of 280 µmol kg −1 , under the premise that the marine SO 2− 4 inventory should scale positively with atmospheric pO 2 . With this much lower initial SO 2− 4 inventory (i.e., 10 2 times less than the modern ocean), steady-state [SO 2− 4 ] distributions are significantly more heterogeneous than in the modern, high-pO 2 case (Fig. 7). Ocean  ] is approximately homogeneous spatially in surface waters, even with a significantly reduced seawater inventory (Fig. 7a), but is strongly variable within the ocean interior (Fig. 7e, h). Indeed, in our low-pO 2 simulations SO 2− 4 serves as the principal oxidant for organic matter remineralization in the ocean interior, with the result that its distribution effectively mirrors that of [O 2 ] in the modern case in both spatial texture and overall magnitude (compare Fig. 7e, h with Fig. 5d, g). Dissolved SO 2− 4 in this simulation never drops to zero, a consequence of our initial 280 µmol kg −1 concentration of SO 2− 4 representing the oxidative potential of 560 µmol kg −1 of O 2 , some 3 times higher than the mean [O 2 ] value in the modern ocean interior (∼ 170 µmol kg −1 ).
Dissolved CH 4 concentrations in the surface and shallow subsurface ocean are variable but much higher than in our modern simulations, typically on the order of ∼ 1-2 µmol kg −1 (Fig. 7c). The benthic [CH 4 ] distribution shows concentrations up to ∼ 8 µmol kg −1 , with concentrations in excess of 1 µmol kg −1 pervasively distributed across the seafloor. In general, the benthic [CH 4 ] distribution inversely mirrors that of [SO 2− 4 ] (Fig. 7f), which results from the fact that in the low-pO 2 case SO 2− 4 again serves as the principal oxidant of methane. Concentrations of CH 4 in the ocean interior can approach ∼ 10 µmol kg −1 but in the zonal average are typically less than 5 µmol kg −1 (Fig. 7i). Overall, the oceanic CH 4 inventory increases dramatically in the low-pO 2 case relative to the modern simulation, from ∼ 4.5 Tmol CH 4 to ∼ 1900 Tmol CH 4 .
The major metabolic fluxes within the ocean's microbial CH 4 cycle for our "ancient" configuration are shown in Fig. 8. Column-integrated rates of microbial methanogenesis are greater than in the high-pO 2 case by up to a factor of ∼ 10 2 (Fig. 8a), with methanogenesis also showing a much broader areal distribution. Within the ocean interior, rates of methanogenesis are most elevated in the upper ∼ 1 km (Fig. 8d) as a consequence of elevated rates of organic carbon remineralization combined with a virtual absence of dissolved O 2 beneath the upper ∼ 200 m. Rates of aerobic methanotrophy, which is effectively absent in the ocean interior (Fig. 8e), are elevated relative to those observed in the high-pO 2 simulation by less than an order of magnitude and are concentrated in the tropical surface ocean near the equatorial divergence (Fig. 8b). In contrast, AOM is strongly coupled spatially to microbial methanogenesis, with rates that are often well over ∼ 10 2 times higher than those observed in the high-pO 2 case (Fig. 8c, f). Once again, AOM dominates the consumption of CH 4 produced in the ocean interior and is extremely effective at reducing CH 4 fluxes to the atmosphere. Despite a significant increase in overall oceanic CH 4 burden relative to our high-pO 2 simulation (see above and Fig. 7i), atmospheric pCH 4 increases only modestly from ∼ 0.8 ppm to 6 ppm [v/v], equivalent to an additional radiative forcing of only ∼ 2 W m −2 , due to efficient microbial consumption in the upper ocean.

Atmospheric carbon injection
To illustrate the capabilities of the model in exploring the time-dependent (perturbation) behavior of the CH 4 cycle, we perform a simple carbon injection experiment in which 3000 PgC is injected directly into the atmosphere either as CH 4 or as CO 2 , starting from our modern steady state. The injection is spread over 1000 years, with an instantaneous initiation and termination of carbon input to the atmosphere. This is meant only to illustrate the time-dependent behavior of the model in the face of an idealized carbon cycle perturbation, rather than to evaluate any particular scenario for explaining previous climate transients in Earth's history. However, the magnitude and duration of this carbon injection, corresponding to 3 PgC yr −1 , are meant to roughly mimic the upper end of estimates for the Paleocene-Eocene Thermal Maximum (PETM), a transient global warming event at ∼ 56 Ma hypothesized to have been driven by emissions of CO 2 and/or CH 4 (Kirtland Turner, 2018). This flux is much lower than the current anthropogenic carbon input of ∼ 10 PgC yr −1 . For simplicity, and because we focus on only the first 3000 years following carbon injection, we treat the ocean-atmosphere system as closed, with the result that all injected carbon ultimately accumulates within the ocean and atmosphere rather than being removed through carbonate compensation and silicate weathering.
Following a carbon release to the atmosphere in the form of CH 4 , there is an immediate and significant increase in atmospheric pCH 4 to values greater than 10 ppmv, followed by a gradual increase to a maximum of ∼ 12 ppmv throughout the duration of the CH 4 input (Fig. 9a). Much of this methane is exchanged with the surface ocean and consumed by aerobic methanotrophy, while some is photochemically oxidized directly in the atmosphere, both of which lead to a significant but delayed increase in atmospheric pCO 2 (Fig. 9b). This increase in atmospheric pCH 4 and pCO 2 leads to an increase in global average surface air temperature (SAT) of ∼ 7 • C (Fig. 9d) and an increase in mean ocean temperature (MOT) of ∼ 2 • C (Fig. 9e), along with significant acidification of the surface ocean (Fig. 9c).
The increase in atmospheric pCO 2 and drop in ocean pH are nearly identical if we instead inject the carbon as CO 2 rather than CH 4 (Fig. 9b, c). However, when carbon is injected as CH 4 , there is an additional transient increase in global surface air temperature of ∼ 2 • C and roughly 0.5 • C Figure 8. Major biological fluxes in the marine methane cycle for our "ancient" configuration. Panels show column-integrated (a-c) and zonally averaged (d-f) rates of methanogenesis, aerobic methanotrophy, and anaerobic methane oxidation (AOM) in the ocean interior. Figure 9. Response to a 3000 PgC release directly to the atmosphere spread over 1000 years, assuming carbon is injected as either CH 4 or CO 2 . Atmospheric pCH 4 (a), pCO 2 (b), mean surface ocean pH (c), mean surface air temperature (SAT; d), and mean ocean temperature (MOT; e) are shown for a CH 4 injection (grey) and a CO 2 injection (black). Panel (f) shows the difference in SAT and MOT between the CH 4 and CO 2 injection scenarios ( T i = T CH 4 ,i -T CO 2 ,i ) through time.
of additional whole-ocean warming for the same carbon input and duration (Fig. 9f). This results from the fact that, mole for mole, CH 4 is a much more powerful greenhouse gas than is CO 2 , and oxidation of CH 4 to CO 2 is not instantaneous during the carbon release interval. Combined, these factors result in a disequilibrium situation in which a proportion of carbon released to the atmosphere remains in the form of CH 4 rather than CO 2 , providing an enhancement of warming, especially during the duration of carbon input. This warming enhancement should be considered in past events during which CH 4 release is suspected as a key driver of warming. For instance, additional warming due to CH 4 forcing may help explain the apparent discrepancy between the amount of warming reconstructed by proxy records and proposed carbon forcing during the PETM (Zeebe et al., 2009).

Atmospheric pCH 4 on the early Earth
Using our low-pO 2 steady state as a benchmark case (Sect. 5.2), we briefly explore the sensitivity of atmospheric pCH 4 to a subset of model variables. All model ensembles are initially configured with globally homogeneous marine SO 2− 4 and CH 4 inventories, as well as a background geologic CH 4 flux of 3 Tmol yr −1 , and are spun up for 20 kyr with a fixed pO 2 and pCO 2 . We report atmospheric pCH 4 from the final model year. Our purpose here is not to be exhaustive or to elucidate any particular period of Earth's history, but to demonstrate some of the major factors controlling the atmospheric abundance of CH 4 on a low-oxygen Earth-like planet. We present results from individual sensitivity ensembles from our benchmark low-pO 2 case over the following parameter ranges: (1) atmospheric pO 2 between 10 −4 and 10 −1 times the present atmospheric level (PAL), equivalent to roughly 2×10 −5 and 2×10 −2 atm, respectively; (2) initial marine SO 2− 4 inventories corresponding to globally uniform seawater concentrations between 0 and 1000 µmol kg −1 ; and (3) biological energy quanta (BEQ) for anaerobic methane oxidation between 5 and 30 kJ mol −1 .
Results for our low-pO 2 sensitivity ensembles are shown in Fig. 10. We find a similar sensitivity of atmospheric pCH 4 to atmospheric pO 2 to that observed by Olson et al. (2016). In particular, atmospheric CH 4 abundance initially increases as atmospheric pO 2 drops below modern values to roughly 2 %-3 % PAL, after which decreasing pO 2 causes pCH 4 to drop. This behavior is well-known from previous 1-D photochemical model analysis and arises principally from increasing production of OH via water vapor photolysis as shielding of H 2 O by ozone (O 3 ) decreases at low atmospheric pO 2 (Pavlov et al., 2003;Claire et al., 2006;Goldblatt et al., 2006). However, peak atmospheric pCH 4 is significantly reduced in our models relative to those of Olson et al. (2016). For example, at an "optimal" atmospheric pO 2 of ∼ 2.5 % PAL, Olson et al. (2016) predict a steady-state atmospheric pCH 4 of ∼ 35 ppmv, while we predict a value of ∼ 10 ppmv (Fig. 10a). This difference can be attributed to our updated O 2 -O 3 -CH 4 photochemistry parameterization together with a significant upward revision in the rate constant for aerobic methanotrophy. Nevertheless, our results strongly reinforce the arguments presented in Olson et al. (2016) and, taken at face value, further marginalize the role of CH 4 as a significant climate regulator at steady state during most of the Proterozoic eon (between ∼ 2.5 and 0.5 Ga).
Atmospheric CH 4 abundance is also strongly sensitive to the marine SO 2− 4 inventory (Fig. 10b). The scaling we observe between the initial SO 2− 4 inventory and steady-state atmospheric pCH 4 is very similar to that reported by , with a sharp drop in the marine CH 4 inventory and atmospheric CH 4 abundance as marine SO 2− 4 drops below ∼ 100 µmol kg −1 (Fig. 10b). The implication is that for most of Earth's history anaerobic oxidation of CH 4 in the ocean interior has served as an important inhibitor of CH 4 fluxes from the ocean biosphere. However, during much of the Archean eon (between 4.0 and 2.5 Ga), sulfur isotope analysis indicates that marine SO 2− 4 concentrations may instead have been on the order of ∼ 1-10 µmol kg −1 (Crowe et al., 2014), while atmospheric pO 2 would also have been much lower than the values examined here (Pavlov and Kasting, 2002). The impact of the ocean biosphere and redox chemistry on atmospheric pCH 4 and Earth's climate system may thus have been much more important prior to ∼ 2.5 billion years ago.
Atmospheric CH 4 is significantly impacted by the value chosen for the biological energy quantum (BEQ). With all other parameters held constant, we observe an increase in steady-state atmospheric pCH 4 from ∼ 7 to ∼ 25 ppmv when increasing the BEQ value from 20 to 30 kJ mol −1 (Fig. 10c). This effect is mediated primarily by the importance of anaerobic methanotrophy when atmospheric pO 2 is low and the ocean interior is pervasively reducing. The standard free energy of AOM is of the same order of magnitude as the BEQ (see above), which elevates the importance of thermodynamic drive in controlling global rates of AOM. We would expect this effect to be much less important when aerobic methanotrophy is the predominant CH 4 -consuming process within the ocean biosphere, as the standard free energy of this metabolism is over an order of magnitude greater than typical BEQ values for microbial metabolism (e.g., Hoehler, 2004). In any case, our results suggest that the role of thermodynamics should be borne in mind in scenarios for which AOM is an important process in the CH 4 cycle and seawater [SO 2− 4 ] is relatively low.

Discussion and conclusions
The global biogeochemical cycling of CH 4 is central to the climate and redox state of planetary surface environments and responds to the internal dynamics of other major biogeochemical cycles across a very wide range of spatial and temporal scales. There is thus strong impetus for the ongoing development of a spectrum of models designed to explore planetary CH 4 cycling, from simple box models to more computationally expensive 3-D models with dynamic and interactive ocean circulation. Our principal goal here is the development of a mechanistically realistic but simple and flexible representation of CH 4 biogeochemical cycling in Earth's ocean-atmosphere system, with the hope that this can be further developed to explore steady-state and time-dependent changes to the global CH 4 cycle in Earth's past and future and ultimately to constrain CH 4 cycling dynamics on Earthlike planets beyond our solar system. To accomplish this, we have refined the organic carbon remineralization scheme in the cGENIE Earth system model to reflect the impact of anaerobic organic matter recycling in sinking aggregates within oxygenated waters and to include the carbon cycling and isotopic effects of microbial CH 4 production. We have also incorporated revised schemes for microbial CH 4 consumption that include both kinetic and thermodynamic constraints, and we have updated the parameterized atmospheric O 2 -O 3 -CH 4 photochemistry to improve accuracy and for use across a wider range of atmospheric pO 2 values than explored in previous work. Simulations of roughly modern (high-O 2 ) and Proterozoic (low-O 2 ) Earth system states demonstrate that the model effectively reproduces the first-order features of the modern oceanatmosphere CH 4 cycle and can be effectively implemented across a wide range of atmospheric O 2 partial pressures and marine SO 2− 4 concentrations. In addition, our results strongly reinforce the conclusions of Olson et al. (2016) for the Proterozoic Earth system, while going beyond this to posit that the thermodynamics of anaerobic CH 4 consumption may have been important in regulating atmospheric CH 4 abundance during the Archean. Finally, our simulation of PETMlike carbon injection demonstrates the importance of explicitly considering CH 4 radiative forcing during transient warming events in Earth history.
We suggest that ongoing and future development work should focus on the following: (1) more rigorous tuning of organic carbon remineralization and CH 4 production and consumption schemes based on data fields from the modern ocean; (2) development and implementation of a more Figure 10. Sensitivity ensembles of our "ancient" configuration compared to the results of Olson et al. (2016). Steady-state atmospheric pCH 4 values as a function of assumed atmospheric pO 2 (a) and initial marine SO 2− 4 inventory (b) are shown for our "ancient" configuration (filled circles; see text) and from Olson et al. (black crosses). Shown below are additional ensembles showing the impact of varying the minimum free energy yield required for microbial methane oxidation (BEQ; c) on atmospheric pCH 4 . All simulations were spun up from cold for 20 kyr, with the results shown from the last model year.
flexible parameterization of atmospheric photochemistry that allows the roles of atmospheric temperature structure, water vapor abundance, and atmospheric pCO 2 to be explored; (3) coupling of deep ocean chemistry with a description of marine methane hydrates and associated sedimentary CH 4 cycling; and (4) developing a representation of the production and consumption of CH 4 by terrestrial ecosystems.
Code availability. A manual describing code installation, basic model configuration, and an extensive series of tutorials is provided. The Latex source of the manual and PDF file can be obtained by cloning (https://github.com/derpycode/muffindoc, last access: 15 November 2020). The user manual contains instructions for obtaining, installing, and testing the code, as well as running experiments. The version of the code used in this paper is tagged as release v0.9.14 and has the DOI https://doi.org/10.5281/zenodo.4002934 (Ridgewell et al., 2020). Configuration files for the specific experiments presented in the paper can be found in cgenie.muffin/genieuserconfigs/MS/reinhardetal.GMD.2020. Details on the different experiments, plus the command line needed to run each, are given in README.txt.
Data availability. The rate data used to construct Fig. 2 are also included in the Supplement and as an accessory dataset with the DOI https://doi.org/10.5281/zenodo.4081700 (Reinhard, 2020).
Author contributions. CTR, SLO, and AR developed new model code. CTR and CP compiled and analyzed empirical data for rates of methanotrophy. CTR performed all model simulations and data analysis. CTR prepared the paper with contributions from all coauthors.
Competing interests. The authors declare that they have no conflict of interest.