Oceanic and atmospheric methane cycling in the cGENIE Earth system model

The methane 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 methane can be produced and destroyed, and with the flexibility to be able to explore not only relatively recent perturbations to Earth's methane cycle but also to probe methane cycling and associated climate impacts under the reducing conditions characteristic of most of Earth history and likely widespread on other Earth-like planets. Here, we present an expansion of the ocean-atmosphere methane cycle in the intermediate-complexity Earth system model cGENIE, including parameterized atmospheric photochemistry and schemes for microbial methanogenesis, aerobic methanotrophy, and anaerobic oxidation of methane. 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 methane levels and marine dissolved methane 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 time-dependent behavior of the methane cycle resulting from transient carbon injection into the atmosphere, and present model ensembles that examine the effects of oceanic chemistry and the thermodynamics of microbial metabolism on steady-state atmospheric methane abundance.


Introduction
The global biogeochemical cycle of methane (CH4) 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 (Valentine, 2011;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 CH4 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). As the most abundant hydrocarbon in Earth's atmosphere CH4 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 CH4 classified as a critical near-term climate forcing (Myhre et al., 2013), but has also resulted in dramatic impacts during certain periods of Earth history. For example, high steady-state atmospheric CH4 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 CH4 inventory have been invoked as drivers of extreme climatic perturbations throughout Earth 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 CH4 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 CH4 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 CH4 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-or 2-D photochemical models while surface fluxes are specified arbitrarily or are based on a simple 1-box ocean-biosphere model. A range of slightly more complex 'box' model approaches have been applied to simulate transient perturbations to Earth's CH4 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 CH4 degassing from gas hydrate reservoirs have been developed and applied to relatively recent periods of Earth 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 CH4 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 focus 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 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 modelling of marine CH4 cycling by developing and employing a 3-D ocean circulation and climate model (CCSM-3) to simulate the impact of injecting clathrate-derived CH4 into the Arctic ocean. However, microbial consumption of CH4 in the ocean interior was parameterized via an empirical log-linear function that implicitly neglects anaerobic oxidation of methane (AOM) via dissolved sulfate (SO4 2-), which on the modern Earth is an enormously important internal CH4 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 modelling approaches by applying an ocean general circulation model (GCM) to examine the role of aerobic methanotrophy in modulating ocean-atmosphere fluxes of CH4 during Archean time (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 CH4 sink in the oceans even at relatively low dissolved SO4 2levels.
Though this represented an important further step forward in understanding marine CH4 cycling on the early Earth, Olson et al. (2016) employed a simplified parameterization of aerobic CH4 consumption, neglected the thermodynamics of CH4-consuming metabolisms under energylimited conditions, and employed a parameterization of atmospheric O2-O3-CH4 photochemistry that is most readily applicable to only a subset of the atmospheric pO2 values characteristic of Earth history (Daines and Lenton, 2016;Olson et al., 2016). While all of these studies provided new modelling innovations and advances in understanding, important facets of global CH4 cycling, particularly as relevant to the evolution of early Earth, were lacking.
Here, we present a new framework for modeling the ocean-atmosphere biogeochemical CH4 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 CH4 cycle on ocean-bearing planets, with an initial focus on periods of Earth history (or other habitable ocean-bearing planets) that lack a robust terrestrial biosphere. We also aim to provide a numerical modeling foundation from which to further develop a more complete CH4 cycle within the cGENIE framework, including, for example, dynamic CH4 hydrate cycling and the production/consumption of CH4 by terrestrial ecosystems.
The outline of the paper is as follows. In Section 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 CH4 cycle. In Section 3 we describe the major microbial metabolisms involved in the oceanic CH4 cycle and compare our parameterizations to data from modern marine environments. In Section 4 we describe two alternative parameterizations of atmospheric O2-O3-CH4 photochemistry incorporated into the model and compare these to modern/recent observations. In Section 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 Section 7, following a brief summary in Section 6.

Ocean physics and climate model -C-GOLDSTEIN
The ocean physics and climate model in cGENIE is 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 dynamic-thermodynamic 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 freshwater 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 x 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 (PO4): where DOP represents dissolved organic phosphorus, u represents the proportion of photosynthetic production that is initially partitioned into a dissolved organic phase, l represents a decay constant (time -1 ) for dissolved organic matter, and represents photosynthetic nutrient uptake following Doney et al. (2006): . (3) Rates of photosynthesis are regulated by terms describing the impact of available light (FI), nutrient abundance (FN), temperature (FT), and fractional sea ice coverage in each grid cell (fice).
Rates of photosynthetic nutrient uptake are further scaled to ambient dissolved PO4 ([PO4]) according to an optimal uptake timescale (tbio).
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' zc in Doney et al. (2006) of 75 m. The terms kI and kP represent half-saturation constants for light and dissolved phosphate, respectively. In addition, the effect of temperature on nutrient uptake is parameterized according to: , where kT0 and keT 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 two change in rate with a temperature change of 10ºC (e.g., a Q10 response of ~2.0). Lastly, the (3), not present in the default parameterization of Ridgwell et al. (2007), allows for biological productivity to scale more directly with available PO4 when dissolved PO4 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 Fz POM is the particulate organic matter flux at a given depth (and zh is the base of the photic zone), z is depth, ri POM and li POM 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 (O2, NO3 -, and SO4 2-) are consumed according to decreasing free energy yield (Froelich et al., 1979), with consumption rates (Ri) 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 (NO3 -). The total consumption of settling POM within each ocean layer is governed by the predetermined remineralization profiles (Equation 7). The Ri 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 (l) determines the total fraction of DOM degraded (and Equations 8-10 again determine how the consumption of electron acceptors is partitioned). The ki terms represent half-saturation constants for each metabolism, ki 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 H2 is assumed to be quantitatively converted to CH4 by hydrogenotrophic methanogens. We thus ignore the scenario in which some fraction of H2 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 CH4 (e.g., where ki and ki i terms are as described above (Table 1). We disable nitrate (NO3) 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) (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 non-trivial anaerobic carbon remineralization within sinking particles even in the presence of relatively high [O2] in the ocean water column (Fig. 1c).
While the model tracks the carbon isotope composition of oceanic and atmospheric CH4 (d 13 C, reported in per mil 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 CH4 during methanogenesis of -35‰ by default (Table 2), which will tend to produce microbial CH4 with a d 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/anaerobic methanotrophy, air-sea gas exchange of CH4, or photochemical breakdown of CH4 in the atmosphere. It does, however, include a comprehensive 13 C scheme associated with ocean-atmosphere cycling of CO2 (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 (RAER) with a mixed kinetic-thermodynamic formulation (Jin & Bethke, 2005;Regnier et al., 2011), in which CH4 oxidation kinetics are controlled by substrate availability, thermodynamic energy yield, and temperature: .
A rate constant for aerobic methanotrophy (y -1 ) is defined as kAER, while Fi terms denote kinetic (k) and thermodynamic (t) factors as defined below and a temperature (T) factor as given in Eq.
The kinetic factor (Fk) for aerobic methanotrophy is controlled by substrate availability according to: , where brackets denote concentration and the k term denotes a half-saturation constant with respect to O2. We employ a hybrid parameterization in which kinetics are first-order with respect to CH4 but also scaled by a Michaelis-Menten-type term for O2. This formulation is based on the rationale that half-saturation constants for CH4 are typically similar to (or greater than) the dissolved CH4 levels attained in anoxic water column environments (Regnier et al., 2011) but is also meant to allow for rapid CH4 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 ∆Gr denotes the Gibbs free energy of reaction under in-situ conditions, ∆GBQ represents the minimum energy required to sustain ATP synthesis (Hoehler et al., 2001;Hoehler, 2004;Jin and Bethke, 2007), c 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, ∆Gr 0 represents the Gibbs free energy of the reaction under standard conditions, and gi values represent activity coefficients. Note that we assume an H2O activity of unity.

Anaerobic oxidation of methane (AOM)
The oxidation of methane can also be coupled to electron acceptors other than O2, including nitrate (NO3 -), sulfate (SO4 2-), 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 SO4 2-: . 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 CH4 at the sulfatemethane transition zone (SMTZ) represents an extremely large sink flux of CH4 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 CH4 ( 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 influence of substrate availability, thermodynamic drive, and temperature: where kAOM is a rate constant for anaerobic methane oxidation (y -1 ), while Fi terms denote kinetic (k) and thermodynamic (t) factors as defined below and a temperature (T) factor as given in Eq.
The kinetics of anaerobic methane oxidation are specified according to: where brackets denote concentration and the k term denotes a half-saturation constant with respect to SO4 2-. We employ a hybrid parameterization in which kinetics are first-order with respect to CH4 but are also scaled by a Michaelis-Menten-type term for SO4 2for reasons discussed above.
The effect of thermodynamic energy yield on anaerobic methane oxidation is specified as follows: As above, ∆Gr denotes the Gibbs free energy of reaction under in-situ conditions, ∆GBQ is the minimum energy required to sustain ATP synthesis (the 'biological quantum'), c 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 ∆Gr 0 again represents the Gibbs free energy of the net AOM reaction given above under standard conditions, and gi values represent activity coefficients. Again, we assume an H2O 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 Supplementary Data), 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 CH4 fluxes from the ocean that are consistent with recent observational constraints (Fig.   2c). It is important to note, however, that these values are not extensively tuned and could be adjusted depending on the application. For example, transient CH4 release experiments could employ rate constants that are scaled upward to reflect transient ('bloom') elevations in microbial community CH4 consumption as observed in field studies (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/mixed culture experiments with aerobic methanotrophs Conrad, 1992, 1993;Hanson and Hanson, 1996;Dunfield and Conrad, 2000; van Bodegom et al., 2001), and 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  , 2006). We assume a default biological quantum (∆GBQ) 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/lower salinity than those characteristic of Earth's modern oceans.

Air-sea gas exchange
Ocean-atmosphere fluxes of CH4 (Jgas) 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), [CH4]cell denotes the ambient dissolved CH4 concentration in a given surface ocean grid cell, [CH4]sat represents the dissolved CH4 concentration at saturation with a given atmospheric pCH4, temperature, and salinity, and kgas represents a gas transfer velocity. Solubility is based on a Bunsen solubility coefficient (b) corrected for ambient temperature (T) and salinity (S) according to: , 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 CH4 follows by default that for other gases accounted for in BIOGEM, such as O2 and CO2, as described in (Ridgwell et al., 2007)

Parameterized O2-O3-CH4 photochemistry
Once degassed to the atmosphere, CH4 becomes involved in a complex series of photochemical reactions initiated by hydroxyl radical (OH) attack on CH4 (Kasting et al., 1983;Prather, 1996;Pavlov et al., 2000;Schmidt and Shindell, 2003 where Mi terms represent the atmospheric inventories of O2 and CH4, respectively, and keff denotes an effective rate constant (Tmol -1 y -1 ) that is itself a complicated function of atmospheric O2, CH4, and CO2 (Claire et al., 2006). At each timestep, the distribution of chemical species (e.g., other than temperature and humidity) in the atmosphere is homogenized (Ridgwell et al., 2007) and keff is estimated based on the resulting instantaneous mean partial pressures of O2 and CH4 according to a bivariate fit to a large suite of 1-D atmospheric photochemical models. These photochemical model results (Claire, personal communication) are derived following Claire et al. (2006). Briefly, values for keff are computed by a 1-D model of atmospheric photochemistry assuming a range of fixed surface mixing ratios of O2 and CH4 and a constant atmospheric CO2 of 10 -2 bar. We then fit a fifth-order polynomial surface to these keff values as a function of atmospheric pO2 and pCH4 (Fig. 3).
Our default parameterization of O2-O3-CH4 chemistry (C06) is fit over a pO2 range of 10 -14 to 10 -1 bar, a pCH4 range of 10 -6 to 2 x 10 -3 bar, and a constant high background pCO2 of 10 -2 bar (Claire et al., 2006). We thus truncate the atmospheric lifetime of CH4 at a lower bound of 7.6 years in our default parametrization, and provide an alternative parameterization of photochemical CH4 destruction at roughly modern pO2 and pCO2 (SS03) derived from the results of Schmidt and Shindell (2003) for use in more geologically recent, high-O2 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 CH4 in the atmosphere, including atmospheric pCO2, the atmospheric profile of H2O, 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 CH4 to the atmosphere (balanced by stoichiometric consumption of CO2 and release of O2), and allow the oceanic and atmospheric CH4 cycle to spin up for 20 kyr. We then compare steady-state atmospheric pCH4 as a function of terrestrial CH4 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 O2 and CH4 inventories. Nevertheless, both the default scheme and the alternative parameterization for recent geologic history (and analogous planetary environments) with high-pO2/low-pCO2 atmospheres accurately reproduce atmospheric pCH4 values given estimated glacial, preindustrial, and modern terrestrial CH4 fluxes (Fig. 4C), and both display the predicted saturation of CH4 sinks at elevated atmospheric CH4 observed in more complex photochemical models. We note, however, the alternative parameterization tends to yield slightly higher atmospheric pCH4 at surface fluxes greater than ~50 Tmol y -1 (Fig. 4C). (In the remainder of the manuscript, we employ the default (C06) parameterization for atmospheric O2-O3-CH4 chemistry and do not discuss the simple high-pO2/low-pCO2 alternative further.)

High-pO2 ('modern') steady state
We explore a roughly modern steady state with appropriate continental geography and simulated overturning circulation (as in Cao et al., 2009)  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 [O2] 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 between roughly 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 [O2] 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 SO4 2-([SO4 2-]) 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). Benthic [SO4 2-] distributions show some similarity to those of [O2] (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 salinity-normalized (not shown). In the ocean interior, [SO4 2-] is largely spatially invariant with a value of approximately 28 mmol kg -1 (Fig. 5h).
Dissolved CH4 concentrations ([CH4]) 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 [CH4] 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 CH4 cycling (Jayakumar et al., 2001). Within the ocean interior, dissolved CH4 can accumulate in the water column in excess of ~100 nmol kg -1 in association with relatively low-O2 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  and is particularly vigorou s 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 CH4 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, 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 CH4 produced via microbial methanogenesis is consumed via AOM, either near the seafloor or within the ocean interior.

Low-pO2 ('ancient') steady state
Next, we explore a low-pO2 steady state, similar to the Proterozoic Earth (Reinhard et al., 2017) but played out in a modern continental configuration and overturning circulation, by initializing for our low-pO2 simulation. Dissolved O2 concentrations are now extremely heterogeneous throughout the surface ocean, ranging over an order of magnitude from less than 1 µmol kg -1 to over 10 µmol kg -1 , with concentrations that are regionally well in excess of air saturation at the prescribed pO2 of 2.1 x 10 -4 atm (Fig. 7a). Previous studies have shown that these features are not unexpected at very low atmospheric pO2 (Olson et al., 2013;Reinhard et al., 2016). We note, however, that the distribution and maximum [O2] in our low-pO2 simulation are both somewhat different from those presented in Olson et al. (2013) and Reinhard et al. (2016). We attribute this primarily to the different 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 PO4 3than 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 [O2] dynamics at low atmospheric pO2 (Olson et al., 2013;Reinhard et al., 2016), our regional [O2] patterns still generally track the localized balance between photosynthetic O2 release and consumption through respiration and reaction with inorganic reductants, rather than temperature-dependent solubility patterns (Fig. 5a). Within the ocean interior, O2 is consumed within the upper few hundred meters and is completely absent in benthic settings (Fig. 7d, g).
In our low-pO2 simulations we initialize the ocean with a globally uniform [SO4 2-] of 280 µmol kg -1 , under the premise that marine SO4 2inventory should scale positively with atmospheric pO2.
With this much lower initial SO4 2inventory (i.e., 10 2 times less than the modern ocean), steady state [SO4 2-] distributions are significantly more heterogeneous than in the modern, high-pO2 case (Fig. 7). Ocean [SO4 2-] 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-pO2 simulations SO4 2serves as the principal oxidant for organic matter remineralization in the ocean interior, with the result that its distribution effectively mirrors that of [O2] in the modern case in both spatial texture and overall magnitude (compare Fig. 7e, h with Fig. 5d, g). Dissolved SO4 2in this simulation never drops to zero, a consequence of our initial 280 µmol kg -1 concentration of SO4 2representing the oxidative potential of 560 µmol kg -1 of O2, some 3 times higher than the mean [O2] value in the modern ocean interior (~170 µmol kg -1 ).
Dissolved CH4 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) (Fig. 7f), which results from the fact that in the low-pO2 case SO4 2again serves as the principal oxidant of methane. Concentrations of CH4 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 CH4 inventory increases dramatically in the low-pO2 case relative to the modern simulation, from ~4.5 Tmol CH4 to ~1900 Tmol CH4. Figure 8 shows the major metabolic fluxes within the ocean's microbial CH4 cycle for our 'ancient' configuration. Column-integrated rates of microbial methanogenesis are greater than in the high-pO2 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 O2 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 the high-pO2 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-pO2 case (Fig. 8c, f). Once again, AOM dominates the consumption of CH4 produced in the ocean interior and acts as an extremely effective throttle on CH4 fluxes to the atmosphere. Despite a significant increase in overall oceanic CH4 burden relative to our high-pO2 simulation (see above and Fig. 7i), atmospheric pCH4 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 CH4 cycle, we perform a simple carbon injection experiment in which 3,000 PgC are injected directly into the atmosphere either as CH4 or as CO2, starting from our modern steady state. The injection is spread over 1,000 years, with an instantaneous initiation and termination of carbon input to the atmosphere. The magnitude and duration of this carbon injection, corresponding to 3 PgC y -1 , is meant to roughly mimic the upper end of estimates for the Paleocene-Eocene Thermal Maximum, a transient global warming event at ~56 Ma hypothesized to have been driven by emissions of CO2 and/or CH4 (Kirtland Turner, 2018). This flux is much lower than the current anthropogenic carbon input of ~10 PgC y -1 . For simplicity, and because we focus on only the first 3,000 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 CH4, there is an immediate and significant increase in atmospheric pCH4 to values greater than 10 ppmv, followed by a gradual increase to a maximum of ~12 ppmv throughout the duration of the CH4 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 pCO2 (Fig. 9b). This increase in atmospheric pCH4 and pCO2 leads to an increase in global average surface air temperature (SAT) of ~7ºC (Fig. 9d), 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 pCO2 and drop in ocean pH are nearly identical if we instead inject the carbon as CO2 rather than CH4. (Fig. 9b, c). However, when carbon is injected as CH4, there is an additional transient increase in global surface air temperature of ~2ºC and roughly 0.5ºC of additional whole ocean warming for the same carbon input and duration (Fig. 9f). This results from the fact that mole-for-mole, CH4 is a much more powerful greenhouse gas than is CO2, and oxidation of CH4 to CO2 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 CH4 rather than CO2, providing an enhancement of warming, especially during the duration of carbon input. This warming enhancement should be considered in past events during which CH4 release is suspected as a key driver of warming. For instance, additional warming due to CH4 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 pCH4 on the early Earth
Using our low-pO2 steady state as a benchmark case (Section 5.2), we briefly explore the sensitivity of atmospheric pCH4 to a subset of model variables. All model ensembles are initially configured with globally homogeneous marine SO4 2and CH4 inventories and a background geologic CH4 flux of 3 Tmol y -1 , and are spun up for 20 kyr with a fixed pO2 and pCO2. We report atmospheric pCH4 from the final model year. Our purpose here is not to be exhaustive or to elucidate any particular period of Earth history, but to demonstrate some of the major factors controlling the atmospheric abundance of CH4 on a low-oxygen Earth-like planet. We present results from individual sensitivity ensembles from our benchmark low-pO2 case over the following parameter ranges: (1) atmospheric pO2 between 10 -4 to 10 -1 times the present atmospheric level (PAL), equivalent to roughly 2 x 10 -5 and 2 x 10 -2 atm, respectively; (2) initial marine SO4 2inventories corresponding to globally uniform seawater concentrations between 0 and 1,000 µmol kg -1 ; and (3) biological energy quanta (BEQ) for anaerobic methane oxidation between 5 and 30 kJ mol -1 .
Results for our low-pO2 sensitivity ensembles are shown in Figure 10. We find a similar sensitivity of atmospheric pCH4 to atmospheric pO2 to that observed by . In particular, atmospheric CH4 abundance initially increases as atmospheric pO2 drops below modern values to roughly 2-3% PAL, after which decreasing pO2 causes pCH4 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 H2O by ozone (O3) decreases at low atmospheric pO2 (Pavlov et al., 2003;Claire et al., 2006;Goldblatt et al., 2006). However, peak atmospheric pCH4 is significantly reduced in our models relative to those of Olson et al. (2016).
For example, at an 'optimal' atmospheric pO2 of ~2.5% PAL Olson et al. (2016) predict a steady state atmospheric pCH4 of ~35 ppmv, while we predict a value of ~10 ppmv (Fig. 10a). This difference can be attributed to our updated O2-O3-CH4 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 CH4 as a significant climate regulator at steady state during most of the Proterozoic Eon (between ~2.5 and 0.5 Ga).
Atmospheric CH4 abundance is also strongly sensitive to the marine SO4 2inventory (Fig. 10b).
The scaling we observe between initial SO4 2inventory and steady state atmospheric pCH4 is very similar to that reported by Olson et al. (2016), with a sharp drop in the marine CH4 inventory and atmospheric CH4 abundance as marine SO4 2drops below ~100 µmol kg -1 (Fig. 10b). The implication is that for most of Earth history anaerobic oxidation of CH4 in the ocean interior has served as an important inhibitor of CH4 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 SO4 2concentrations may instead have been on the order of ~1-10 µmol kg -1 (Crowe et al., 2014), while atmospheric pO2 would also have been much lower than the values examined here (Pavlov & Kasting, 2002). The impact of the ocean biosphere and redox chemistry on atmospheric pCH4 and Earth's climate system may thus have been much more important prior to ~2.5 billion years ago.
Interestingly, atmospheric CH4 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 pCH4 from ~7 ppmv 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 pO2 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 CH4 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 CH4 cycle and seawater [SO4 2-] is relatively low.

Discussion and Conclusions
The global biogeochemical cycling of CH4 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 CH4 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 CH4 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 global CH4 cycle in Earth's past and future and ultimately to constrain CH4 cycling dynamics on Earth-like 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 CH4 production. We have also incorporated revised schemes for microbial CH4 consumption that include both kinetic and thermodynamic constraints, and have updated the We suggest that ongoing and future development work should focus on: (1) more rigorous tuning of organic carbon remineralization and CH4 production/consumption schemes based on data fields from the modern ocean; (2) development and implementation of a more flexible parameterization of atmospheric photochemistry that allows the roles of atmospheric temperature structure, water vapor abundance, and atmospheric pCO2 to be explored; (3) coupling of deep ocean chemistry with a description of marine methane hydrates and associated sedimentary CH4 cycling; and (4) developing a representation of the production/consumption of CH4 by terrestrial ecosystems.

Model 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     Table 2. Shown in (c) are globally integrated diffusive fluxes of CH4 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., 2019) Figure 9. Response to a 3,000 PgC release directly to the atmosphere spread over 1,000 years, assuming carbon is injected as either CH4 or CO2. Atmospheric pCH4 (a), pCO2 (b), mean surface ocean pH (c), mean surface air temperature (SAT; d), and mean ocean temperature (MOT; e) are shown for a CH4 injection (grey) and a CO2 injection (black). Panel (f) shows the difference in SAT and MOT between the CH4 and CO2 injection scenarios (∆Ti = TCH4,i -TCO2,i) through time. c.