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

The methane (CH₄) 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 CH₄ can be produced and destroyed, and with the flexibility to be able to explore not only relatively recent perturbations to Earth's CH₄ cycle but also to probe CH₄ cycling and associated climate impacts under the very low-O₂ 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 CH₄ cycle in the intermediate-complexity Earth system model cGENIE, including parameterized atmospheric O₂–O₃–CH₄ 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 CH₄ levels and marine dissolved CH₄ 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 CH₄ cycle resulting from transient carbon injection into the atmosphere, and we present model ensembles that examine the effects of atmospheric pO₂, oceanic dissolved SO${}_{\mathrm{4}}^{\mathrm{2}-}$, and the thermodynamics of microbial metabolism on steady-state atmospheric CH₄ abundance. Future model developments will address the sources and sinks of CH₄ associated with the terrestrial biosphere and marine CH₄ gas hydrates, both of which will be essential for comprehensive treatment of Earth's CH₄ cycle during geologically recent time periods.

1 Introduction

The global biogeochemical cycle of methane (CH₄) 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 CH₄ 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 CH₄ 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₄ 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₄ 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₄ 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₄ 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₄ 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₄ 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 ocean–biosphere model. A range of slightly more complex “box” model approaches have been applied to simulate transient perturbations to Earth's CH₄ cycle and attendant climate impacts on timescales ranging from ∼ 10⁵ years (Dickens et al., 1997; Dickens, 2003) to ∼ 10⁸ years (Daines and Lenton, 2016). In addition, offline and/or highly parameterized approaches toward simulating the impact of transient CH₄ 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 (Archer et al., 2009; Hunter et al., 2013). However, the most sophisticated and mechanistic models of global CH₄ 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 (Shindell et al., 2013).

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⁶ 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₄ cycling by developing and employing a 3-D ocean circulation and climate model (CCSM-3) to simulate the impact of injecting clathrate-derived CH₄ into the Arctic Ocean. However, microbial consumption of CH₄ in the ocean interior was parameterized via an empirical log-linear function that implicitly neglects anaerobic oxidation of methane (AOM) via dissolved sulfate (SO${}_{\mathrm{4}}^{\mathrm{2}-}$), which on the modern Earth is an enormously important internal CH₄ 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₄ 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₄ sink in the oceans even at relatively low dissolved SO${}_{\mathrm{4}}^{\mathrm{2}-}$ levels. Though this represented an important further step forward in understanding marine CH₄ cycling on the early Earth, Olson et al. (2016) employed a simplified parameterization of aerobic CH₄ consumption, neglected the thermodynamics of CH₄-consuming metabolisms under energy-limited conditions, and employed a parameterization of atmospheric O₂–O₃–CH₄ photochemistry that is most readily applicable to only a subset of the atmospheric pO₂ 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₄ 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₄ 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₄ 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₄ cycle within the cGENIE framework, including, for example, dynamic CH₄ hydrate cycling and the production and consumption of CH₄ 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₄ cycle. In Sect. 3 we describe the major microbial metabolisms involved in the oceanic CH₄ cycle and compare our parameterizations to data from modern marine environments. In Sect. 4 we describe two alternative parameterizations of atmospheric O₂–O₃–CH₄ 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

2.1 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 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 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.

2.2 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₄):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial {\mathrm{PO}}_{\mathrm{4}}}{\partial t}}=-\mathrm{\Gamma }+\mathit{\lambda }\mathrm{DOP},\text{(2)}& {\displaystyle \frac{\partial \mathrm{DOP}}{\partial t}}=\mathit{\upsilon }\mathrm{\Gamma }-\mathit{\lambda }\mathrm{DOP},\end{array}$$

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⁻¹) for dissolved organic matter, and Γ represents photosynthetic nutrient uptake following Doney et al. (2006):

$$\begin{array}{}\text{(3)}& \mathrm{\Gamma }=F_{\mathrm{I}}\cdot F_{\mathrm{N}}\cdot F_{\mathrm{T}}\cdot (\mathrm{1}-f_{\mathrm{ice}})\cdot {\displaystyle \frac{\left[{\mathrm{PO}}_{\mathrm{4}}\right]}{{\mathit{\tau }}_{\mathrm{bio}}}}.\end{array}$$

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₄ ([PO₄]) 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:

$$\begin{array}{}\text{(4)}& {\displaystyle }{F}_{\mathrm{I}}={\displaystyle \frac{I}{I+{\mathit{\kappa }}_{\mathrm{I}}}},\text{(5)}& {\displaystyle }{F}_{\mathrm{N}}={\displaystyle \frac{\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]}{{\mathit{\kappa }}_{\mathrm{P}}+\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]}},\end{array}$$

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

$$\begin{array}{}\text{(6)}& F_{\mathrm{T}}=k_{\mathrm{T}\mathrm{0}}\cdot \exp \left[{\displaystyle \frac{T}{k_{\mathrm{eT}}}}\right],\end{array}$$

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₁₀ 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₄ when dissolved PO₄ 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:

$$\begin{array}{}\text{(7)}& F_z^{\mathrm{POM}}=F_{z=z_h}^{\mathrm{POM}}\cdot \left(\sum _i{r}_i^{\mathrm{POM}}\cdot \exp \left({\displaystyle \frac{z_h-z}{l_i^{\mathrm{POM}}}}\right)\right),\end{array}$$

where $F_z^{\mathrm{POM}}$ 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_i^{\mathrm{POM}}$ and $l_i^{\mathrm{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⁶ m (Table 1).

Table 1 Default parameters for organic matter production and water column remineralization.

¹ Ridgwell et al. (2007); ² Meyer et al. (2016); ³ Doney et al. (2006); ⁴ Olson et al. (2016).

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₂, NO${}_{\mathrm{3}}^{-}$, and SO${}_{\mathrm{4}}^{\mathrm{2}-}$) 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:

$$\begin{array}{}\text{(8)}& {\displaystyle }{R}_{{\mathrm{O}}_{\mathrm{2}}}={\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{{\mathit{\kappa }}_{{\mathrm{O}}_{\mathrm{2}}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}},\text{(9)}& {\displaystyle }{R}_{{\mathrm{NO}}_{\mathrm{3}}}={\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{{\mathit{\kappa }}_{{\mathrm{NO}}_{\mathrm{3}}}+\left[{\mathrm{NO}}_{\mathrm{3}}\right]}}\cdot {\displaystyle \frac{{\mathit{\kappa }}_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{i}}}{{\mathit{\kappa }}_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{i}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}},\text{(10)}& {\displaystyle }{R}_{{\mathrm{SO}}_{\mathrm{4}}}={\displaystyle \frac{\left[{\mathrm{SO}}_{\mathrm{4}}\right]}{{\mathit{\kappa }}_{{\mathrm{SO}}_{\mathrm{4}}}+\left[{\mathrm{SO}}_{\mathrm{4}}\right]}}\cdot {\displaystyle \frac{{\mathit{\kappa }}_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{i}}}{{\mathit{\kappa }}_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{i}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\cdot {\displaystyle \frac{{\mathit{\kappa }}_{{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{i}}}{{\mathit{\kappa }}_{{\mathrm{NO}}_{\mathrm{3}}}+\left[{\mathrm{NO}}_{\mathrm{3}}\right]}},\end{array}$$

with the exception that in the biogeochemical configuration used here we do not consider nitrate (NO${}_{\mathrm{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 ${\mathit{\kappa }}_i^{\mathrm{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.

3 Oceanic methane cycling

3.1 Microbial methanogenesis

Methanogenesis represents the terminal step in our remineralization scheme and follows the overall stoichiometry:

$${\mathrm{2}\mathrm{CH}}_{\mathrm{2}}\mathrm{O}\to {\mathrm{CH}}_{\mathrm{4}}{+\mathrm{CO}}_{\mathrm{2}}.$$

This can be taken to implicitly include fermentation of organic matter to acetate followed by acetoclastic methanogenesis:

$$\begin{array}{c}{\displaystyle }{\mathrm{2}\mathrm{CH}}_{\mathrm{2}}\mathrm{O}\to {\mathrm{CH}}_{\mathrm{3}}\mathrm{COOH},\\ {\displaystyle }{\mathrm{CH}}_{\mathrm{3}}\mathrm{COOH}\to {\mathrm{CH}}_{\mathrm{4}}{+\mathrm{CO}}_{\mathrm{2}},\end{array}$$

or the fermentation of organic matter to acetate followed by anaerobic acetate oxidation and hydrogenotrophic methanogenesis:

$$\begin{array}{c}{\displaystyle }{\mathrm{2}\mathrm{CH}}_{\mathrm{2}}\mathrm{O}\to {\mathrm{CH}}_{\mathrm{3}}\mathrm{COOH},\\ {\displaystyle }{\mathrm{CH}}_{\mathrm{3}}{\mathrm{COOH}+\mathrm{2}\mathrm{H}}_{\mathrm{2}}\mathrm{O}\to {\mathrm{4}\mathrm{H}}_{\mathrm{2}}+{\mathrm{2}\mathrm{CO}}_{\mathrm{2}},\\ {\displaystyle }{\mathrm{4}\mathrm{H}}_{\mathrm{2}}{+\mathrm{CO}}_{\mathrm{2}}\to {\mathrm{CH}}_{\mathrm{4}}{+\mathrm{2}\mathrm{H}}_{\mathrm{2}}\mathrm{O},\end{array}$$

both of which have the same overall net stoichiometry provided that H₂ is assumed to be quantitatively converted to CH₄ by hydrogenotrophic methanogens. We thus ignore the scenario in which some fraction of H₂ 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₄ (e.g., $R_{{\mathrm{CH}}_{\mathrm{4}}}=\mathrm{1}-R_{{\mathrm{O}}_{\mathrm{2}}}-R_{{\mathrm{NO}}_{\mathrm{3}}}-R_{{\mathrm{SO}}_{\mathrm{4}}}$):

$$\begin{array}{}\text{(11)}& R_{{\mathrm{CH}}_{\mathrm{4}}}={\displaystyle \frac{{\mathit{\kappa }}_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{i}}}{{\mathit{\kappa }}_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{i}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\cdot {\displaystyle \frac{{\mathit{\kappa }}_{{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{i}}}{{\mathit{\kappa }}_{{\mathrm{NO}}_{\mathrm{3}}}+\left[{\mathrm{NO}}_{\mathrm{3}}\right]}}\cdot {\displaystyle \frac{{\mathit{\kappa }}_{{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{i}}}{{\mathit{\kappa }}_{{\mathrm{SO}}_{\mathrm{4}}}+\left[{\mathrm{SO}}_{\mathrm{4}}\right]}},\end{array}$$

where the κ_i and ${\mathit{\kappa }}_i^{\mathrm{i}}$ terms are as described above (Table 1). We disable nitrate (NO₃) 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₂] values significantly above 1 µmol kg⁻¹ (Fig. 1a), while anaerobic remineralization is dominated by methanogenesis at [SO${}_{\mathrm{4}}^{\mathrm{2}-}$] values significantly below 1 mmol kg⁻¹ (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₂] in the ocean water column (Fig. 1c).

Figure 1 Fractional organic carbon remineralization by aerobic respiration, sulfate reduction, and methanogenesis in our modified organic matter remineralization scheme. In (a), relative rates of aerobic (O₂) and anaerobic (SO${}_{\mathrm{4}}^{\mathrm{2}-}+$ CH₄) remineralization are plotted as a function of dissolved [O₂]. In (b), relative anaerobic remineralization rates are partitioned between sulfate reduction and methanogenesis as a function of dissolved [SO${}_{\mathrm{4}}^{\mathrm{2}-}$] (dissolved [O₂] is held constant at 10⁻¹⁰ mol kg⁻¹). 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).

While the model tracks the carbon isotope composition of oceanic and atmospheric CH₄ (δ¹³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₄ during methanogenesis of −35 ‰ by default (Table 2), which will tend to produce microbial CH₄ with a δ¹³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₄, or photochemical breakdown of CH₄ in the atmosphere. It does, however, include a comprehensive ¹³C scheme associated with ocean–atmosphere cycling of CO₂ (Kirtland Turner and Ridgwell, 2016; Ridgwell, 2001).

Table 2 Default kinetic and thermodynamic parameters for oceanic methane cycling. Activity coefficients are estimated for T=25 ^∘C and S=35 ‰.

¹ Olson et al. (2016); ² Regnier et al. (2011); ³ Schink (1997); ⁴ Hoehler et al. (2001); ⁵ Hoehler (2004); ⁶ Stoessell and Byrne (1982); ⁷ Cramer (1984); ⁸ Duan et al. (1992); ⁹ Johnson (1982); ¹⁰ Clegg and Brimblecombe (1990); ¹¹ Ulfsbo et al. (2015); ¹² Berner (1965); ¹³ Helz et al. (2011); ¹⁴ Dale et al. (2008).

3.2 Aerobic methanotrophy

Microbial aerobic methanotrophy proceeds according to

$${\mathrm{CH}}_{\mathrm{4}}+{\mathrm{2}\mathrm{O}}_{\mathrm{2}}\to {\mathrm{CO}}_{\mathrm{2}}+{\mathrm{2}\mathrm{H}}_{\mathrm{2}}\mathrm{O}.$$

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 (Jin and Bethke, 2005, 2007; Regnier et al., 2011), in which CH₄ oxidation kinetics are controlled by substrate availability, thermodynamic energy yield, and temperature:

$$\begin{array}{}\text{(12)}& R_{\mathrm{AER}}=k_{\mathrm{AER}}\cdot F_k^{\mathrm{AER}}\cdot F_t^{\mathrm{AER}}\cdot F_{\mathrm{T}}.\end{array}$$

A rate constant for aerobic methanotrophy (yr⁻¹) 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.

The kinetic factor (F_k) for aerobic methanotrophy is controlled by substrate availability according to

$$\begin{array}{}\text{(13)}& F_k^{\mathrm{AER}}=\left[{\mathrm{CH}}_{\mathrm{4}}\right]\cdot {\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{{\mathit{\kappa }}_{\mathrm{O}}^{\mathrm{AER}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}},\end{array}$$

where brackets denote concentration and the κ term denotes a half-saturation constant with respect to O₂. We employ a hybrid parameterization in which kinetics are first-order with respect to CH₄ but also scaled by a Michaelis–Menten-type term for O₂. This formulation is based on the rationale that half-saturation constants for CH₄ are typically similar to (or greater than) the dissolved CH₄ levels attained in anoxic water column environments (Regnier et al., 2011) but is also meant to allow for rapid CH₄ consumption under “bloom” conditions with an appropriately scaled rate constant (see below).

The effect of thermodynamic energy yield on aerobic methanotrophy is given by

$$\begin{array}{}\text{(14)}& F_t^{\mathrm{AER}}=\mathrm{1}-\exp \left[{\displaystyle \frac{\mathrm{\Delta }{G}_{\mathrm{r},\mathrm{AER}}+\mathrm{\Delta }{G}_{\mathrm{BQ},\mathrm{AER}}}{\mathit{\chi }RT}}\right],\end{array}$$

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 rate-determining 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

$$\begin{array}{}\text{(15)}& \mathrm{\Delta }{G}_{\mathrm{r},\mathrm{AER}}=\mathrm{\Delta }{G}_{\mathrm{r},\mathrm{AER}}^{\mathrm{0}}+RT\cdot \ln {\displaystyle \frac{{\mathit{\gamma }}_{{\mathrm{CO}}_{\mathrm{2}}}\left[{\mathrm{CO}}_{\mathrm{2}}\right]}{{\mathit{\gamma }}_{{\mathrm{O}}_{\mathrm{2}}}\left[{\mathrm{O}}_{\mathrm{2}}\right]\cdot {\mathit{\gamma }}_{{\mathrm{CH}}_{\mathrm{4}}}\left[{\mathrm{CH}}_{\mathrm{4}}\right]}},\end{array}$$

where, in addition to the terms defined above, $\mathrm{\Delta }{G}_{\mathrm{r}}^{\mathrm{0}}$ represents the Gibbs free energy of the reaction under standard conditions, and γ_i values represent activity coefficients. Note that we assume an H₂O activity of unity.

3.3 Anaerobic oxidation of methane (AOM)

The oxidation of methane can also be coupled to electron acceptors other than O₂, including nitrate (NO${}_{\mathrm{3}}^{-}$), sulfate (SO${}_{\mathrm{4}}^{\mathrm{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 SO${}_{\mathrm{4}}^{\mathrm{2}-}$:

$$\begin{array}{}\text{(16)}& {\mathrm{CH}}_{\mathrm{4}}+{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\to {\mathrm{HCO}}_{\mathrm{3}}^{-}+{\mathrm{HS}}^{-}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}.\end{array}$$

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₄ at the sulfate–methane transition zone (SMTZ) represents an extremely large sink flux of CH₄ 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₄ (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

$$\begin{array}{}\text{(17)}& R_{\mathrm{AOM}}=k_{\mathrm{AOM}}\cdot F_k^{\mathrm{AOM}}\cdot F_t^{\mathrm{AOM}}\cdot F_{\mathrm{T}},\end{array}$$

where k_AOM is a rate constant for anaerobic methane oxidation (yr⁻¹), 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

$$\begin{array}{}\text{(18)}& F_k^{\mathrm{AOM}}=\left[{\mathrm{CH}}_{\mathrm{4}}\right]\cdot {\displaystyle \frac{\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]}{{\mathit{\kappa }}_{\mathrm{S}}^{\mathrm{AOM}}+\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]}},\end{array}$$

where brackets denote concentration and the κ term denotes a half-saturation constant with respect to SO${}_{\mathrm{4}}^{\mathrm{2}-}$. We employ a hybrid parameterization in which kinetics are first-order with respect to CH₄ but are also scaled by a Michaelis–Menten-type term for SO${}_{\mathrm{4}}^{\mathrm{2}-}$ for reasons discussed above.

The effect of thermodynamic energy yield on anaerobic methane oxidation is specified as follows:

$$\begin{array}{}\text{(19)}& F_t^{\mathrm{AOM}}=\mathrm{1}-\exp \left[{\displaystyle \frac{\mathrm{\Delta }{G}_{\mathrm{r},\mathrm{AOM}}+\mathrm{\Delta }{G}_{\mathrm{BQ},\mathrm{AOM}}}{\mathit{\chi }RT}}\right].\end{array}$$

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

$$\begin{array}{}\text{(20)}& \begin{array}{rl}\mathrm{\Delta }{G}_{\mathrm{r},\mathrm{AOM}}=& \mathrm{\Delta }{G}_{\mathrm{r},\mathrm{AOM}}^{\mathrm{0}}+RT\\ & \cdot \ln {\displaystyle \frac{{\mathit{\gamma }}_{{\mathrm{HCO}}_{\mathrm{3}}^{-}}\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]\cdot {\mathit{\gamma }}_{{\mathrm{HS}}^{-}}\left[{\mathrm{HS}}^{-}\right]}{{\mathit{\gamma }}_{{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}}\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]\cdot {\mathit{\gamma }}_{{\mathrm{CH}}_{\mathrm{4}}}\left[{\mathrm{CH}}_{\mathrm{4}}\right]}},\end{array}\end{array}$$

where $\mathrm{\Delta }{G}_{\mathrm{r}}^{\mathrm{0}}$ 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₂O activity of unity.

3.4 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₄ 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 CH₄ release experiments could employ rate constants that are scaled upward to reflect transient (bloom) elevations in microbial community CH₄ 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 or mixed culture experiments with aerobic methanotrophs (Bender and 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.

Figure 2 Compilation of rate constants for aerobic (AER; a) and anaerobic (AOM; b) methane oxidation. Rate constants are corrected for in situ temperature using a Q₁₀ of 2 (see the Supplement). Vertical red lines show our default values as reported in Table 2. Shown in (c) are globally integrated diffusive fluxes of CH₄ 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).

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⁻¹ 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₄ (J_gas) are governed by temperature- and salinity-dependent solubility and surface wind speed above a given grid cell:

$$\begin{array}{}\text{(21)}& J_{\mathrm{gas}}=A\cdot k_{\mathrm{gas}}\cdot \left([{\mathrm{CH}}_{\mathrm{4}}{]}_{\mathrm{sat}}-[{\mathrm{CH}}_{\mathrm{4}}{]}_{\mathrm{cell}}\right),\end{array}$$

where A denotes the area available for gas exchange (e.g., the area of ice-free surface ocean), [CH₄]_cell denotes the ambient dissolved CH₄ concentration in a given surface ocean grid cell, [CH₄]_sat represents the dissolved CH₄ concentration at saturation with a given atmospheric pCH₄, 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

$$\begin{array}{}\text{(22)}& \begin{array}{rl}\ln \mathit{\beta }=& a_{\mathrm{1}}+a_{\mathrm{2}}(\mathrm{100}/T)+a_{\mathrm{3}}\ln (T/\mathrm{300})\\ & +S\left[b_{\mathrm{1}}+b_{\mathrm{2}}(T/\mathrm{100})+b_{\mathrm{3}}(T/\mathrm{100}{)}^{\mathrm{2}}\right].\end{array}\end{array}$$

Note that the Henry's law constant K₀ is related to the Bunsen solubility coefficient by $K_{\mathrm{0}}=\mathit{\beta }/\mathit{\rho }{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 ‰:

$$\begin{array}{}\text{(23)}& k_{\mathrm{gas}}=k\cdot u^{\mathrm{2}}\cdot {\left[Sc/\mathrm{660}\right]}^{-\mathrm{0.5}},\end{array}$$

where k is a dimensionless gas transfer coefficient, u is surface wind speed, and Sc is the temperature-corrected Schmidt number according to

$$\begin{array}{}\text{(24)}& Sc=c_{\mathrm{1}}-c_{\mathrm{2}}T+c_{\mathrm{3}}{T}^{\mathrm{2}}-c_{\mathrm{4}}{T}^{\mathrm{3}}.\end{array}$$

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₄ follows by default that for other gases accounted for in BIOGEM, such as O₂ and CO₂, as described in Ridgwell et al. (2007)

Table 3 Default constants and coefficients for CH₄ gas exchange. All default parameter values are derived from Wanninkhof (1992). Schmidt number coefficients are for S=35 ‰.

4.2 Parameterized O₂–O₃–CH₄ photochemistry

Once degassed to the atmosphere, CH₄ becomes involved in a complex series of photochemical reactions initiated by hydroxyl radical (OH) attack on CH₄ (Kasting et al., 1983; Prather, 1996; Pavlov et al., 2000; Schmidt and Shindell, 2003). Following Claire et al. (2006) and Goldblatt et al. (2006), we parameterize O₂–O₃–CH₄ photochemistry according to a bimolecular “rate law”:

$$\begin{array}{}\text{(25)}& J_{{\mathrm{CH}}_{\mathrm{4}}}=k_{\mathrm{eff}}\cdot M_{{\mathrm{O}}_{\mathrm{2}}}\cdot M_{{\mathrm{CH}}_{\mathrm{4}}},\end{array}$$

where M_i terms represent the atmospheric inventories of O₂ and CH₄, respectively, and k_eff denotes an effective rate constant (Tmol⁻¹ yr⁻¹) that is itself a complicated function of atmospheric O₂, CH₄, and CO₂ (Claire et al., 2006). We note that in this parameterization, O₃ abundance is not calculated explicitly, but rather the photochemical destruction rate of CH₄ 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₂ and CH₄ 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₂ and CH₄ and a constant atmospheric CO₂ of 10⁻² bar. We then fit a fifth-order polynomial surface to these k_eff values as a function of atmospheric pO₂ and pCH₄ (Fig. 3).

Figure 3 Shown in (a) is the bivariate fit to a suite of 1-D atmospheric photochemical runs for the effective rate constant (k_eff) parameterizing O₂–O₃–CH₄ photochemistry in ATCHEM. Shown in (b) is a frequency distribution of the residuals on k_eff from the underlying photochemical model output.

Our default parameterization of O₂–O₃–CH₄ chemistry (C06) is fit over a pO₂ range of 10⁻¹⁴ to 10⁻¹ bar, a pCH₄ range of 10⁻⁶ to $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$ bar, and a constant high background pCO₂ of 10⁻² bar (Claire et al., 2006). We thus truncate the atmospheric lifetime of CH₄ at a lower bound of 7.6 years in our default parametrization and provide an alternative parameterization of photochemical CH₄ destruction at roughly modern pO₂ and pCO₂ (SS03) derived from the results of Schmidt and Shindell (2003) for use in more geologically recent, high-O₂ 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₄ in the atmosphere, including atmospheric pCO₂, the atmospheric profile of H₂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.

Figure 4 Comparison of steady-state atmospheric pCH₄ as a function of terrestrial CH₄ 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₂=0.1 atm (black crosses; see text). Shown in (c) are steady-state atmospheric CH₄ values as a function of imposed terrestrial CH₄ flux in our “modern” configuration (circles) compared to estimates for the glacial, preindustrial, and modern CH₄ cycles (Kirschke et al., 2013; Bock et al., 2017; Paudel et al., 2016).

As a basic test of our photochemical parameterization, we impose a terrestrial (wetland) flux of CH₄ to the atmosphere (balanced by stoichiometric consumption of CO₂ and release of O₂) and allow the oceanic and atmospheric CH₄ cycle to spin up for 20 kyr. We then compare steady-state atmospheric pCH₄ as a function of terrestrial CH₄ 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₂ and CH₄ inventories. Nevertheless, both the default scheme and the alternative parameterization for recent geologic history (and analogous planetary environments), with high-pO₂,low-pCO₂ atmospheres, accurately reproduce atmospheric pCH₄ values given estimated glacial, preindustrial, and modern terrestrial CH₄ fluxes (Fig. 4c), and both display the predicted saturation of CH₄ sinks at elevated atmospheric CH₄ observed in more complex photochemical models. We note, however, that the alternative parameterization tends to yield slightly higher atmospheric pCH₄ at surface fluxes greater than ∼ 50 Tmol yr⁻¹ (Fig. 4c). (In the remainder of the paper, we employ the default (C06) parameterization for atmospheric O₂–O₃–CH₄ chemistry and do not discuss the simple high-pO₂−low-pCO₂ alternative further.)

5 Example applications of the new capabilities in the cGENIE model

5.1 High-pO₂ (“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₂, pCO₂, and pCH₄ of [v∕v] 20.95 %, 278 ppm, and 700 ppb, respectively, and with globally uniform oceanic concentrations of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ (28 mmol kg⁻¹) and CH₄ (1 nmol kg⁻¹). We fix globally averaged solar insolation at the modern value (1368 W m⁻²) with seasonally variable forcing as a function of latitude and set radiative forcing for CO₂ and CH₄ 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₂ and pCO₂ (and δ¹³C of atmospheric CO₂) restored to preindustrial values at every time step and with an imposed wetland flux of CH₄ to the atmosphere of 20 Tmol yr⁻¹ that has a δ¹³C value of −60 ‰. Atmospheric pCH₄ and all oceanic tracers are allowed to evolve freely.

Surface, benthic, and ocean interior distributions of dissolved oxygen (O₂), sulfate (SO${}_{\mathrm{4}}^{\mathrm{2}-}$), and methane (CH₄) are shown in Fig. 5 for our roughly modern simulation. Dissolved O₂ ([O₂]) 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₂] 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₂] 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.

Figure 5 Tracer distributions in surface (a–c) and benthic (d–f) grid cells and in the zonally averaged ocean interior (g–i) for O₂ (a, d, g), SO${}_{\mathrm{4}}^{\mathrm{2}-}$ (b, e, h), and CH₄ (c, f, i) in our “modern” configuration. Note the different concentration units for each tracer.

Concentrations of dissolved SO${}_{\mathrm{4}}^{\mathrm{2}-}$ ([SO${}_{\mathrm{4}}^{\mathrm{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 [SO${}_{\mathrm{4}}^{\mathrm{2}-}$] distributions show some similarity to those of [O₂] (Fig. 5e), though again the differences are very small relative to the overall prescribed initial tracer inventory of 28 mmol kg⁻¹ and disappear almost entirely when salinity-normalized (not shown). In the ocean interior, [SO${}_{\mathrm{4}}^{\mathrm{2}-}$] is largely spatially invariant with a value of approximately 28 mmol kg⁻¹ (Fig. 5h).

Dissolved CH₄ concentrations ([CH₄]) in the surface and shallow subsurface ocean are much more variable but typically on the order of ∼ 1–2 nmol kg⁻¹, 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₄] distribution shows locally elevated values up to ∼ 300–400 nmol kg⁻¹ 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₄ cycling (Jayakumar et al., 2001). Within the ocean interior, dissolved CH₄ can accumulate in the water column in excess of ∼ 100 nmol kg⁻¹ in association with relatively low-O₂ conditions at intermediate depths, with zonally averaged values as high as ∼ 70 nmol kg⁻¹ but more typically in the range of ∼ 20–40 nmol kg⁻¹ (Fig. 5i). These concentrations are comparable to those observed locally in low-O₂ 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₄ cycle for our “modern” configuration are shown in Fig. 6. Methanogenesis is focused in regions characterized by relatively low [O₂] 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₄ 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⁻¹ d⁻¹ compare well with field measurements of AOM within oceanic OMZs (Thamdrup et al., 2019). In general, the bulk of CH₄ produced via microbial methanogenesis is consumed via AOM, either near the seafloor or within the ocean interior.

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.

5.2 Low-pO₂ (“ancient”) steady state

Next, we explore a low-pO₂ 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₂, pCO₂, and pCH₄ of [v∕v] $\mathrm{2.1}\times {\mathrm{10}}^{-\mathrm{4}}$ atm (equivalent to a value 10⁻³ times the present atmospheric level, PAL), 278 ppm, and 500 ppm, respectively, as well as globally uniform oceanic concentrations of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ (280 µmol kg⁻¹) and CH₄ (50 µmol kg⁻¹). We again fix globally averaged solar insolation at the modern value (1368 W m⁻²), with seasonally variable forcing as a function of latitude, and set radiative forcing for CO₂ and CH₄ 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₂ and pCO₂ (and δ¹³C of atmospheric CO₂) restored to the initial values specified above at every time step, with an imposed “geologic” flux of CH₄ to the atmosphere of 3 Tmol yr⁻¹ at a δ¹³C value of −60 ‰. Atmospheric pCH₄ and all oceanic tracers are allowed to evolve freely.

Surface, benthic, and ocean interior distributions of [O₂], [SO${}_{\mathrm{4}}^{\mathrm{2}-}$], and [CH₄] are shown in Fig. 7 for our low-pO₂ simulation. Dissolved O₂ concentrations are now extremely heterogeneous throughout the surface ocean, ranging over an order of magnitude from less than 1 to over 10 µmol kg⁻¹, with concentrations that are regionally well in excess of air saturation at the prescribed pO₂ of $\mathrm{2.1}\times {\mathrm{10}}^{-\mathrm{4}}$ atm (Fig. 7a). Previous studies have shown that these features are not unexpected at very low atmospheric pO₂ (Olson et al., 2013; Reinhard et al., 2016). We note, however, that the distribution and maximum [O₂] in our low-pO₂ 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 PO${}_{\mathrm{4}}^{\mathrm{3}-}$ 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₂] dynamics at low atmospheric pO₂ (Olson et al., 2013; Reinhard et al., 2016), our regional [O₂] patterns still generally track the localized balance between photosynthetic O₂ release and consumption through respiration and reaction with inorganic reductants, rather than temperature-dependent solubility patterns (Fig. 5a). Within the ocean interior, O₂ is consumed within the upper few hundred meters and is completely absent in benthic settings (Fig. 7d, g).

Figure 7 Tracer distributions in surface (a–c) and benthic (d–f) grid cells and in the zonally averaged ocean interior (g–i) for O₂ (a, d, g), SO${}_{\mathrm{4}}^{\mathrm{2}-}$ (b, e, h), and CH₄ (c, f, i) in our “ancient” configuration (see text). Note the different concentration units for each tracer and the differing scales relative to Fig. 5.

In our low-pO₂ simulations we initialize the ocean with a globally uniform [SO${}_{\mathrm{4}}^{\mathrm{2}-}$] of 280 µmol kg⁻¹, under the premise that the marine SO${}_{\mathrm{4}}^{\mathrm{2}-}$ inventory should scale positively with atmospheric pO₂. With this much lower initial SO${}_{\mathrm{4}}^{\mathrm{2}-}$ inventory (i.e., 10² times less than the modern ocean), steady-state [SO${}_{\mathrm{4}}^{\mathrm{2}-}$] distributions are significantly more heterogeneous than in the modern, high-pO₂ case (Fig. 7). Ocean [SO${}_{\mathrm{4}}^{\mathrm{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-pO₂ simulations SO${}_{\mathrm{4}}^{\mathrm{2}-}$ serves as the principal oxidant for organic matter remineralization in the ocean interior, with the result that its distribution effectively mirrors that of [O₂] in the modern case in both spatial texture and overall magnitude (compare Fig. 7e, h with Fig. 5d, g). Dissolved SO${}_{\mathrm{4}}^{\mathrm{2}-}$ in this simulation never drops to zero, a consequence of our initial 280 µmol kg⁻¹ concentration of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ representing the oxidative potential of 560 µmol kg⁻¹ of O₂, some 3 times higher than the mean [O₂] value in the modern ocean interior (∼ 170 µmol kg⁻¹).

Dissolved CH₄ 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⁻¹ (Fig. 7c). The benthic [CH₄] distribution shows concentrations up to ∼ 8 µmol kg⁻¹, with concentrations in excess of 1 µmol kg⁻¹ pervasively distributed across the seafloor. In general, the benthic [CH₄] distribution inversely mirrors that of [SO${}_{\mathrm{4}}^{\mathrm{2}-}$] (Fig. 7f), which results from the fact that in the low-pO₂ case SO${}_{\mathrm{4}}^{\mathrm{2}-}$ again serves as the principal oxidant of methane. Concentrations of CH₄ in the ocean interior can approach ∼ 10 µmol kg⁻¹ but in the zonal average are typically less than 5 µmol kg⁻¹ (Fig. 7i). Overall, the oceanic CH₄ inventory increases dramatically in the low-pO₂ case relative to the modern simulation, from ∼ 4.5 Tmol CH₄ to ∼ 1900 Tmol CH₄.

The major metabolic fluxes within the ocean's microbial CH₄ cycle for our “ancient” configuration are shown in Fig. 8. Column-integrated rates of microbial methanogenesis are greater than in the high-pO₂ case by up to a factor of ∼ 10² (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₂ 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₂ 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² times higher than those observed in the high-pO₂ case (Fig. 8c, f). Once again, AOM dominates the consumption of CH₄ produced in the ocean interior and is extremely effective at reducing CH₄ fluxes to the atmosphere. Despite a significant increase in overall oceanic CH₄ burden relative to our high-pO₂ simulation (see above and Fig. 7i), atmospheric pCH₄ increases only modestly from ∼ 0.8 ppm to 6 ppm [v∕v], equivalent to an additional radiative forcing of only ∼ 2 W m⁻², due to efficient microbial consumption in the upper ocean.

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.

5.3 Atmospheric carbon injection

To illustrate the capabilities of the model in exploring the time-dependent (perturbation) behavior of the CH₄ cycle, we perform a simple carbon injection experiment in which 3000 PgC is injected directly into the atmosphere either as CH₄ or as CO₂, 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⁻¹, 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₂ and/or CH₄ (Kirtland Turner, 2018). This flux is much lower than the current anthropogenic carbon input of ∼ 10 PgC yr⁻¹ (Ciais et al., 2013). 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₄, there is an immediate and significant increase in atmospheric pCH₄ to values greater than 10 ppmv, followed by a gradual increase to a maximum of ∼ 12 ppmv throughout the duration of the CH₄ 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₂ (Fig. 9b). This increase in atmospheric pCH₄ and pCO₂ 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).

Figure 9 Response to a 3000 PgC release directly to the atmosphere spread over 1000 years, assuming carbon is injected as either CH₄ or CO₂. Atmospheric pCH₄ (a), pCO₂ (b), mean surface ocean pH (c), mean surface air temperature (SAT; d), and mean ocean temperature (MOT; e) are shown for a CH₄ injection (grey) and a CO₂ injection (black). Panel (f) shows the difference in SAT and MOT between the CH₄ and CO₂ injection scenarios ($\mathrm{\Delta }{T}_{\mathrm{i}}=T_{{\mathrm{CH}}_{\mathrm{4}},\mathrm{i}}$ – $T_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{i}}$) through time.

The increase in atmospheric pCO₂ and drop in ocean pH are nearly identical if we instead inject the carbon as CO₂ rather than CH₄ (Fig. 9b, c). However, when carbon is injected as CH₄, 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, CH₄ is a much more powerful greenhouse gas than is CO₂, and oxidation of CH₄ to CO₂ 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₄ rather than CO₂, providing an enhancement of warming, especially during the duration of carbon input. This warming enhancement should be considered in past events during which CH₄ release is suspected as a key driver of warming. For instance, additional warming due to CH₄ 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).

5.4 Atmospheric pCH₄ on the early Earth

Using our low-pO₂ steady state as a benchmark case (Sect. 5.2), we briefly explore the sensitivity of atmospheric pCH₄ to a subset of model variables. All model ensembles are initially configured with globally homogeneous marine SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and CH₄ inventories, as well as a background geologic CH₄ flux of 3 Tmol yr⁻¹, and are spun up for 20 kyr with a fixed pO₂ and pCO₂. We report atmospheric pCH₄ 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₄ on a low-oxygen Earth-like planet. We present results from individual sensitivity ensembles from our benchmark low-pO₂ case over the following parameter ranges: (1) atmospheric pO₂ between 10⁻⁴ and 10⁻¹ times the present atmospheric level (PAL), equivalent to roughly $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{2}}$ atm, respectively; (2) initial marine SO${}_{\mathrm{4}}^{\mathrm{2}-}$ inventories corresponding to globally uniform seawater concentrations between 0 and 1000 µmol kg⁻¹; and (3) biological energy quanta (BEQ) for anaerobic methane oxidation between 5 and 30 kJ mol⁻¹.

Results for our low-pO₂ sensitivity ensembles are shown in Fig. 10. We find a similar sensitivity of atmospheric pCH₄ to atmospheric pO₂ to that observed by Olson et al. (2016). In particular, atmospheric CH₄ abundance initially increases as atmospheric pO₂ drops below modern values to roughly 2 %–3 % PAL, after which decreasing pO₂ causes pCH₄ 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₂O by ozone (O₃) decreases at low atmospheric pO₂ (Pavlov et al., 2003; Claire et al., 2006; Goldblatt et al., 2006). However, peak atmospheric pCH₄ is significantly reduced in our models relative to those of Olson et al. (2016). For example, at an “optimal” atmospheric pO₂ of ∼ 2.5 % PAL, Olson et al. (2016) predict a steady-state atmospheric pCH₄ of ∼ 35 ppmv, while we predict a value of ∼ 10 ppmv (Fig. 10a). This difference can be attributed to our updated O₂–O₃–CH₄ 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₄ as a significant climate regulator at steady state during most of the Proterozoic eon (between ∼ 2.5 and 0.5 Ga).

Figure 10 Sensitivity ensembles of our “ancient” configuration compared to the results of Olson et al. (2016). Steady-state atmospheric pCH₄ values as a function of assumed atmospheric pO₂ (a) and initial marine SO${}_{\mathrm{4}}^{\mathrm{2}-}$ 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₄. All simulations were spun up from cold for 20 kyr, with the results shown from the last model year.

Atmospheric CH₄ abundance is also strongly sensitive to the marine SO${}_{\mathrm{4}}^{\mathrm{2}-}$ inventory (Fig. 10b). The scaling we observe between the initial SO${}_{\mathrm{4}}^{\mathrm{2}-}$ inventory and steady-state atmospheric pCH₄ is very similar to that reported by Olson et al. (2016), with a sharp drop in the marine CH₄ inventory and atmospheric CH₄ abundance as marine SO${}_{\mathrm{4}}^{\mathrm{2}-}$ drops below ∼ 100 µmol kg⁻¹ (Fig. 10b). The implication is that for most of Earth's history anaerobic oxidation of CH₄ in the ocean interior has served as an important inhibitor of CH₄ 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${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations may instead have been on the order of ∼ 1–10 µmol kg⁻¹ (Crowe et al., 2014), while atmospheric pO₂ 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₄ and Earth's climate system may thus have been much more important prior to ∼ 2.5 billion years ago.

Atmospheric CH₄ 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₄ from ∼ 7 to ∼ 25 ppmv when increasing the BEQ value from 20 to 30 kJ mol⁻¹ (Fig. 10c). This effect is mediated primarily by the importance of anaerobic methanotrophy when atmospheric pO₂ 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₄-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₄ cycle and seawater [SO${}_{\mathrm{4}}^{\mathrm{2}-}$] is relatively low.

6 Discussion and conclusions

The global biogeochemical cycling of CH₄ 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₄ 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₄ 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₄ cycle in Earth's past and future and ultimately to constrain CH₄ 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 CH₄ production. We have also incorporated revised schemes for microbial CH₄ consumption that include both kinetic and thermodynamic constraints, and we have updated the parameterized atmospheric O₂–O₃–CH₄ photochemistry to improve accuracy and for use across a wider range of atmospheric pO₂ values than explored in previous work. Simulations of roughly modern (high-O₂) and Proterozoic (low-O₂) Earth system states demonstrate that the model effectively reproduces the first-order features of the modern ocean–atmosphere CH₄ cycle and can be effectively implemented across a wide range of atmospheric O₂ partial pressures and marine SO${}_{\mathrm{4}}^{\mathrm{2}-}$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₄ consumption may have been important in regulating atmospheric CH₄ abundance during the Archean. Finally, our simulation of PETM-like carbon injection demonstrates the importance of explicitly considering CH₄ 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₄ production and 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 pCO₂ to be explored; (3) coupling of deep ocean chemistry with a description of marine methane hydrates and associated sedimentary CH₄ cycling; and (4) developing a representation of the production and consumption of CH₄ by terrestrial ecosystems.