Articles | Volume 15, issue 20
Model description paper
20 Oct 2022
Model description paper |  | 20 Oct 2022

CANOPS-GRB v1.0: a new Earth system model for simulating the evolution of ocean–atmosphere chemistry over geologic timescales

Kazumi Ozaki, Devon B. Cole, Christopher T. Reinhard, and Eiichi Tajika

A new Earth system model of intermediate complexity – CANOPS-GRB v1.0 – is presented for use in quantitatively assessing the dynamics and stability of atmospheric and oceanic chemistry on Earth and Earth-like planets over geologic timescales. The new release is designed to represent the coupled major element cycles of C, N, P, O, and S, as well as the global redox budget (GRB) in Earth's exogenic (ocean–atmosphere–crust) system, using a process-based approach. This framework provides a mechanistic model of the evolution of atmospheric and oceanic O2 levels on geologic timescales and enables comparison with a wide variety of geological records to further constrain the processes driving Earth's oxygenation. A complete detailed description of the resulting Earth system model and its new features are provided. The performance of CANOPS-GRB is then evaluated by comparing a steady-state simulation under present-day conditions with a comprehensive set of oceanic data and existing global estimates of bio-element cycling. The dynamic response of the model is also examined by varying phosphorus availability in the exogenic system. CANOPS-GRB reliably simulates the short- and long-term evolution of the coupled C–N–P–O2–S biogeochemical cycles and is generally applicable across most period of Earth's history given suitable modifications to boundary conditions and forcing regime. The simple and adaptable design of the model also makes it useful to interrogate a wide range of problems related to Earth's oxygenation history and Earth-like exoplanets more broadly. The model source code is available on GitHub and represents a unique community tool for investigating the dynamics and stability of atmospheric and oceanic chemistry on long timescales.

1 Introduction

A quarter century has passed since the first discovery of exoplanets (Mayor and Queloz, 1995). In the next quarter century, a full-scale search for signs of life – biosignatures– on Earth-like exoplanets is one of the primary objectives of the next generation of exoplanetary observational surveys (National Academies of Sciences and Medicine, 2019; The LUVOIR Team, 2019). The definition of biosignatures includes a variety of signatures that require biological activity for their origin (Des Marais et al., 2002; Lovelock, 1965; Sagan et al., 1993; Schwieterman et al., 2018; National Academies of Sciences and Medicine, 2019), but atmospheric composition has received the most interdisciplinary attention since the dawn of the search for life beyond our own planet (Hitchcock and Lovelock, 1967; Lovelock, 1965, 1972, 1975; Sagan et al., 1993) because of its potential for remote detectability. Indeed, it is likely that deciphering of exoplanetary atmospheric composition based on spectroscopic information will, at least for the foreseeable future, be our only promising means for life detection beyond our solar system. However, the detection of atmospheric composition cannot immediately answer the question of the presence or absence of a surface biosphere because significant gaps remain in our understanding of the relationships between atmospheric composition and biological activity occurring at the surface on life-bearing exoplanets. Many of these gaps arise from a lack of robust theoretical and quantitative frameworks for the emergence and maintenance of remotely detectable atmospheric biosignatures in the context of planetary biogeochemistry.

It is also important to emphasize that the abundance of atmospheric biosignature gases of living planets will evolve via an intimate interaction between life and global biogeochemical cycles of bio-essential elements across a range of timescales. Indeed, the abundances of biosignature gases such as molecular oxygen (O2) and methane (CH4) in Earth's atmosphere have evolved dramatically through coevolutionary interaction with Earth's biosphere for nearly 4 billion years – through remarkable fluctuations in atmospheric chemistry and climate (Catling and Kasting, 2017; Lyons et al., 2014; Catling and Zahnle, 2020). To the extent that the coupled evolution of life and the atmosphere is a universal property of life-bearing planets that maintain robust atmospheric biosignatures, the construction of a biogeochemical framework for diagnosing atmospheric biosignatures should be a subject of urgent interdisciplinary interest.

Establishing a mechanistic understanding of our own planet's evolutionary history is also an important milestone for the construction of a search strategy for life beyond our solar system, as it provides the first step towards understanding how remotely detectable biosignatures emerge and are maintained on a planetary scale. While numerous atmospheric biosignature gases have been proposed, the most promising candidates have been “redox-based” species, such as O2, ozone (O3), and CH4 (Meadows, 2017; Meadows et al., 2018; Reinhard et al., 2017a; Krissansen-Totton et al., 2018). In particular, O2 is of great interest to astrobiologists because of its crucial role in metabolism on Earth. Thus, considerable effort has been devoted over recent decades toward quantitatively and mechanistically understanding Earth's oxygenation history. In particular, a recent surge in the generation of empirical records for Earth's redox evolution has yielded substantial progress in our “broad stroke” understanding of Earth's oxygenation history and has shaped our view of biological evolution (Kump, 2008; Lyons et al., 2014). One of the intriguing insights obtained from the accumulated geochemical records is that atmospheric O2 levels might have evolved more dynamically than previously thought – our current paradigm of Earth's oxygenation history suggests that atmospheric O2 levels may have risen and then plummeted during the early Proterozoic and then remained low (probably <10 % of the present atmospheric level; PAL) for much of the ∼1 billion years leading up to the catastrophic climate system perturbations and the initial diversification of complex life during the late Proterozoic.

The possibility of low but “post-biotic” atmospheric O2 levels during the mid-Proterozoic has important ramifications not only for our basic theoretical understanding of long-term O2 cycle stability on a planet with biological O2 production but also for biosignature detectability (Reinhard et al., 2017a). However, our quantitative and mechanistic understanding of the Earth's O2 cycle in deep time is still rudimentary at present. For example, one possible explanation for low atmospheric O2 levels during the mid-Proterozoic is simply a less active or smaller biosphere (Crockford et al., 2018; Derry, 2015; Laakso and Schrag, 2014; Ozaki et al., 2019a). However, mechanisms for regulating biotic O2 generation rates and stabilizing atmospheric O2 levels at low levels on billion-year timescales remain obscure. As a result, the level of atmospheric O2 and its stability during the early mid-Proterozoic are the subject of vigorous debate (Bellefroid et al., 2018; Canfield et al., 2018; Cole et al., 2016; Planavsky et al., 2018, 2016; Tang et al., 2016; Zhang et al., 2016). Perhaps even more importantly, a relatively rudimentary quantitative framework for probing the dynamics and stability of the oxygen cycle leads to the imprecision of geochemical reconstructions of ocean–atmosphere O2 levels.

Planetary atmospheric O2 levels are governed by a kinetic balance between sources and sinks. Feedback arises because the response of source–sink fluxes to changes in atmospheric O2 levels is intimately interrelated to each other. Since the biogeochemical cycles of C, N, P, and S exert fundamental control on the redox budget through nonlinear interactions and feedback mechanisms, a mechanistic understanding of these biogeochemical cycles is critical for understanding Earth's O2 cycle. However, the wide range of timescales that characterize C, N, P, O2, and S cycling through the reservoirs of the Earth system makes it difficult to fully resolve the mechanisms governing the dynamics and stability of atmospheric O2 levels from geologic records. From this vantage, developing new quantitative tools that can explore biogeochemical cycles under conditions very different from those of the present Earth is an important pursuit.

This study is motivated by the conviction that an ensemble of “open” Earth system modeling frameworks with explicit and flexible representation of the coupled C–N–P–O2–S biogeochemical cycles will ultimately be required to fully understand the dynamics and stability of Earth's O2 cycle and its controlling factors. In particular, a coherent mechanistic framework for understanding the global redox (O2) budget (GRB) is critical for filling remaining gaps in our understanding of Earth's oxygenation history and the cause-and-effect relationships with an evolving biosphere. Here, we develop a new Earth system model, named CANOPS-GRB, which implements the coupled biogeochemical cycles of C–N–P–O2–S within the Earth's surface system (ocean–atmosphere–crust). The core of this model is an ocean biogeochemical model, CANOPS (Ozaki et al., 2011; Ozaki and Tajika, 2013; Ozaki et al., 2019a). This model has been used to examine conditions for the development of widespread oceanic anoxia and euxinia during the Phanerozoic (Ozaki et al., 2011; Ozaki and Tajika, 2013; Kashiyama et al., 2011) and to quantitatively constrain biogeochemical cycles during the Precambrian (Cole et al., 2022; Ozaki et al., 2019a, b; Reinhard et al., 2017b). In this study, we extend this model to simulate the biogeochemical dynamics of the coupled ocean–atmosphere–crust system. The model design (such as the complexity of the processes and spatial-temporal resolution of the model) is constrained by the requirement of simulation length (>100 million years) and actual model run time. A lack of understanding of biogeochemistry in deep time and availability and quality of geologic records also limit the model structure. With this in mind, we aim for a comprehensive, simple, yet realistic representation of biogeochemical processes in the Earth system, yielding a unique tool for investigating coupled biogeochemical cycles within the Earth system over a wide range of timescales. We have placed particular emphasis on the development of a global redox budget in the ocean–atmosphere–crust system given its importance in the secular evolution of atmospheric O2 levels. CANOPS-GRB is an initial step towards developing the first large-scale biogeochemistry evolution model suited for the wide range of redox conditions, including explicit consideration of the coupled C–N–P–O2–S cycles and the major biogenic gases in planetary atmospheres (O2 and CH4).

2 Model description

Here we present a full description of a new version of the Earth system model CANOPS – CANOPS-GRB v1.0 – which is designed to facilitate simulation for a wide range of biogeochemical conditions so as to permit quantitative examination of evolving ocean–atmosphere chemistry throughout Earth's history. Below we first describe the concept of model design (Sect. 2.1). Next, we describe the overall structure of the model and the basic design of global biogeochemical cycles (Sect. 2.2 and 2.3). That is followed by a detailed description of each sub-model.

2.1 CANOPS-GRB in the hierarchy of biogeochemical models

A full understanding of Earth's evolving O2 cycle requires a quantitative framework that includes mechanistic links between biological metabolism, ocean–atmosphere chemistry, and geologic processes. Such a framework must also represent the feedbacks between ocean–atmosphere redox state and biogeochemical cycles of redox-dependent bio-essential elements. Over recent decades, considerable progress has been made in quantifying the feedbacks between atmospheric O2 levels and the coupled C–N–P–O2–S biogeochemical cycles over geological timescales (Berner, 2004a; Lasaga and Ohmoto, 2002; Betts and Holland, 1991; Holland, 1978; Bolton et al., 2006; Slomp and Van Cappellen, 2007; Van Cappellen and Ingall, 1994; Colman and Holland, 2000; Belcher and McElwain, 2008). Refinements to our understanding of mechanisms regulating Earth's surface redox state have been implemented in low-resolution box models where the ocean–atmosphere system is expressed by a few boxes (Bergman et al., 2004; Laakso and Schrag, 2014; Lenton and Watson, 2000a, b; Van Cappellen and Ingall, 1996; Handoh and Lenton, 2003; Petsch and Berner, 1998; Claire et al., 2006; Goldblatt et al., 2006; Alcott et al., 2019). These models offer insights into basic system behavior and can illuminate the fundamental mechanisms that exert the most leverage on biogeochemical cycles because of their simplicity, transparency, and low computational demands. However, these model architectures also have important quantitative limitations. For example, with low spatial resolution the modeler needs to assume reasonable (but a priori) relationships relating to internal biogeochemical cycles in the system. For instance, because of a lack of high vertical resolution, oceanic box models (Knox and McElroy, 1984; Sarmiento and Toggweiler, 1984; Siegenthaler and Wenk, 1984) usually overestimate the sensitivity of atmospheric CO2 levels to biological activity at high-latitude surface ocean relative to projections by general circulation models (Archer et al., 2000). Oceanic biogeochemical cycles and chemical distributions are also characterized by strong vertical and horizontal heterogeneities, which have the potential to affect the strength of feedback processes (Ozaki et al., 2011). In other words, the low-resolution box modeling approach might overlook the strength and response of the internal feedback loops. Thus, the development of an ocean model with high resolution of ocean interior and reliable representation of water circulation is preferred to investigate the mechanisms controlling atmospheric O2 levels under conditions very different from those of the modern Earth.

In the last decade, comprehensive Earth system models of intermediate complexity (EMICs) have also been developed and extended to include ocean sediments and global C cycling (Ridgwell and Hargreaves, 2007; Lord et al., 2016). Such models can be integrated over tens of thousands of years, allowing experimentation with hypothetical dynamics of global biogeochemical cycles in the geological past (Reinhard et al., 2020; Olson et al., 2016). However, a key weakness of existing EMICs is the need to parameterize (or ignore) boundary (input/output) fluxes – either due to the computational expense of explicitly specifying boundary conditions or due to poorly constrained parameterizations. For example, the oceanic P cycle is usually treated as a closed system, limiting the model's applicability to timescales less than the oceanic P residence time (∼15–20 kyr). Further, boundary conditions such as continental configuration and oceanic bathymetry are variable or poorly constrained in deep time, and the use of highly complex models is difficult because of the computational cost. Finally, exploration of hypotheses concerning the biogeochemical dynamics in deep time often require large model ensembles across broad parameter space given the scope of uncertainty. This makes the computational cost of EMICs intractable at present for many key questions.

The CANOPS-GRB model is designed to capture the major components of Earth system biogeochemistry on timescales longer than ∼103 years but is simple enough to allow for runs on the order of 109 model years. The model structure is also designed so that the model captures the essential biogeochemical processes regulating the global O2 budget, while keeping the calculation cost as moderate as possible. For example, the simple relationships of biogeochemical transport processes at the interface of the Earth system (hydrogen escape to space, early diagenesis in marine sediments, and weathering) are employed based on the systematic application of 1-D models in previous studies (Bolton et al., 2006; Daines et al., 2017; Middelburg et al., 1997; Claire et al., 2006; Wallmann, 2003), providing a powerful, computationally efficient means for exploring the Earth system under a wide range of conditions. The resultant CANOPS-GRB model can be run on a standard personal computer on a single CPU with an efficiency of approximately 6 million model years per CPU hour. In other words, model runs in excess of 109 model years are tractable with modest wall times (approximately 7 d). The model is thus not as efficient as simple box models but is highly efficient relative to EMICs, making sensitivity experiments and exploration of larger parameter space over a billion years feasible, particularly with implementation on a high-performance computing cluster (see Cole et al., 2022). CANOPS-GRB thus occupies a unique position within the hierarchy of global biogeochemical cycle models, rendering it a useful tool for the development of more comprehensive, low- to intermediate-complexity models of Earth system on very long timescales.

2.2 Overall model structure

The overall structure of the model is shown in Fig. 1. The model consists of ocean, atmosphere, and sedimentary reservoirs. The core of the model is an ocean model, comprising a high-resolution 1-D intermediate-complexity box model of the global ocean (Sect. 2.4). The ocean model is coupled to a parameterized marine sediment module (Sect. 2.4.4) and a one-box model of the atmosphere (Sect. 2.6). The atmospheric model includes O2 and CH4 as chemical components, and abundances of these molecules are calculated based on the mass balance between sources and sinks (e.g., biogenic fluxes of O2 and CH4 from the ecosystems and photochemical reactions). The net air–sea gas exchange of chemical species (O2, H2S, NH3 and CH4) is quantified according to the stagnant film model (Liss and Slater, 1974; Kharecha et al., 2005) (Sect. 2.4.5). The ocean and atmosphere models are embedded in a “rock cycle” model that simulates the evolution of sedimentary reservoir sizes on geologic timescales (Sect. 2.5). Three sedimentary reservoirs (organic carbon, ORG; pyrite sulfur, PYR; and gypsum sulfur, GYP) are considered in the CANOPS-GRB model. These reservoirs interact with the ocean–atmosphere system through weathering, outgassing, and burial.

Figure 1CANOPS-GRB model configuration.(a) The schematic of material cycles in the surface (ocean–atmosphere–crust) system. Three sedimentary reservoirs, organic carbon (ORG), pyrite sulfur (PYR), and gypsum sulfur (GYP), are considered. Sedimentary reservoirs interact with the ocean–atmosphere system via weathering, volcanic degassing, and burial. No interaction with the mantle is included, except for the input of reduced gases from the mantle. Total mass of sulfur is conserved in the surface system. (b) Schematic of ocean and atmosphere modules. “L” and “H” denote the low- to mid-latitude mixed surface layer and high-latitude surface layer, respectively. An ocean area of 10 % is assumed for H. River flux for each region is proportional to the areal fraction. Ocean interior is divided into two sectors, high-latitude deep water (HD) and low- to mid-latitude deep water (LD), which are vertically resolved. The area of HD is 25 % of the whole ocean. The deep overturning circulation, V˙, equals the poleward flow in the model surface layer (from L to H). Kvl(z) and Kvh(z) are the vertical eddy diffusion coefficients in the LD and HD regions, respectively. Khor and V˙h are the horizontal diffusion coefficient and polar convection, respectively. The black hatch represents the seafloor topography assumed. The parameters regarding geometry and water transport are tabulated in Table 3.


The ocean model is a vertically resolved transport reaction model of the global ocean, which was originally developed by Ozaki et al. (2011) and Ozaki and Tajika (2013). The model consists of 122 boxes across two regions – a low- to mid-latitude region and a high-latitude region (Fig. 1b). The ocean model describes water transport processes as exchange fluxes between boxes and via eddy diffusion terms. More specifically, ocean circulation is modeled as an advection–diffusion model of the global ocean – a general and robust scheme that is capable of producing well-resolved modern profiles of circulation tracers using realistic parameter values (the physical setup of the model can be found in Sect. 2.4.1 and 2.4.2). The biogeochemical sub-model provides a mechanistic description of the marine biogeochemical cycles of C, P, N, O2, and S (Sect. 2.4.3). This includes explicit representation of a variety of biogeochemical processes such as biological productivity in the sunlit surface oceans; a series of respiration pathways and secondary redox reactions under oxic and anoxic conditions (Sect. 2.4.3); and deposition, decomposition, and burial of biogenic materials in marine sediments (Sect. 2.4.4), allowing a mechanistically based examination of biogeochemical processes. The suite of metabolic reactions included in the model is listed in Table 1. Ocean biogeochemical tracers considered in the CANOPS-GRB model are phosphate (PO43-), nitrate (NO3-), total ammonia (ΣNH3), dissolved oxygen (O2), sulfate (SO42-), total sulfide (ΣH2S), and methane (CH4). Note that biogeochemical cycling of trace metals (e.g., Fe and Mn) is not included in the current version of the model. All H2S and NH3 degassing from the ocean to the atmosphere is assumed to be completely oxidized by O2 to SO42- and NO3- and returns to the ocean surface. These simplifications limit application of the model to very poorly oxygenated Earth system states (pO2<10-5 PAL). Ocean model performance was tested for the modern-day ocean field observational data (Sect. 3). Simulation results were also compared to previously published integrated global flux estimates.

The CANOPS model has been extended and altered a number of times since first publication. The description of biogeochemical cycles in the original version of CANOPS (Ozaki and Tajika, 2013; Ozaki et al., 2011) does not include the S and CH4 cycles because of their aims to investigate the conditions for the development of oceanic anoxia and euxinia on timescales less than a million years during the Phanerozoic. More recently, Ozaki et al. (2019a) implemented an open-system modeling approach for the global S and CH4 cycles, enabling quantitative analysis of global redox budget for given atmospheric O2 levels and crustal reservoir sizes. In this version of CANOPS atmospheric O2 levels and sedimentary reservoirs are treated as boundary conditions because imposing them simplifies the model and significantly reduces computing time. However, this approach does not allow exploration of the dynamic behavior of atmospheric O2 in response to other boundary conditions. In the newest version presented here, significant improvements in the representation of global biogeochemistry were achieved by (1) an explicit calculation of atmospheric O2 levels based on atmospheric mass balance (Sect. 2.6), (2) expansion of the model framework to include secular evolution of sedimentary reservoirs (Sect. 2.5.5), and (3) simplification of the global redox budget between the surface (ocean–atmosphere–crust) system and the mantle (Sect. 2.3.5). These improvements are in line with the requirement of an “open” Earth system model, which is necessary for a systematic, quantitative understanding of Earth's oxygenation history.

Table 1Biogeochemical reactions considered in the CANOPS-GRB model.

a OM denotes organic matter, (CH2O)α(NH4+)βH3PO4. b ΣH2S = H2S + HS.

Download Print Version | Download XLSX

2.3 Global biogeochemical cycles

We construct a comprehensive biogeochemical model in order to investigate the interaction between dynamic behaviors of Earth's oxygenation history and its biogeochemical processes, as well as redox structure of the ocean. Here we provide the basic implementation of global biogeochemical cycles of C, P, N, and S, with particular emphasis on processes of mass exchange between reservoirs that play a critical role in global redox budget (Fig. 2). Our central aim here is to describe the overall design of the biogeochemical cycles. The details of each sub-model are provided in the following sections.

Figure 2Schematics of global biogeochemical cycling in CANOPS-GRB. (a) Global C cycle. The primary source of C for the ocean–atmosphere system is volcanic degassing and oxidative weathering of sedimentary organic carbon, whereas primary sink is burial of marine and terrigenous organic matter into sediments. Inorganic carbon reservoirs (depicted as dashed boxes) and DOC are not considered. NPPocn represents marine net primary production. DIC stands for dissolved inorganic carbon. POC stands for particulate organic carbon. MSR stands for microbial sulfate reduction. AOM stands for anaerobic oxidation of methane. CANOPS-GRB includes CH4 generation via methanogenesis and its oxidation reactions via methanotrophy and AOM in the ocean interior, as well as CH4 degassing flux to the atmosphere and its photooxidation. The rates of CH4 photooxidation and hydrogen escape to space are calculated based on parameterizations proposed by previous studies (Goldblatt et al., 2006; Claire et al., 2006). Note that CH4 flux from land biosphere is not shown here. (b) Global P cycle schematic. Weathering of reactive P (Preac) is the ultimate source, whereas burial in sediments is the primary sink. A part of the weathered P is buried as terrigenous organic P, and the remaining amount is delivered to the ocean. The redox-dependent P burial in marine sediments is modeled by considering three phases (organic P, Fe-bound P, and authigenic P). DIP stands for dissolved inorganic phosphorus. POP stands for particulate organic phosphorus. The hypothetical P scavenging via Fe species in anoxic-ferruginous waters is depicted, but it is not modeled in our standard model configuration. (c) Global N cycle schematic. The abundance of inorganic nitrogen species (ammonium and nitrate), which are lumped into DIN (dissolved inorganic nitrogen), is affected by denitrification and nitrification. The primary source is nitrogen fixation and riverine flux, whereas primary sink is denitrification and burial in marine sediments. PON stands for particulate organic nitrogen. The nitrogen weathering and riverine flux is assumed to be equal to the burial flux so that there is no mass imbalance in global N budget. Aeolian delivery of N from continent to the ocean is not included. (d) Global S cycle schematic. The reservoir sizes of sedimentary sulfur (pyrite sulfur, PYR, and gypsum sulfur, GYP) and two sulfur species (SO42- and ΣH2S) in the ocean are controlled by volcanic outgassing, weathering, burial, MSR, AOM, and sulfide oxidation reactions. Weathering and volcanic inputs are the primary source of S to the ocean, and burial of pyrite and gypsum in marine sediments is the primary sink. It is assumed that hydrogen sulfide escaping from the ocean to the atmosphere is completely oxidized and returns to the ocean as sulfate. The organic sulfur cycle is ignored in this study.

2.3.1 Carbon cycle

The CANOPS-GRB model includes particulate organic carbon (POC), atmospheric CH4, dissolved CH4 in the ocean, and sedimentary organic carbon (ORG) as carbon reservoirs (Fig. 2a). The primary sources of carbon for the ocean–atmosphere system are volcanic degassing and oxidative weathering of sedimentary organic carbon, while the primary sink is burial of marine and terrigenous organic matter in sediments. Atmospheric CO2, dissolved inorganic carbon (DIC), and dissolved organic carbon (DOC) are not explicitly modeled in the current version of the model, and the full coupling of the inorganic carbon cycle within CANOPS-GRB is left as an important topic for future work. Neglecting the inorganic carbon cycle means that there are no climatic feedbacks in the system, and because of this simplification, the CANOPS-GRB model cannot be applied to problems such as those in which the Earth's climate and redox states of the ocean–atmosphere system are closely related each other or to validate model predictions based on geologic records (such as δ13C), but this allows us to avoid introducing additional complexities and uncertainties into the model.

Organic carbon cycle

The biogeochemical model is driven by the cycling of the primary nutrient phosphorus, which is assumed to be the ultimate limiting factor for biological productivity (see Sect. 2.4.3). Previous versions of CANOPS do not take into account the impact of the activity of terrestrial ecosystem on the global O2 budget. In the CANOPS-GRB model, we improve on this by evaluating the activity levels of terrestrial and marine ecosystems separately. The global net primary production (NPP), JNPP (in terms of organic C), is given as a sum of the oceanic (JNPPocn) and terrestrial (JNPPlnd) NPP:

(1) J NPP = J NPP ocn + J NPP lnd .

Biological production in the ocean surface layer depends on P availability, while nutrient assimilation efficiency is assumed to be lower in the high-latitude region (Sect. 2.4.3). Terrestrial NPP is affected by the atmospheric O2 level (Sect. 2.5.1). In this study, the flux (in terms of moles per year) is expressed with a capital J, whereas the flux density (in terms of moles per square meter per year) is expressed with a lowercase j.

In our standard model configuration, oceanic primary production follows canonical Redfield stoichiometry (C:N:P=106:16:1) (Redfield et al., 1963). Flexible C:N:P stoichiometry of particulate organic matter (POM) can be explored by changing a user flag. Nutrients (P and N) are removed from seawater in the photic zone via biological uptake and exported as POM to deeper aphotic layers. The exported POM sinks through the water column with a speed of vPOM (the reference value is 100 m d−1). As it settles through the water column, POM is subject to decomposition via a series of respiration pathways dependent on the redox state of proximal seawater (Sect. 2.4.3). This gives rise to the release of dissolved constituent species back into seawater. Within each layer a fraction of POM is also intercepted by a sediment layer at the bottom of each water depth. Fractional coverage of every ocean layer by seafloor is calculated based on the prescribed bathymetry (Sect. 2.4.1). Settling POM reaching the seafloor undergoes diagenetic alteration (releasing additional dissolved species into seawater) and/or permanent burial. The ocean model has 2×60 sediment segments (HD and LD have 60 layers, respectively), and for each segment the rates of organic matter decomposition and burial are calculated by semi-empirical relationships extracted from ocean sediment data and 1-D modeling of early diagenesis (Sect. 2.4.4). Specifically, the organic C (Corg) burial at each water depth is calculated based on the burial efficiency (BEorg), which is defined as the fraction of POC buried in sediments relative to that deposited on the seafloor at each water depth and is also a function of sedimentation rate and bottom-water O2 levels. Organic matter not buried is subject to decomposition.

The key biogeochemical fluxes of our reference state (mimicking the present condition) are summarized in Table 2. The reference value for burial rate of terrigenous Corg is set at 3 Tmol C yr−1, assuming that burial of terrigenous organic matter accounts for ∼20 % of the total burial. Combined with the burial rate of marine Corg in our standard run, the total burial rate is 14.3 Tmol C yr−1, representing the dominant O2 source flux to the modern ocean–atmosphere system. At steady state, this is balanced by oxidative weathering and volcanic outgassing of sedimentary Corg: The reference value of oxidative weathering of organic matter is determined as 13.0 Tmol C yr−1 based on the global O2 budget (Sect. 2.3.5). Previous versions of CANOPS (Ozaki et al., 2019a) treat sedimentary reservoirs as a boundary condition. This model limitation is removed in the CANOPS-GRB model – the reservoir size of sedimentary Corg (ORG) freely evolves based on the mass balance through burial, weathering, and volcanic outgassing (Sect. 2.5.5). We adopted an oft-quoted value of 1250 Emol (E=1018) for our reference value of the ORG, based on literature survey (Berner, 1989; Garrels and Lerman, 1981).

Table 2Key biogeochemical fluxes obtained from the reference run. * denotes the reference value. Tmol = 1012 mol, Tg = 1012 g.

Download Print Version | Download XLSX

Methane cycle

The ocean model includes biogenic CH4 generation via methanogenesis and its oxidation reactions via methanotrophy and anaerobic oxidation of methane (AOM) in the ocean interior (Reactions R10 and R11 in Table 1), as well as a CH4 degassing flux to the atmosphere. The land model also calculates the biogenic CH4 flux from the terrestrial ecosystem to the atmosphere using a transfer function (Sect. 2.5.2). The abundance of CH4 in the atmosphere is explicitly modeled as a balance of its source (degassing from marine and terrestrial ecosystems) and sink (photooxidation and hydrogen escape), where CH4 sink fluxes are calculated according to parameterized O2-dependent functions proposed by previous studies. More specifically, the oxidation rate of CH4 in the upper atmosphere is calculated based on the empirical parameterization obtained from a 1-D photochemistry model (Claire et al., 2006). The rate of hydrogen escape to space is evaluated with the assumption that it is diffusion limited and that CH4 is a major H-containing chemical compound carrying hydrogen to the upper atmosphere (Goldblatt et al., 2006). No continental abiotic or thermogenic CH4 fluxes are taken into account because previous estimates of the modern fluxes are negligible relative to the biogenic flux, although we realize that it could have played a role in the global redox budget (<0.3 Tmol yr−1; Fiebig et al., 2009). We also note that the current version of the model does not include the possibility of aerobic CH4 production in the sea (Karl et al., 2008). Our reference run calculates atmospheric CH4 to be 0.16 ppmv (Fig. 17), slightly lower than that of the preindustrial level of 0.7 ppmv (Raynaud et al., 1993; Etheridge et al., 1998), but we consider this to be within reasonable error given unknowns in the CH4 cycle.

2.3.2 Phosphorus cycle

Phosphorus is an essential element for all life on Earth and it is regarded as the “ultimate” bio-limiting nutrient for primary productivity on geologic timescales (Tyrrell, 1999). Thus, the P cycle plays a prominent role in regulating global O2 levels. In the CANOPS-GRB model, we model the reactive (i.e., bioavailable) P (Preac) cycling in the system and ignore non-bioavailable P. Specifically, dissolved inorganic P (DIP) and particulate organic P (POP) are explicitly modeled (Fig. 2b), whereas dissolved organic P (DOP) is ignored.

On geologic timescales, the primary source of P to the ocean–atmosphere system is continental weathering: a phosphorus is released through the dissolution of apatite which exists as a trace mineral in silicate and carbonate rocks (∼0.1 wt %; Föllmi, 1996). The total Preac flux via weathering, JPw, is given as follows:

(2) J P w = f P f R J P w , ,

where * denotes the reference value and fP and fR are parameters that control the availability of P in the system. Specifically, fR is a global erosion factor representing the impact of tectonic activity on total terrestrial weathering rate, and fP represents the availability of Preac, which is used in a sensitivity experiment to assess the response of atmospheric O2 levels to changing Preac availability (Sect. 4.1). A fraction of the weathering flux JPw is removed via burial on land, while the remainder is transported to the ocean:


where JPb,lnd and JPr denote the burial rate of terrigenous organic P and riverine Preac flux to the ocean, respectively, k11 is a reference value for the fraction of the total P flux removed by the terrestrial biosphere, and V denotes the vegetation mass normalized to the modern value. These treatments are based on the Earth system box model COPSE (Bergman et al., 2004; Lenton et al., 2016, 2018; Lenton and Watson, 2000b), which has been extensively tested and validated against geologic records during the Phanerozoic. In the CANOPS-GRB model, JPr is tuned so that modeled oceanic P inventory of the reference state is consistent with modern observations of the global ocean (Sect. 3.2.3). Our resulting tuned value is 0.155 Tmol P yr−1, falling in the mid-range of published estimates of 0.11–0.33 Tmol P yr−1, although previous estimates of the riverine Preac flux show large uncertainty (Sect. 3.2.3).

Note that our representation of P weathering ignores the effect of climate (Eq. 2). In the current version of the model the rate of P weathering is treated as one of the model forcings. Although ignoring the climate feedback on P mobility makes interpretation of the model results more straightforward, the incorporation of a climate-sensitive crustal P cycle is an important avenue for future work.

Since atmospheric P inputs are equivalent to less than 10 % of the continental P supply to the modern oceans and much of this flux is not bioavailable (Graham and Duce, 1979), we neglect the aeolian flux in this study. Therefore, riverine input is the primary source of Preac to the ocean. We highlight that open-system modeling is crucial for realistic simulations of ocean biogeochemistry on timescales longer than the residence time of P in the ocean (15–20 kyr for the modern ocean) (Hotinski et al., 2000), and in this framework the riverine input of Preac must be balanced over the long term by loss to sediments via burial. The change in total marine Preac inventory, MP, is given as follows:

(5) d M P d t = J P r - J P b , ocn ,

where JPb,ocn denotes the total burial flux of Preac in the marine system, which is the sum of the burial fluxes of three reactive phases, i.e., organic P (Porg), Fe-bound P (P-Fe), and Ca-bound P (P-Ca) (Sect. 2.4.4):

(6) J P b , ocn = J Porg b + J P - Fe b + J P - Ca b .

O2-dependent P burial is taken into account using empirical relationships from previous studies (Slomp and Van Cappellen, 2007; Van Cappellen and Ingall, 1996, 1994). The burial of Porg at each water depth is a function of burial efficiency, which is controlled by the burial efficiency of organic matter, C/P stoichiometry of POM, sedimentation rate, and bottom-water [O2]. We note that the strength of anoxia-induced P recycling in marine sediments is very poorly constrained, especially in the Precambrian oceans (Reinhard et al., 2017b). Recent studies also suggest that the P retention potential in marine sediments could be affected not only by bottom-water O2 levels but also by redox states (sulfidic vs. ferruginous) and the Ca2+ concentration of bottom waters, as well as various environmental conditions such as temperature and pH (Zhao et al., 2020; Algeo and Ingall, 2007; Papadomanolaki et al., 2022). These are fruitful topics for future research.

We do not explicitly account for P removal via hydrothermal processes because it is estimated that this contribution is secondary in the modern marine P cycle (0.014–0.036 Tmol P yr−1; Wheat et al., 1996, 2003). We note, however, that the hydrothermal contribution to the total P budget in the geologic past remains poorly constrained. We also note that in anoxic, ferruginous oceans, P scavenging by Fe minerals could also play an important role in controlling P availability and the overall budget (Reinhard et al., 2017b; Derry, 2015; Laakso and Schrag, 2014). Modern observations (Dellwig et al., 2010; Turnewitsch and Pohl, 2010; Shaffer, 1986) and modeling efforts (Yakushev et al., 2007) of the redoxcline in the Baltic Sea and the Black Sea suggest an intimate relationship between Mn, Fe, and P cycling. Trapping efficiencies of DIP by settling authigenic Fe and Mn-rich particles were found to be as high as 0.63 (the trapping efficiency is defined as the downward flux of P in Mn and Fe oxides divided by the upward flux of DIP) (Turnewitsch and Pohl, 2010). Although coupled Mn–Fe–P dynamics might have been a key aspect of the biogeochemical dynamics in the Precambrian oceans, we exclude this process in our standard model due to poor constraints and provide a clear and simplified picture of basic model behavior. The key features between the P availability and atmospheric O2 levels are explored by changing fP in this study (Sect. 4).

2.3.3 Nitrogen cycle

In the CANOPS-GRB model, two dissolved inorganic nitrogen (DIN) species (total ammonium ΣNH4+ and nitrate NO3-) and particulate organic nitrogen (PON) are explicitly calculated (Fig. 2c). Atmospheric nitrogen gas is assumed to never limit biospheric carbon fixation and is not explicitly calculated. Dissolved organic N (DON) and terrestrial N cycling (e.g., N fixation by terrestrial ecosystems and riverine-terrigenous organic N transfer) are ignored.

In the surface ocean, N assimilation via nitrate and ammonium depends on the availability of these compounds. If the N required for sustaining a given level of biological productivity is not available, the additional N required is assumed to be provided by atmospheric N2 via nitrogen fixers. The ocean model explicitly calculates denitrification and nitrification reactions in the water column and marine sediments (Reactions R5 and R8 in Table 1). The benthic denitrification rate is estimated using a semi-empirical parameterized function obtained from a 1-D early diagenetic model (see Sect. 2.4.4), while nitrification is modeled as a single step reaction (Reaction R8). N2O and its related reactions, such as anammox, are not currently included.

The oceanic N cycle is open to external inputs of nitrogen. While the ultimate source of N to the ocean–atmosphere system is weathering of organic N, nitrogen fixation represents the major input flux to the ocean with the capacity to compensate for N loss due to denitrification. The time evolution of DIN inventory, MN, in the ocean can be written as follows:

(7) d M N d t = J N fix - J deni wc - J deni sed + J N org w - J Norg b ,

where JNfix denotes the N fixation rate and Jdeniwc and Jdenised are denitrification rates in the water column and sediments, respectively. The first set of terms on the right-hand side represents the internal N cycle in the ocean–atmosphere system, while the second set of terms represents the long-term N budget which interacts with sedimentary reservoir. Ultimately, loss of fixed N from the ocean–atmosphere system only occurs via burial of organic N (Norg) in sediments, JNorgb. This loss is compensated for by continental weathering, JNorgw, which is assumed to be equal to the burial rate of Norg so that the N cycle has no impact on the global redox budget. In the current version of the model, we ignore aeolian flux and all riverine N fluxes other than weathering since these are minor relative to N fixation (Wang et al., 2019). As a result, modeled N fixation required for oceanic N balance can be regarded as an upper estimate.

2.3.4 Sulfur cycle

The original CANOPS ocean model (Ozaki and Tajika, 2013; Ozaki et al., 2011) treated two sulfur species, SO42- and ΣH2S, in a closed system. Neither inputs to the ocean from rivers, hydrothermal vents, and submarine volcanoes nor outputs due to evaporite formation and sedimentary pyrite burial were simulated. This simplification can be justified when the timescale of interest is less than the residence time of the S cycle (∼10–20 Myr). The recently revised CANOPS model (Ozaki et al., 2019a) extends the framework by incorporating the S budget in the ocean. In their model framework, the sedimentary S reservoirs are treated as boundary conditions. The size of sedimentary gypsum and pyrite reservoirs are prescribed, and no explicit calculations of mass balance are performed. In CANOPS-GRB, we removed this model limitation, and the sedimentary reservoirs are explicitly evaluated based on mass balance which is controlled by burial, outgassing and weathering (see Sect. 2.5.5). Specifically, seawater SO42, ΣH2S, and sedimentary sulfur reservoirs of pyrite sulfur (PYR) and gypsum sulfur (GYP) are explicitly evaluated in the current version of the model. No atmospheric sulfur species are calculated – all H2S degassing from the ocean to the atmosphere is assumed to be oxidized to sulfate and to return to the ocean. The organic sulfur cycle is not considered in this study.

Sulfur enters the ocean mainly from river runoff, JSr, with minor contributions from volcanic outgassing of sedimentary pyrite, Jpyrm, and gypsum, Jgypm. The reference value for the riverine flux is set at 2.6 Tmol S yr−1, consistent with the published estimate of 2.6±0.6 Tmol S yr−1 (Raiswell and Canfield, 2012). The riverine flux is written as the sum of gypsum weathering and oxidative weathering of pyrite: JSr=Jgypw+Jpyrw. Sulfur weathering fluxes are also assumed to be proportional to the sedimentary reservoir size. Estimates of modern volcanic input fall within the range of 0.3–3 Tmol S yr−1 (Catling and Kasting, 2017; Kagoshima et al., 2015; Raiswell and Canfield, 2012; Walker and Brimblecombe, 1985). We adopted a value of 0.8 Tmol S yr−1 for this flux (Kagoshima et al., 2015). Our total input of 3.4 Tmol S yr−1 is also within the range of the previous estimate of 3.3±0.7 Tmol S yr−1 (Raiswell and Canfield, 2012). Sulfur is removed from the ocean either via pyrite burial, Jpyrb, or gypsum deposition, Jgypb (Fig. 2d). The time evolution of the inventory of total S in the ocean can thus be written as follows:

(8) d M SO 4 + M H 2 S d t = J S r + J pyr m + J gyp m - J pyr b + J gyp b ,

where MSO4 and MH2S denote the inventory of sulfate and hydrogen sulfide in the ocean, respectively. Two sulfur species (SO42- and ΣH2S) are transformed via microbial sulfate reduction (MSR) (Reaction R6), AOM (Reaction R11), and aerobic sulfide oxidation reactions (Reaction R9). The above equation thus can be divided into following equations:


where JH2Sox denotes the oxidation of hydrogen sulfide and JMSR&AOM is sulfate reduction via MSR and AOM. Pyrite burial is represented as the sum of pyrite precipitation in the water column and sediments: Jpyrb=Jpyrb,wc+Jpyrb,sed, where the pyrite burial rate in marine sediments is assumed to be proportional to the rate of benthic sulfide production. The proportional coefficient, pyrite burial efficiency (epyr), is one of the tunable constants of the model. For normal (oxic) marine sediments, epyr is tuned such that the seawater SO42- concentration for our reference run is consistent with modern observations (Sect. 2.4.4). Pyrite precipitation in the water column is assumed to be proportional to the concentration of ΣH2S.

Although the present-day marine S budget is likely out of balance because of a lack of major gypsum formation, the S cycle can be considered to operate at steady state on timescales longer than the residence time of sulfur in the ocean. According to S isotope mass balance calculations, ∼10 %–45 % of the removal flux is accounted for by pyrite burial, and the remainder is removed via formation of gypsum and anhydrite in the near-modern oceans (Tostevin et al., 2014). Although gypsum deposition would have been strongly influenced by tectonic activity (Halevy et al., 2012), we assume that the rate of gypsum deposition on geologic timescales is proportional to the ion product of Ca2+ and SO42- (Berner, 2004b) in the low- to mid-latitude surface layer (L) and is defined as follows:

(11) J gyp b = [ Ca 2 + ] l [ SO 4 2 - ] l [ Ca 2 + ] [ SO 4 2 - ] J gyp b , = f Ca [ SO 4 2 - ] l [ SO 4 2 - ] J gyp b , ,

where “l” denotes the low- to mid-latitude surface layer and fCa is a parameter that represents the seawater Ca2+ concentration normalized by the present value (fCa=1 for the reference run). The reference value of gypsum burial Jgypb, is determined by assuming that gypsum deposition accounts for ∼60 % of the total S removal from the near-modern ocean.

2.3.5 Global redox budget

In the previous version of the CANOPS model (Ozaki et al., 2019a), the atmospheric O2 level was prescribed as a boundary condition, rather than modeled in order to limit computational demands. In this study, we remove this model limitation by introducing an explicit mass balance calculation of atmospheric O2 (Sect. 2.6.3). This improvement allows us to explore the dynamic response of O2 levels in the ocean-atmosphere system (Sect. 4).

Figure 3Schematics of global redox (O2) budget. Arrows represent the O2 flux. The primary source is burial of organic carbon and pyrite sulfur in sediments and hydrogen escape to space. The primary sink is volcanic outgassing and weathering of crustal organic matter and pyrite. PYR is the sedimentary reservoir of pyrite sulfur. ORG is the sedimentary reservoir of organic carbon. CANOPS-GRB tracks the global redox (O2) budget for each simulation.


CANOPS-GRB v1.0 is designed to be a part of a comprehensive global redox budget (GRB) framework (Fig. 3) (Catling and Kasting, 2017; Ozaki and Reinhard, 2021). Here GRB is defined for the combined ocean–atmosphere system. In this study we track GRB in terms of O2 equivalents. The ultimate source of O2 is the activity of oxygenic photosynthesis (and subsequent burial of reduced species, such as organic matter and pyrite sulfur, in sediments), whereas the primary sink of O2 is the oxidative weathering of organic carbon and pyrite which are assumed to be O2 dependent (Sect. 2.5.3). On timescales longer than the residence time of O2 in the ocean–atmosphere system, O2 source fluxes should be balanced by sink fluxes. Specifically, the O2 budget in the coupled ocean–atmosphere system can be expressed as follows:

(12) GRB = J org b , ocn + J org b , lnd - J org w - J org m + 2 J pyr b - J pyr w - J pyr m + J Hesc - J man ,

where the first and second set of terms on the right-hand side represent the redox balance via organic carbon and pyrite sulfur subcycles, respectively. JHesc in the third term denotes hydrogen escape to space, representing the irreversible oxidation of the system. For well-oxygenated atmospheres this process plays a minor role in the redox budget, but for less-oxygenated atmospheres with high levels of CH4 this flux could lead to redox imbalance. In this study we include the input of reducing power (e.g., H2 and CO) from the Earth's interior to the surface, Jman, which is assumed to be equal to the value of JHesc (Jman=JHesc) to avoid redox imbalance in the exogenic system. In reality, mantle degassing and the rate of hydrogen escape are not necessarily equal, resulting in redox imbalance that may exert a fundamental control on atmospheric redox chemistry on geologic timescales (Hayes and Waldbauer, 2006; Ozaki and Reinhard, 2021; Canfield, 2004; Eguchi et al., 2020); however, to maintain simplicity we have left this as a topic for future work. As a result, the terms on the right-hand side must be balanced at steady state. Our model can meet this criterion. Note that the effects of the Fe cycle on the O2 budget (e.g., the oxidative weathering of Fe(II)-bearing minerals; Ozaki et al., 2019a) are not included in the core version of the CANOPS-GRB v1.0 code and in the analyses presented here for the sake of simplicity.

The CANOPS-GRB model also tracks the O2 budgets for the atmosphere (ARB) and ocean (ORB) independently:


where Φexair-sea represents the net exchange of oxidizing power between the ocean and atmosphere via gas exchange (O2 with minor contributions of NH3, H2S and CH4). These separate redox budgets are also tracked in order to validate global budget calculations.

For our reference condition, we obtain the reference value for the oxidative weathering rate of Corg (Jorgw,) using the redox budget via Corg subcycle:

(15) J org w , = J org b , ocn , + J org b , lnd , - J org m , .

Given flux values based on the calculated (Jorgb,ocn,=11.28 Tmol C yr−1) and prescribed (Jorgb,lnd,=3 Tmol C yr−1, Jorgm,=1.25 Tmol C yr−1) values on the right-hand side, Jorgw, is estimated as 13.03 Tmol C yr−1 (Table 2).

2.4 Ocean model

Here we undertake a thorough review, reconsideration, and revision (where warranted) of all aspects of the ocean model, including bringing together developments of the model following the original papers describing the CANOPS ocean model framework (Ozaki and Tajika, 2013; Ozaki et al., 2011).

The ocean model includes exchange of chemical species with external systems via several processes such as air–sea exchange, riverine input, and sediment burial. The biogeochemical model also includes a series of biogeochemical processes, such as the ocean biological pump and redox reactions under oxic, anoxic, and sulfidic conditions. Our ocean model is convenient for investigating Earth system changes on timescales of hundreds of years or longer and can be relatively easily integrated, rendering the model unique in terms of biogeochemical cycle models. CANOPS is also well suited for sensitivity studies and can be used to obtain useful information upstream of more complex models.

Development of the ocean model included two initial goals: First, to adopt a general and robust ocean circulation scheme capable of producing well-resolved modern distributions of circulation tracers using realistic ventilation rates with a limited number of free parameters. The model outputs for circulation tracers are validated by comparison with modern observations (see Sect. 3). This confirms that our ocean circulation scheme is adequate for representing the global patterns of water mass transport. A second key goal was to couple the circulation model with an ocean biogeochemical model and to evaluate performance by comparison with modern ocean biogeochemical data (see Sect. 3.2). Examination of the distributions and globally integrated fluxes of C, N, P, S, and O2 for the modern ocean reveals that the ocean model can capture the fundamentals of marine biogeochemical cycling.

2.4.1 Structure

CANOPS ocean model is a 1-D (vertically resolved) intermediate-complexity box model of ocean biogeochemistry (see Fig. 1b for the schematic structure) originally developed by Ozaki and Tajika (2013) and Ozaki et al. (2011). Our model structure is an improved version of the HILDA model (Joos et al., 1991; Shaffer and Sarmiento, 1995). Unlike simple one-dimensional global ocean models (e.g., Southam et al., 1982), the HILDA-type model includes explicit high-latitude dynamics whereby the high-latitude surface layer exchanges properties with the deep ocean. This treatment is crucial for simulating preformed properties and observed chemical distributions, especially for phosphate and dissolved O2 in a self-consistent manner. Unlike simple box-type global ocean models, the model has high vertical resolution. This is needed for representing proper biogeochemical processes that show strong depth dependency. Furthermore, HILDA type models (Arndt et al., 2011; Shaffer et al., 2008), unlike multi-box-type global ocean models (Hotinski et al., 2000), use a small number of free parameters to represent ocean physics and biology. The simple and adaptable structure of the model should make it applicable to a wide range of paleoceanographic problems. The ocean model couples a diffusion-advection model of the global ocean surface and interior with a biogeochemical model (Sect. 2.4.3) and a parameterized sediment model (Sect. 2.4.4).

The ocean surface consists of a mixed layer at low- to mid-latitude (L) and high-latitude (H). Below the surface layers, we adopt the present-day averaged seafloor topography of Millero (2006), with the hypsometric profile shown in Fig. 4a. Below the surface water layers, the ocean interior comprises two regions: the high- to mid-latitude region (HD) and the low- to mid-latitude region (LD). Each region is subdivided vertically with high resolution (Δz=100 m). Each of the 60 ocean layers in each latitude region (120 total) is assigned ocean sediment properties. The cross-sectional area, volume, and sediment surface area of each box is calculated from the benthic hypsometry. Inclusion of the bathymetry allows evaluation of the flux of biogenic materials which settle on, and are buried in, seafloor sediments at each water depth (Sect. 2.4.3 and 2.4.4).

Figure 4Ocean bathymetry and water transport. (a) Seafloor topography (cumulative seafloor area) (Millero, 2006) adopted in the CANOS-GRB model. (b) Lateral water advection from HD to LD section assumed in the standard run (in Sv). Total advection rate V˙ was set at 20 Sv. (c) Upwelling rate in the LD region (in m yr−1) of the standard run.

2.4.2 Transport

The ocean circulation model represents a general and robust scheme that is capable of producing well-resolved modern profiles of circulation tracers using realistic parameter values and the coupled biogeochemical model (Sect. 2.4.3) or parameterized sediment model (Sect. 2.4.4).

The time–space evolution of model variables in the ocean is described by a system of horizontally integrated vertical diffusion equations for non-conservative substances. The tracer conservation equation establishes the relationship between change of tracer concentration at a given grid point and the processes that can change that concentration. These processes include water transport by advection and mixing and sources and sinks due to biological and chemical transformations. The temporal and spatial evolution of the concentration of a dissolved component in the aphotic zone is described by a horizontally integrated vertical diffusion equation, which relates the rate of change of tracer concentration at a given point to the processes that act to change the tracer concentration:

(16) [ X ] t = [ X ] t trans + Θ bio + Θ react ,

where [X] represents horizontally integrated physical variables (such as potential temperature, salinity, or 14C) or concentration of a chemical component, t denotes time, and Θbio and Θreact represent internal sources and sinks associated with the biological pump and chemical reactions, respectively. An external source or sink term Θex, which represents riverine input and/or air–sea gas exchange, is added to the surface layers. The first term on the right-hand side of Eq. (16) represents the physical transport:

(17) [ X ] t trans = - A l , h ( z ) w l , h ( z ) [ X ] z + z A l , h ( z ) K v l , h ( z ) [ X ] z + K hor 2 [ X ] y 2 .

The terms on the right-hand side express (from left to right) the advection, vertical diffusion, and horizontal diffusion. Here, l and h indicate the LD and HD, respectively. The factors Kvl,h(z), Khor, Al,h(z), and wl,h(z) denote the vertical and horizontal diffusion coefficients, the areal fraction of the water layer at water depth z to the sea surface area, and upwelling (for LD) or downwelling (for HD) velocity, respectively.

In the CANOPS ocean model, ocean circulation and mixing are characterized by five physical parameters: (1) water transport via thermohaline circulation, V˙, associated with high-latitude sinking and low- to mid-latitude upwelling, (2) constant horizontal diffusion between the aphotic zones, Khor, (3) strong, depth-dependent vertical diffusion between the aphotic zones in the high-latitude region, Kvh(z), (4) high-latitude convection, V˙h, and (5) depth-dependent vertical diffusion in the low- to mid-latitude region, Kvl(z). These parameters are tuned to give tracer distributions consistent with present-day observations. Thermohaline circulation and high-latitude convection are considered to be general physical modes on any rotating planet, and our simplified water transport scheme allows us to represent them with a limited number of free parameters. However, we emphasize that the water transport scheme explored here is designed to represent the modern ocean circulation on Earth. As a result, some of these parameterizations may need to be modified when being applied to ancient oceans or oceans on exoplanets. Nevertheless, given our simple design, our water transport scheme is relatively flexible to modify the water circulations that are markedly different from the modern ocean.


Advective water transport in the ocean model represents the major features of modern meridional overturning circulation. The rate of production of ventilated ocean waters ranges from 14 to 27 Sv (1 Sv = 106 m3 s−1) in the North Atlantic and from 18 to 30 Sv in the Southern Ocean (e.g., Doney et al., 2004; Lumpkin and Speer, 2007). The formation of deepwater effectively supplies “fresh” ventilated water to the abyss. We choose V˙=20 Sv as a reference value, giving a mean overturning time of about 2140 years, consistent with the ventilation time estimated from observations (Broecker and Peng, 1982).

The downwelling of the surface waters at H forms HD that flows into the intermediate to deep oceanic layers of LD, which in turn upwells over L (Fig. 1b). In many one-dimensional ocean models, downwelling water enters the ocean interior via the deepest model layer (e.g., Southam et al., 1982; Shaffer and Sarmiento, 1995; Volk and Hoffert, 1985). In the real ocean, downwelling waters are transported along isopycnal layers below approximately 1000 m (e.g., Doney et al., 2004; Lumpkin and Speer, 2007; Shaffer and Sarmiento, 1995; Volk and Hoffert, 1985). Hence, we assume that high-latitude deep water flows into each ocean layer below 1100 m. While there is some uncertainty in the pattern of lateral advection, the flow is determined in our model assuming a constant upwelling rate below a depth of 1100 m in the LD region. The upwelling and downwelling rate wl,h(z) is then determined by the seafloor topography and the deep-water lateral inflow, assuming continuity. Figure 4b shows the lateral advection of deep waters with a reference circulation rate V˙ of 20 Sv. This assumption provides a plausible upwelling rate, which is consistent with the oft-quoted value of 2–3 m yr−1 (Broecker and Peng, 1982) (Fig. 4c).

Vertical mixing

Ocean circulation is dominated by turbulent processes driven by wind and tidal mixing. These processes occur as eddies, which occur at a wide range of spatial scales, from centimeters to whole ocean basins. In numerical models of ocean circulation, turbulent mixing in the ocean interior is commonly represented as a diffusion process, characterized by an eddy diffusion coefficient. The vertical eddy diffusion coefficient Kv(z) is typically on the order of 10−5 to 10−4 m2 s−1, and it is common to assume a depth dependence that smoothly increases from the thermocline (10-5 m2 s−1) to the abyss (10-4 m2 s−1) using an inverse or hyperbolic tangent function (e.g., Shaffer et al., 2008; Yakushev et al., 2007). To account for thermocline ventilation, we assumed a relatively high vertical diffusion coefficient in mid-water depth (Kl=6.3×10-5 m2 s−1 for water depth 500–1500 m). We also adopted a higher value for the vertical diffusion coefficient (Ku=1.6×10-4 m2 s−1) in the uppermost 500 m of the ocean in order to represent the highly convective Ekman layer in the upper part of the ocean.

(18) K v l ( z ) = K u ( z - 500 m ) K l ( - 500 z - 1500 m ) κ s + κ d - κ s 2 1 + tanh z - z l z l ( otherwise ) ,

where κs and κd are vertical mixing coefficients and zl is the transition length scale (Romaniello and Derry, 2010). In the high-latitude region where no permanent thermocline exists, more rapid communication with deep waters can occur. Previous studies have pointed out that the vertical diffusivities at high latitude can be very high (up to O(10−2 m2 s−1)) (e.g., Sloyan, 2005). To account for this we include high-latitude convection between H and YD (V˙h=57.4 Sv) and higher vertical diffusion (Kvh(z)=2×Kvl(z)).

Horizontal diffusion

The horizontal diffusivity is included according to Romaniello and Derry (2010). On basin scales, the horizontal (isopycnal) eddy diffusivity is 107–108 times larger than the vertical (diapycnal) eddy diffusivity due to anisotropy of the density field. For a spatial scale of 1000 km, horizontal eddy diffusion is estimated to be O(103 m2 s−1) (e.g., Ledwell et al., 1998). We adopt this value. As Romaniello and Derry (2010) did previously, we assume horizontal mixing follows the pathways of advective fluxes between laterally adjacent regions. The reciprocal exchange fluxes may be written as

(19) J hor ex = K hor A [ X ] y = K hor A L Δ [ X ] ,

where Jhorex denotes the exchange fluxes between the layers (in mol yr−1), A represents the cross-sectional area separating two adjacent reservoirs, L is a characteristic spatial distance separating the reservoirs, and Δ[X] is the difference in concentration between two reservoirs (Romaniello and Derry, 2010). By assuming that L is of the same order as the length of the interface separating the two regions, we can approximate AΔz×O(L), where Δz is the thickness of the interface separating the two regions. Then we obtain

(20) J hor ex = K hor Δ z Δ [ X ] .

Therefore, when we discretize the ocean interior at 100 m spacing approximately 0.1 Sv of reciprocal mixing occurs between adjacent layers.

Ocean circulation tracers

We use potential temperature θ, salinity S, and radioactive carbon 14C, as physical tracers. Distributions of these tracers are determined by the transport mechanisms described above. In this study, we adopt the values at the surface layers (L and H) as upper boundary conditions: θl=15C, θh=0C, Sl=35 psu, Sh=34 psu, Δ14Cl=-40 ‰, and Δ14Ch=-100 ‰. The radioactive decay rate for 14C is 1.21×10-4 yr−1. Although 14C can be incorporated in the biogenic materials and transported into deep water, we ignore this biological effect for simplicity. The associated error is ∼10 % of the profiles produced by circulation and radioactive decay (Shaffer and Sarmiento, 1995). The parameter values used in the ocean circulation model are listed in Table 3.

Table 3Physical setup of the ocean circulation model.

Download Print Version | Download XLSX

2.4.3 Ocean biogeochemical framework

The ocean circulation model is coupled to a biogeochemical model, which includes an explicit representation of a variety of biogeochemical processes in the ocean. The parameters used in the oceanic biogeochemical model are listed in Table 4.

Table 4Parameter values used in the oceanic biogeochemistry module of CANOPS-GRB.

Download Print Version | Download XLSX

Biological production

The overall biogeochemical cycling scheme is based on the cycling of primary nutrient (phosphate; PO43-), which limits biological productivity – export production is related to the availability of P within the euphotic zone (Maier-Reimer, 1993; Yamanaka and Tajika, 1996; Shaffer et al., 2008):

(21) j exp l , h = α l , h h m ε l , h [ PO 4 3 - ] l , h [ PO 4 3 - ] l , h [ PO 4 3 - ] l , h + K P ,

where jexp represents new and export production of POC (in unit of mol C m−2 yr−1), α denotes C:P stoichiometry of POM, hm is the mixed-layer depth, ε denotes the assimilation efficiency factor for P uptake, and KP denotes the half-saturation constant. The value of ε for the low- to mid-latitude region is assumed to be 1. In contrast, we assume a lower efficiency for the high-latitude region because biological production tends to be limited by environmental factors other than phosphate availability (e.g., amount of solar radiation, mixed-layer thickness, sea ice formation, and iron availability). This is used as one of the fitting parameters in the model. Downwelling waters contain a certain level of nutrients (i.e., preformed nutrients).

In our standard run, the stoichiometry of organic matter is parameterized using the canonical Redfield ratio (C:N:P=106:16:1) (Redfield et al., 1963). However, we note that flexible C:N:P stoichiometry has been the subject of recent discussion. In the modern oceans, C:N:P ratios of exported POM vary across latitude, reflecting ecosystem structure (Galbraith and Martiny, 2015). Local observations (and laboratory experiments) suggest that the C:N:P ratio of cyanobacteria is a function of seawater PO43- concentration (Larsson et al., 2001). The evolutionary perspective has also been discussed (Quigg et al., 2003; Sharoni and Halevy, 2022). In the previous version of the CANOPS model, the C–N–P stoichiometry of primary producers responds dynamically to P availability in the surface layer (Reinhard et al., 2017b):


where α and β represent the C/P ratio and N/P ratio of POM, * denotes the canonical Redfield ratios, max denotes the maximum value (αmax=400 and βmax=60), and γP0 and γP1 are tunable constants (γP0=0.1µM and γP1=0.03µM) (Kuznetsov et al., 2008). In the CANOPS-GRB model, this dynamic response of POM stoichiometry can be explored by changing the user flag from the standard static response. In this study, we do not explore the impacts of flexible POM stoichiometry on global biogeochemistry (i.e., αmax=α and βmax=β).

Biological production in the surface mixed layer increases the concentration of dissolved O2 and reduces the concentrations of DIP and DIN according to the stoichiometric ratio (Reactions R1 and R2; Table 1). DIN consumption is partitioned between nitrate and ammonium, assuming that ammonium is preferentially assimilated. CANOPS-GRB evaluates the availability of fixed N in the surface ocean, and any N deficiency required for a given level of productivity is assumed to be compensated for on geologic timescales by N fixers. In other words, it is assumed that biological N fixation keeps pace with P availability so that P (not N) ultimately determines oceanic biological productivity.

To date, models of varying orders of complexity have been developed to simulate oceanic primary production and nutrient cycling in the euphotic layer, from a single nutrient and single phytoplankton component system to the inclusion of multiple nutrients and trophic levels in the marine ecosystem, usually coupled to physical models (e.g., Yakushev et al., 2007; Oguz et al., 2000). To avoid this level of complexity, we introduce a parameter, fexp, called export ratio (Sarmiento and Gruber, 2006), which relates the flux densities of export production and NPP, as follows:

(24) j NPP ocn = j exp f exp ,

where jNPPocn denotes the NPP in terms of mol C m−2 yr−1. In the modern ocean, the globally averaged value of fexp is estimated at 0.2 (Laws et al., 2000), and we assumed this value in this study. The rate of recycling of organic matter in the photic zone is thus given by

(25) j recy = j NPP ocn - j exp = 1 - f exp f exp j exp .

The respiration pathway of jrecy depends on the availability of terminal electron acceptors (O2, NO3- and SO42-). Following exhaustion of these species as terminal electron acceptors, organic matter remineralization occurs by methanogenesis (Reaction R7). See below for the treatment of organic matter remineralization in the water column.

Biological pump

Most POM exported to the deep sea is remineralized in the water column before reaching the seafloor (e.g., Broecker and Peng, 1982). Nutrients returning to seawater at intermediate depths may rapidly return to the surface ocean and support productivity. The remaining fraction of POM that reaches the sediment ultimately exerts an important control on oceanic inventories of nutrients and O2. An adequate representation of the strength of biological pump is therefore critical to any descriptions of global biogeochemical cycles.

The governing equation of the concentration of biogenic particles G is

(26) G t + v POM G z = - r G ,

where r is a decomposition rate and vPOM is the settling velocity of POM in the water column. We assume a settling velocity of 100 m d−1 for our reference value (e.g., Suess, 1980), although a very wide range of values and depth dependencies have been reported (e.g., Berelson, 2001a). Therefore, the settling velocity is fast enough to neglect advective and diffusive transport of biogenic particles. Note that the settling velocity would affect the intensity of biological pump and chemical distribution in the ocean interior. Considering the ballast hypothesis in the modern ocean (Armstrong et al., 2001; Francois et al., 2002; Ittekkot, 1993; Klaas and Archer, 2002), the settling velocity of POM in the geological past would very likely have been different from the modern ocean. As Kashiyama et al. (2011) pointed out, there would be a critical aspect among sinking rate of POM, intensity of biological pump and chemical distribution in the ocean. The quantitative and comprehensive evaluation of their effect is an important issue for the future work (Fakhraee et al., 2020).

In order to solve Eq. (26) explicitly, a relatively small time step (∼1 d) would be required. However, because the sinking velocity and remineralization of biogenic material are fast processes, we assume that the POM export and remineralization occurs in the same time step (ignoring the term G/t). The concentration of biogenic particles can be solved as follows:

(27) G ( z + Δ z ) = G ( z ) exp - r Δ z v POM ,

where Δz is a spatial resolution of the model.

Organic matter decomposition

As POM settles through the water column, it is nearly entirely decomposed back to dissolved tracers. Therefore, decomposition of POM is a key process for modeling biogeochemistry in the ocean. To avoid the complex treatment of this process (such as repackaging and aggregation or dispersal of particles), various empirical schemes for POM sinking flux have been proposed, such as exponential (Volk and Hoffert, 1985) or power law (Martin et al., 1987) functions (Fig. 5). However, the estimation of Volk and Hoffert generally tends to overestimate in the upper water column (<1.5 km) and underestimate at depth. It is important to note that data series of sediment trap measurements were obtained from a limited geographic and depth range. Berelson (2001b) and Lutz et al. (2002) conducted further estimates of the sediment flux and found regional variability in the sinking flux. Broadly, these data indicate that commonly applied flux relationships generally tend to overestimate flux to depth.

Figure 5Empirical relationships between POC settling flux normalized to export production (Lutz et al., 2002) and water depth (Archer et al., 1998; Berelson, 2001b; Martin et al., 1987; Volk and Hoffert, 1985). The profile of the CANOPS-GRB model is depicted as a red line. The black dots represent observational data (Honjo and Manganini, 1993; Lutz et al., 2002; Tsunogai and Noriki, 1991; Honjo, 1980, and references therein).

The microbial degradation of different groups of organic matter differs over timescales ranging from hours to millions of years. In order to represent the decrease in POM lability with time and water depth, we adopt the so-called multi-G model (Westrich and Berner, 1984) that describes the detailed kinetics of organic matter decomposition (Ozaki and Tajika, 2013; Ozaki et al., 2011). In the CANOPS model, POM is described using two degradable fractions (G1 and G2) and one inert (G3) fraction using different rate constants ki (i=1, 2, 3) for each component. Rate constants are tuned on the basis of consistency with the typical profile of the POM sinking flux estimated from sediment trap studies (Fig. 5). In this study, constant stoichiometries between C, N, and P during the remineralization of POM are assumed throughout the water column, taking values equal to those characterizing mean export production.

The electron acceptor used in the respiration reaction changes from dissolved O2 to other oxidants (e.g., NO3- and SO42-) as O2 becomes depleted. The respiration pathway is controlled by the free energy change per mole of organic carbon oxidized. The organic matter decomposition is performed by the oxidant that yields the greatest free energy change per mole of organic carbon oxidized. When the oxidant is depleted, further decomposition will proceed utilizing the next most efficient (i.e., the most energy producing) oxidant until either all oxidants are consumed or oxidizable organic matter is depleted (e.g., Froelich et al., 1979; Berner, 1989). In oxic waters, organic matter is remineralized by an aerobic oxidation process (Reaction R4). As dissolved O2 is depleted, NO3- and/or SO42- will be used (Reactions R5 and R6). Denitrification is carried out by heterotrophic bacteria under low concentrations of dissolved O2 if there is sufficient nitrate. For anoxic, sulfate-lean oceans, the methanogenic degradation of organic matter will occur (Reaction R7). In the CANOPS-GRB model, we parameterized the dependence of decomposition of POM with a Michaelis–Menten type relationship with respect to the terminal electron acceptors:


where KO2, KNO3, and KMSR are Monod constants and KO2, KNO3, and KMSR are inhibition constants. The Monod-type expressions are widely used in mathematical models of POM decomposition processes (e.g., Boudreau, 1996). The oxidants for organic matter decomposition change with the availability of each oxidant, which vary with time and water depth. The parameter values are based on previous studies on early diagenetic processes in marine sediments (Boudreau, 1996; Van Cappellen and Wang, 1996). SO42- has been one of the major components of the Phanerozoic oceans and has been an important oxidizing agent in anaerobic systems. In the original CANOPS model (Ozaki and Tajika, 2013; Ozaki et al., 2011), it was assumed that the saturation constant KMSR is zero, meaning that the SO42- is never a limiting factor. In contrast, during the Precambrian, seawater SO42- could have been extremely low (Lyons and Gill, 2010). The half-saturation constant for MSR (KMSR) determines the degree to which MSR contributes to the total respiration rates. However, estimates for KMSR in natural environments and pure cultures vary over several orders of magnitude (∼0.002–3 mM) (Tarpgaard et al., 2011; Pallud and Van Cappellen, 2006). We assume a reference value of 0.2 mM for this study.

Finally, temperature may also have played an important role in organic matter decomposition rates. The dependence of ammonification on temperature is sometimes described by an exponential function or Q10 function (e.g., Yakushev et al., 2007). While we recognize that the temperature dependency of organic matter decomposition might have played an important role in oceanic biogeochemical cycles in the geological past (Crichton et al., 2021), these dynamics are not included in CANOPS-GRB v1.0.

Secondary redox reactions

Total ammonia (ΣNH3), total sulfide (ΣH2S), and methane (CH4), produced during organic matter degradation, are subject to oxidation to NO3-, SO42-, and CO2 via a set of secondary redox reactions (Table 1). Rate constants for these reactions are taken from the literature. The ocean model includes nitrification (Reaction R8), total sulfide oxidation by O2 (Reaction R9), aerobic oxidation of CH4 by O2 (Reaction R10), and AOM by SO42- (Reaction R11). Nitrification, the oxidation of ammonium to nitrate, occurs in several stages and is accomplished mainly by chemolithotrophic bacteria (Sarmiento and Gruber, 2006). In this study, we treat all nitrification reactions as a combined reaction (Reaction R8). The rate of this process is assumed to depend on the concentration of both oxygen and ammonia as follows:

(32) R 8 = k R 8 [ NH 4 + ] [ O 2 ] .

The oxidation of sulfide formed in anoxic waters by MSR can also be written as a series of reactions (e.g., Yakushev and Neretin, 1997), but we treat it as an overall reaction (Reaction R9). The rate of this secondary redox reaction is also formulated using a bimolecular rate law:

(33) R 9 = k R 9 [ Σ H 2 S ] [ O 2 ] .

The rate constant for this process has been shown to vary significantly as a function of several redox-sensitive trace metals that act as catalysts (Millero, 1991). Here we assume kR9=3650 mM−1 yr−1 based on the observations of the chemocline of the Black Sea (Oguz et al., 2001).

In the original CANOPS model (Ozaki et al., 2019a; Ozaki and Tajika, 2013), syngenetic pyrite formation in the water column was not considered. In a more recent revision of the model, this process was added (Cole et al., 2022) and parameterized such that iron sulfide formation is assumed to be proportional to the hydrogen sulfide concentration:

(34) R pyr wc = k pyr wc [ Σ H 2 S ] ,

where kpyrwc is a model constant (its reference value is set at 0.01 yr−1). This constant is a function of the ferrous iron concentration in seawater, but it is the subject of large uncertainty. The total flux (in mol S yr−1) can be obtained by integrating the precipitation flux density over the whole ocean:

(35) J pyr wc = R pyr wc d V d z d z .

The aerobic oxidation of CH4 is formulated using a bimolecular rate law:

(36) R 10 = k R 10 [ CH 4 ] [ O 2 ] .

The rate of AOM is formulated using a Monod-type law (Beal et al., 2011):

(37) R 11 = k R 11 [ CH 4 ] [ SO 4 2 - ] K AOM + [ SO 4 2 - ] .

Rate constants for above reactions are taken from the literature (Table 4). Secondary redox reactions were calculated implicitly with an operator splitting scheme (Steefel and Macquarrie, 1996) so as to maintain numerical stability.

2.4.4 Sediment–water exchange

The burial of biogenic material in marine sediments plays a critical role in global biogeochemical cycles, especially with respect to the marine budgets of nutrients, carbon, and sulfur. This is intimately linked to atmospheric O2 levels on geologic timescales. Specifically, the burial rate of Corg in marine sediments exerts a primary control on the evolution of atmospheric O2 levels throughout Earth's history. Given the complexity of biogeochemical processes within sediments and our limited knowledge on many of the early diagenetic processes, we adopt some semi-empirical relationships extracted from ocean sediment data. This approach, rather than explicit modeling, is also required to reduce the computational cost of the simulation on timescales >100 Myr. The related parameter values are listed in Table 5.

Table 5Parameters used in the sediment–water interface module of CANOPS-GRB.

Download Print Version | Download XLSX

POM deposition

The fraction of settling POM that reaches the sediment surface, Jorgdep (in mol C yr−1) is a function of both the settling flux density, jorgdep (in mol C m−2 yr−1), and topography (Fig. 4a):

(38) J org dep = z 1 z 2 j org dep ( z ) d A d z d z ,

where the settling flux density can be written as follows:

(39) j org dep = v POM G ,

where G is the concentration of POM and vPOM denotes the sinking velocity.

Carbon cycling

Interactions between the ocean and underlying sediments play an important role in influencing whole-ocean chemical and nutrient inventories on geologic timescales. POM deposited to the seafloor is subject to decomposition during diagenetic processes associated with burial in marine sediments. Only a small fraction of organic matter will ultimately be buried and removed from the surface environment. However, understanding what factors control the preservation of organic matter in marine sediments has been a controversial topic, and we still lack a robust understanding of this process. With this issue in mind, we adopt an empirical approach obtained using the observational data from previous studies.

The burial flux density of Corg at each water depth, jorgb,ocn (in terms of mol C m−2 yr−1), is calculated based on burial efficiency, BEorg:

(40) j org b , ocn = BE org j org dep .

Burial efficiency is defined as the fraction of organic matter buried in sediments relative to the total depositional flux. Burial efficiency is described by simplified parametric laws based on empirical relationships from modern-day observations. Previous studies demonstrate strong dependency of this term on total sedimentation rate, SR (e.g., Henrichs and Reeburgh, 1987). Figure 6 demonstrates the relationship between BEorg and SR compiled from literature surveys. The sedimentation rate in the modern ocean varies over about 5 orders of magnitude, with a primary dependence on material supplied from the continents. There is a strong relationship, especially for SR less than 0.01 cm yr−1. In contrast to the strong SR dependence under oxic conditions, anoxic settings show a much weaker dependence of BEorg on SR (Betts and Holland, 1991; Henrichs and Reeburgh, 1987) (Fig. 6). In this study, the following relationship proposed by Henrichs and Reeburgh (1987) is adopted for sediments underlying well-oxygenated bottom water (O2 concentration of bottom water, [O2]bw>200µM):

(41) BE org = SR 0.4 2.1 .

Given that BEorg depends on the [O2]bw (Lasaga and Ohmoto, 2002; Katsev and Crowe, 2015), we adopt the following formulation for sediments underlying less-oxygenated bottom waters ([O2]bw<30µM) (Dale et al., 2012):

(42) BE org = b 2 + b 1 - b 2 1 + SR / a ,

where a=0.019, b1=0.05, and b2=0.7, respectively. For intermediate [O2]bw levels, BEorg is evaluated as a function of [O2]bw with a log-linear interpolation method. Note that the original CANOPS model (Ozaki and Tajika, 2013; Ozaki et al., 2011) adopted Eq. (41) without considering the O2 dependency, whereas more recent versions employ Eq. (42) for both oxic and anoxic sediments with different values of a, b1, and b2. In CANOPS-GRB v1.0, we adopted both equations because this provided a more accurate reproduction of Corg burial distribution in the modern ocean (Sect. 3.2.2).

Figure 6Burial efficiency of organic carbon (BEorg) as a function of sedimentation rate (SR). The dots denote the observational data compiled from literature survey (Betts and Holland, 1991; Canfield, 1993; Henrichs and Reeburgh, 1987; Tromp et al., 1995; Hartnett et al., 1998). The color represents the O2 concentration of bottom water, [O2]bw, with gray dots for the unknown [O2]bw value. Blue and red lines are the relationship for well-oxygenated ([O2]bw>200µM) (Henrichs and Reeburgh, 1987) and anoxic ([O2]bw<30µM) marine sediments adopted in the CANOPS-GRB model, whereas the solid gray line is a previously proposed empirical relationship of Betts and Holland (1991).

Sedimentation rate depends strongly on water depth and distance from shore (Hedges et al., 1999), and we apply the relationship between water depth, z, and the reference value of SR shown by (Tromp et al., 1995) (Fig. 7).

(43) z = 2700 × erfc 2.1 + log SR .

Using these formulas with seafloor topography (Fig. 4a) and jorgdep (Eq. 39), we can calculate jorgb,ocn for each ocean depth. In the CANOPS-GRB model, we also introduce an erosion factor, fR, representing the global weathering and sedimentation rate (Sect. 2.5.3). Given the intimate coupling between global erosion rate and mass transfer from continents to the ocean, SR scales with the erosion factor (fR=1 for our reference run):

(44) SR ( z ) = f R SR ( z ) .

Figure 7Sedimentation rate as a function of water depth. Data (black dots) were compiled from literature survey (Colman and Holland, 2000; Baturin, 2007; Betts and Holland, 1991; Tromp et al., 1995; Cha et al., 2005; Reimers et al., 1992). Black line represents the relationship assumed in the CANOPS-GRB model. Previously estimated empirical relationships (Middelburg et al., 1997; Tromp et al., 1995) are also shown.

Organic matter that is not buried is subject to decomposition. The benthic decomposition rate at each water depth is given as follows:

(45) j recy sed = j org dep - j org b , ocn = 1 - BE org j org dep .

The respiration pathway used in the benthic decomposition is evaluated based on semi-empirical relationships obtained by 1-D early diagenesis models (see below). The fraction of aerobic degradation in total sedimentary respiration, faero, is calculated based on oxygen exposure time (τOET):

(46) f aero = 1 - f deni 1 - e - k τ OET ,

where fdeni denotes the fraction of denitrification and k is an empirical constant. τOET is given by

(47) τ OET = OPD SR ,

where OPD is the oxygen penetration depth (cm) and SR denotes a linear sedimentation rate (cm yr−1). In the CANOPS-GRB model OPD is calculated by a simplified parametric law obtained from a 1-D early diagenetic model of C and O2. We performed a series of experiments (n=5652) in order to parameterize OPD as a polynomial function with the following variables: sedimentation rate SR (cm yr−1), bottom-water O2 concentration [O2]bw (µM), depositional flux of POC jorgdep (mmol C cm−2 yr−1), and bottom water temperature Tbw (C). The variables are allowed to vary over a parameter space spanning 10−4 cm yr-1<SR<101 cm yr−1, 100µM< [O2]bw<103µM, 10−4 mmol C cm−2 yr-1<jorgdep<101 mmol C cm−2 yr−1, and 0 C <Tbw<30C.

(48) log OPD = a 0 + a 1 log SR + a 2 log [ O 2 ] bw + a 3 log j org dep + a 4 ( log SR ) 2 + a 5 ( log [ O 2 ] bw ) 2 + a 6 ( log j org dep ) 2 + a 7 ( log SR ) ( log [ O 2 ] bw ) + a 8 ( log [ O 2 ] bw ) ( log j org dep ) + a 9 ( log S R ) ( log j org dep ) + a 10 T bw ,

where a0=-2.24869, a1=0.110645, a2=1.12569, a3=-0.281005, a4=0.014827, a5=-0.124721, a6=0.0894604, a7=0.00279531, a8=-0.127797, a9=0.0017995, and a10=0.0085171. This parametric fit provides a rapid means of obtaining OPD from a 1-D early diagenetic model of C and O2 (Fig. 8). Note that Eq. (48) is verified for [O2]bw>1µM. When bottom water O2 concentration is lower than 1 µM, OPD is set at zero.

Figure 8The correlation between the simulated OPD and the OPD obtained from an empirical relationship of Eq. (48) (n=5652). The gray line denotes the 1:1 line (r2=0.9595).


Phosphorus cycling

Marine P inventory is controlled not only by the riverine P input flux from land but also by the efficiency of P recycling in marine sediments (Van Cappellen and Ingall, 1994). Because the estimated P diffusive flux from seafloor sediments is much greater than the riverine P flux (Delaney, 1998; Hensen et al., 1998; Ruttenberg, 2003; McManus et al., 1997; Wallmann, 2003, 2010), changes in diagenetic processes affecting P recycling and burial in marine sediments could have a significant impact on global oceanic biogeochemical cycles.

Figure 9Schematic of P burial in marine sediments. The primary source of P to the sediment is the deposition of organic matter, which represents the C:P ratio of primary producers, α. Phosphorus deposited in sediments is a subject of decomposition and sink switching. Most of the deposited organic P is decomposed before burial and the DIP released to pore waters diffuses to the bottom water. A fraction of the liberated P is trapped by iron hydroxides or buried as authigenic minerals (e.g., carbonate fluorapatite). Three reactive phases, organic P (Porg), Fe-bound P (PFe), and Ca-bound P (PCa), are considered in the CANOPS-GRB model. The burial of these species are redox dependent: burial efficiency is affected by bottom water O2 concentration. Because of the sink switching, sedimentary Corg/Preac, rather than Corg/Porg ratios, provides a correct measure of the retention vs. diffusive loss of remineralized P.


A schematic of benthic P cycling is shown in Fig. 9. The majority of organic matter delivered to the sediment–water interface is regenerated (Jahnke, 1996), but a fraction of DIP released via respiration to pore waters is redistributed to other phases such as iron hydroxide or carbonate fluorapatite within the sediments. This mechanism is known as “sink-switching” (e.g., Anderson et al., 2001; Filippelli, 2001), and results in P burial other than organic P playing a more important role in the total P sink (Ruttenberg, 1993, 2003; Compton et al., 2000). Three different P pools are considered in the CANOPS-GRB model: organic P (Porg), Fe-hydroxide-bound P (PFe), and authigenic Ca-bound P (PCa). The sum of these pools is defined as biologically reactive P (Preac) (bioavailable in the ocean to fuel primary productivity). The marine C and P cycles are coupled not only through the Corg/Porg ratio of POM (α) but also through the Corg/Preac ratio of marine sediments. It is important to note that, as argued by Anderson et al. (2001), the fundamental measure of the retention vs. diffusive loss of remineralized P is not the sedimentary Corg/Porg ratio but the ratio of Corg/Preac.

Field observations of marine and lacustrine sediments have revealed that the burial efficiency of P depends on the redox conditions of bottom waters (Ingall and Jahnke, 1994): phosphorus retention potential is suppressed under anoxic bottom water conditions. Elevated Corg/Preac ratios observed in permanently anoxic environments suggest preferential regeneration of P relative to C under these conditions (Algeo and Ingall, 2007; Anderson et al., 2001; Colman and Holland, 2000; Filippelli, 2001; Ingall and Jahnke, 1997). In the CANOPS-GRB model, P benthic regeneration rate is calculated at each sediment segment based on the POP depositional flux density jPdep(=jorgdep/α) and P burial efficiency, which is a function of both [O2]bw and SR. We assume the following formulation for the Corg/Porg ratio of the buried organic phase on the basis of previous studies (Slomp and Van Cappellen, 2007; Van Cappellen and Ingall, 1994, 1996):

(49) ( C org / P org ) b = ( C org / P org ) b oxic ( C org / P org ) b anox ( C org / P org ) b anox [ O 2 ] bw oxic + 1 - [ O 2 ] bw oxic ( C org / P org ) b oxic f τ for  O 2 bw < oxic ( C org / P org ) b oxic f τ for  O 2 bw oxic ,

where (Corg/Porg)boxic and (Corg/Porg)banox denote organic C/P ratios for fully oxic and anoxic conditions and oxic (=250µM) is a threshold value of [O2]bw below which preferential P regeneration occurs. (Corg/Porg)boxic is set to twice the value of the Redfield ratio, and (Corg/Porg)banox is an upper limit assumed for organic matter buried under fully anoxic overlying waters, estimated as 40 times larger than the Redfield value on the basis of previous studies on various ancient anoxic basin sediments (Slomp and Van Cappellen, 2007; Ingall et al., 1993). In Eq. (49), we also include the dependence of buried Corg/Porg ratio on SR, expressed as fτ. Modern observations suggest that SR is a one of the major factors influencing the preservation vs. remineralization of sedimentary organic C and P. Organic C preservation in marine sediments tends to be enhanced at higher SR. In contrast, the Corg/Porg ratio of sedimentary organic matter shows a nonlinear relationship with respect to SR (Ingall and Cappellen, 1990) (Sect. 3.2.3), suggesting more complex behavior of benthic P cycling. Specifically, in the pelagic deep ocean, preferential P regeneration is not observed, likely due to the long timescale of diagenesis prior to burial (Ingall and Cappellen, 1990). In the CANOPS model, fτ is formulated as follows:

(50) f τ = 0.5 + 0.5 exp - 0.001 cm SR .

Specifically, the Corg/Porg ratio approaches the Redfield value for oxygenated pelagic sediments.

The burial flux density of Porg can be calculated as the Corg burial flux density divided by (Corg/Porg)b:

(51) j P org b = j org b , ocn C org / P org b .

The burial efficiency of Porg can be written as follows:

(52) BE P org j P org b j P dep = j P org b j org dep / α = α BE org j org b / j P org b = α BE org C org / P org b ,

where jPdep denotes the POP settling flux density to sediments, which is coupled to the C/P stoichiometry of POM (=jorgdep/α).

Under oxic bottom-water conditions, remineralized organic P can be trapped efficiently at the sediment–water interface by ferric iron phases. In contrast, under anoxic bottom water conditions, a lack of ferric iron phases allows most mineralized P to diffuse out of the sediment. This redox-dependent P burial is assumed to be linearly proportional to the [O2]bw (Slomp and Van Cappellen, 2007):

(53) j P - Fe b = BE P org oxic [ O 2 ] bw oxic j P org dep  for  [ O 2 ] bw < oxic BE P org oxic j P org dep  for  [ O 2 ] bw oxic ,

where BEPorgoxic denotes the burial efficiency of Porg under well-oxygenated bottom water conditions ([O2]bw>oxic). We assume that the retention potential of PFe in sediments overlain by oxic bottom waters is comparable to that of Porg (Ruttenberg, 1993).

Some authors have also proposed that authigenic P burial, the dominant process for P burial today (Ruttenberg, 1993), depends on the redox conditions of the bottom water (Slomp and Van Cappellen, 2007; Slomp et al., 2002). In the CANOPS-GRB model, we adopt the following redox dependence used by Slomp and Van Cappellen (2007):

(54) j P - Ca b = 2 BE P org oxic a auth + 1 - a auth [ O 2 ] bw oxic j P org dep for  [ O 2 ] bw < oxic 2 BE P org oxic j P org dep for  [ O 2 ] bw oxic .

We assume that Porg, PFe and PCa account roughly for 25 %, 25 %, and 50 %, respectively, of the total reactive P buried in oxygenated sediments (Ruttenberg, 1993). Therefore, the burial efficiency of authigenic P phases is larger than that of Porg and PFe by a factor of 2. The redox dependency of authigenic P burial is controlled by a parameter, aauth. There is still great uncertainty as to the sensitivity of P retention efficiency of authigenic P phases to changing redox conditions. For instance, a recent modeling study suggests that the burial of authigenic P is influenced not only by the redox state of bottom water but also by seawater chemistry (especially Ca2+) (Zhao et al., 2020), temperature, and pH (Papadomanolaki et al., 2022). In our reference run, we set aauth at 1, i.e., no redox dependency for authigenic P burial.

When above formulations are adopted, the ratio

(55) C org / P reac = Marine C org  burial rate Marine P reac  burial rate ,

varies between 63 and 370 as a function of ocean redox state. This is in the range of an estimation derived from various observations of modern and ancient sediments (Papadomanolaki et al., 2022; Algeo and Ingall, 2007). Given that the continental shelves are a main locus of reactive P burial, the separate treatment of continental shelves and margin sediments from the pelagic ocean could affect the nonlinearity of the redox-dependent P cycle. However, this was left as one of the subjects of future work.

Nitrogen cycling

The benthic denitrification rate is estimated with a semi-empirical relationship (Middelburg et al., 1996). Middelburg and colleagues performed a series of experiments (n=2000) with a 1-D early diagenetic model of C–N–O2 to parameterize benthic denitrification jdenised (µmolCcm-2d-1) as a polynomial function using jorgdep (µmolCcm-2d-1), z (m), and bottom-water concentrations of dissolved O2 and NO3- (µM):

(56) log j deni sed = c 0 + c 1 log j org dep + c 2 log j org dep 2 + c 3 log [ NO 3 - ] bw log [ O 2 ] bw + c 4 log [ NO 3 - ] bw + c 5 log [ O 2 ] bw + c 6 log z + c 7 log j org dep log [ O 2 ] bw ,

where c0=-2.2567, c1=-0.1850, c2=-0.2210, c3=-0.3995, c4=1.2500, c5=0.4721, and c6=-0.0996, c7=0.4256. This polynomial function was obtained by examining a parameter space spanning 50 m <z<6000 m, 10 µM< [O2]bw<350µM, and 1 µM< [NO3-]bw<60µM. jorgdep was allowed to vary within 2 orders of magnitude at each water depth (Middelburg et al., 1996). As pointed out by Romaniello and Derry (2010), the predicted contribution of denitrification to total decomposition fdeni (=jdenitrsed/jrecysed) can sometimes exceed 100 % for [O2]bw<10µM. When the fraction of benthic denitrification to total decomposition exceeds 90 %, benthic denitrification is limited in order to avoid unphysical values (Ozaki and Tajika, 2013; Romaniello and Derry, 2010).

The burial flux density of Norg is calculated by molar ratio of C to N of buried sediments, (Corg/ Norg)b, and the burial flux of Corg:

(57) j N org b = j org b , ocn C org / N org b .

We assumed an average ratio of 10, which is observed in the Washington and Mexico margin (Hedges et al., 1999; Hartnett and Devol, 2003).

Sulfur cycling

The fractions of MSR and methanogenesis to total decomposition of organic matter in marine sediment are given by


The production rate of hydrogen sulfide in sediment, jH2Ssed (mol S m−2 yr−1), is given by

(60) j H 2 S sed = 1 2 f MSR j recy sed + j AOM ,

where jAOM denotes the production rate of sulfide via AOM:

(61) j AOM = 1 2 f meth [ SO 4 2 - ] bw [ SO 4 2 - ] bw + K MSR j recy sed .

Here we assume that AOM is proportional to the CH4 production rate with a sulfate-dependent term.

The rate of pyrite precipitation in sediments would be proportional to the sulfide production rate at the sediment-water interface:

(62) j pyr b , sed = e pyr j H 2 S sed ,

where the proportional coefficient, epyr, is the pyrite burial efficiency. The rate of MSR is a function of the marine redox state, [SO42-], and the availability of degradable organic matter. In well-oxygenated modern oceans most sulfide produced in sediments is re-oxidized and only a few percent of total sulfide is buried as pyrite (Canfield, 1991; Lin and Morse, 1991; Turchyn and Schrag, 2004; Bowles et al., 2014; Jørgensen, 1982). It has been pointed out that efficient oxidation of sulfide is promoted by animal bioturbation (Berner and Westrich, 1985; Canfield and Farquhar, 2009). In contrast, the value of epyr for anoxic sediments is much greater due to the absence of bioturbation and enhanced sulfide production. We assume that epyr asymptotes toward unity with decreasing the bottom water [O2] (Tarhan et al., 2015):

(63) e pyr = e pyr max - e pyr max - e pyr tanh [ O 2 ] bw ,

where epyrmax (=1 in our reference run) denotes the maximum pyrite precipitation efficiency in anoxic sediments. The reference value, epyr, was calibrated using a present-day control simulation such that the present-day seawater [SO42-] is ∼29 mM. The obtained value of 0.117 is generally consistent with modern observations (Bottrell and Newton, 2006; Tarhan et al., 2015; Turchyn and Schrag, 2006) (see Sect. 3.2.5). Although our approach does not provide a mechanistic description of the complex process of pyrite precipitation, it is suitable for many purposes.

Early diagenetic modeling for quantifying the OPD

A simple 1-D early diagenetic model of C and O2 is employed to obtain the parameterization of OPD (Eq. 48). The 100 cm thick sediment is vertically divided into 50 layers with an uneven grid. The grid size increases from the sediment–water interface (Δz=0.25 mm) to the maximum simulated sediment depth (Δz=1.6 cm). The diagenetic model calculates transport and biogeochemical transformation processes at each grid point within these sediment columns as well as the sedimentary burial and recycling fluxes at the model boundaries. The one-dimensional mass conservation equation for POC (wt. %) and dissolved O2 is given by


where DO2 is the diffusion coefficient of O2, SR is the sedimentation rate, and φ is porosity, which is assumed to be constant over the entire sediment column for simplicity. Bioturbation is formulated as a diffusive process with a coefficient Dbio. The effective diffusion coefficient of O2 is then given by

(66) D O 2 = D O 2 T = 0 × 1 + ν O 2 T bw θ 2 + D bio ,

where DO2T=0 denotes a tracer diffusion coefficient in seawater of 0 C and vO2 is a coefficient for temperature dependence of molecular diffusion coefficient. The in situ diffusion coefficient is further corrected for tortuosity θ, which is related to pore water resistivity and porosity via the following expressions (Colman and Holland, 2000; Tromp et al., 1995; Berner, 1980):


where F is the formation factor – defined as the ratio of bulk sediment resistivity to interstitial water resistivity – and m is an empirical constant, varying with sediment type. We assumed the average value for unconsolidated muds (m=2.7) in this work (Tromp et al., 1995). The particle mixing coefficient for bioturbation Dbio is formulated as a function of both sediment accumulation rate and bottom-water O2 concentration (Tromp et al., 1995; Wallmann, 2003):

(69) D bio = 10 1.63 + 0.85 log SR [ O 2 ] bw [ O 2 ] bw + K O 2 .

At the bottom of the sediment column, a no-flux condition was applied. The parameters used in the 1-D early diagenetic model are tabulated in Table 6.

Table 6Parameters used in the 1-D early diagenetic model.

Download Print Version | Download XLSX

2.4.5 Air–sea exchange

To calculate the gas exchange of O2, H2S, NH3, and CH4 across the air–sea interface, we employed a stagnant film model (Liss and Slater, 1974). The flux of a gas X across the air–sea interface is controlled by the difference in partial pressure between the atmosphere and surface waters, which can be described by the following formula:

(70) J X air - sea = v X pis [ X ] aq - [ X ] sat ,

where vXpis, [X]aq, and [X]sat denote piston velocity, the dissolved concentration of species X, and the saturation concentration of species X, respectively. For O2, the saturation concentration is calculated based on solubility (Garcia and Gordon, 1992; Sarmiento and Gruber, 2006) and partial pressure:

(71) [ O 2 ] sat = 1000 22.3916 e l p O 2 p O 2 ,



with T in C. The constants are A0=2.00907, A1=3.22014, A2=4.0501, A3=4.94457, A4=-0.256847, A5=3.88767, B0=-6.24523×10-3, B1=-7.3761×10-3, B2=-1.0341×10-2, B3=-8.17083×10-3, and C0=-4.88682×10-7. The erroneous A3×Ts2 term in the original equation (Garcia and Gordon, 1992) was left out (Sarmiento and Gruber, 2006).

For CH4, H2S, and NH3, [X]sat is given by (Kharecha et al., 2005)

(74) [ X ] sat = K X Henry p X ,

where KXHenry, and pX denote Henry's law coefficient and the partial pressure of species X, respectively. The temperature dependence of X's solubility is expressed as

(75) K X Henry = K X Henry exp K X T 1 T - 1 298.15 ,

where KXHenry denotes the Henry's law coefficient of species X at 25 C and KXT is the temperature-dependence constant.

[X]aq is the dissolved concentration of X. [H2S]aq and [NH3]aq can be written as follows:


where [ΣH2S] = [H2S] + [HS] and [ΣNH3] = [NH4+] + [NH3]. KH2Sdis and KNH3dis are the dissociation constant, defined as follows:


Given values of KH2Sdis, KNH3dis and pH (Millero et al., 1988; Yao and Millero, 1995), [H2S]aq and [NH3]aq can be calculated.

Atmospheric concentrations of H2S and NH3 are set at 0. H2S and NH3 flow past the surface layer of the ocean to the atmosphere are converted to an equal influx of SO42- and NO3- to the surface ocean. The parameters used in the stagnant film model are tabulated in Table 7.

If atmospheric O2 levels are lower than ∼1 %, PAL spatial heterogeneity of the gas exchange flux is expected (Olson et al., 2016); for example, primary productivity (and O2 generation) would be more active in coastal regions than open-ocean gyres. Because our ocean model resolves only two regions for the surface oceans (low- to mid-latitude region L and high-latitude region H), it tends to overestimate the oxidation of reductants in surface mixing layers. To mitigate this model limitation for the CH4 degassing flux, the aerobic oxidation rate of CH4 is decreased to 1×10-7 of the standard value in surface layers (Ozaki et al., 2019a).

Table 7Parameters used in the air–sea exchange module of CANOPS-GRB.

Download Print Version | Download XLSX

2.5 Land model

2.5.1 Net primary productivity

Terrestrial NPP is scaled by global land biomass V normalized to the modern value:

(80) J NPP lnd = V × J NPP lnd , ,

where the present value of terrestrial NPP is set at 60 Gt C yr−1 (Prentice et al., 2001). The global land biomass is a function of atmospheric O2 levels:

(81) V = f UV f fire f O 2 ,

where fO2 represents the direct effect of atmospheric O2 concentration on the C3 plant growth and ffire denotes the effect of fires on land biota (Bergman et al., 2004; Lenton and Watson, 2000b):


Here kfire (=3; Lenton, 2013) is the fire frequency constant and “ignit” is an ignition factor representing the fire frequency as a function of oxygen (Lenton, 2013; Lenton et al., 2018; Lenton and Watson, 2000b):

(84) ignit = min max c 1 p O 2 - c 2 , 0 , c 3 ,

with c1=48, c2=9.08, and c3=5 (Lenton, 2013). CANOPS-GRB also includes an additional factor fUV representing the effect of UV on the terrestrial biosphere as a function of atmospheric O2 levels (Ozaki and Reinhard, 2021):

(85) f UV = tanh p O 2 PAL c UV ,

where cUV is a model parameter, which in our standard model is set at 1 % PAL, meaning that terrestrial plant activity is suppressed when atmospheric O2 is lower than a few % PAL.

2.5.2 Terrestrial biogeochemical cycles

Phosphorus weathering flux, JPw (Eq. 2), is treated as a boundary condition. A fraction of weathered P is ultimately buried as terrigenous organic matter (Eq. 3), whereas the remaining fraction is delivered to the ocean via rivers (Eq. 4). In the CANOPS-GRB model, the reference value of JPr (=0.155 Tmol P yr−1) is tuned so that the oceanic P level of the reference state is consistent with modern observations. The burial rate of terrigenous organic matter (in terms of C) can be written as follows:

(86) J org b , lnd = C org / P org lnd J P b , lnd ,

where (Corg/ Porg)lnd (=1000) is the average C/P burial ratio of terrigenous organic matter (Bergman et al., 2004). In this study, the reference value of Jorgb,lnd was set at 3 Tmol C yr−1. By combining Eqs. (3), (4), and (86) for the reference state, the proportional coefficient k11 of Eq. (3) is determined by the reference state as follows:

(87) k 11 = J org b , lnd , J org b , lnd , + C org / P org lnd J P r , = 0.0189 .

The value of k11 is treated as a constant in this study.

Almost all organic matter produced by terrestrial NPP is decomposed before burial. The total decomposition rate is given by

(88) J org r , lnd = J NPP lnd - J org b , lnd .

CANOPS-GRB includes aerobic respiration and methanogenesis as respiration pathways for terrigenous matter, and the CH4 flux from the terrestrial ecosystem to the atmosphere is evaluated with the assumption that it is proportional to the burial rate of terrigenous organic matter:

(89) J CH 4 lnd = J org b , lnd J org b , lnd , J CH 4 lnd , ,

where the reference value was set at 1 Tmol CH4 yr−1. The net flux of CO2, O2 and CH4 from the terrestrial ecosystem to the atmosphere can be written as follows:


where gO2 and gCH4 denote the fraction of organic matter decomposed by aerobic respiration and methanogenesis, respectively. δ represents the fraction of methane that is consumed by aerobic methanotrophy that is a function of O2:

(93) δ = M O 2 atm M O 2 atm + K O 2 ,

where KO2=0.273×1018 mol (Goldblatt et al., 2006). A fraction of organic matter decomposed by methanogenesis, gCH4, can be calculated based on Eqs. (89) and (92). Following this, gO2 is determined from 1-gCH4.

2.5.3 Weathering

The oxidative weathering of continental crust is a major oxygen sink on geologic timescales, providing a fundamental control on atmospheric O2 levels. The weathering rate in the model is assumed to be proportional to sedimentary reservoir size and a global erosion factor, fR, expressing the effect of continental denudation and erosion on terrestrial weathering:


where Jorgw and Jpyrw denote the oxidative weathering of organic carbon and pyrite, respectively, and forgwO2 and fpyrwO2 represent the O2 dependency. For the oxidative weathering of organic matter, previous biogeochemical models have adapted a (pO2PAL)0.5 relationship (Bergman et al., 2004; Lasaga and Ohmoto, 2002). In this study, we employ alternative empirical relationships based on results obtained from a 1-D weathering model (Bolton et al., 2006; Daines et al., 2017):


where Korgw and Kpyrw denote half-saturation constants (Korgw=0.334 and Kpyrw=0.017) and corgw and cpyrw are normalized constants (corgw=1.334 and cpyrw=1.017), respectively. The Monod-type relationship captures the fact that the rate of oxidative weathering reaches its maximum as determined by the erosion rate under highly oxygenated conditions (i.e., transport-limited regime). For example, due to the fast dissolution kinetics of pyrite, oxidative weathering can be regarded as transport-limited under modern conditions (Bolton et al., 2006) (Fig. 10). In the CANOPS-GRB model, Jorgw is calibrated based on the global redox budget of the reference run (see Sect. 2.3.5).

It is important to note that above equations ignore the possible importance of microbial activity and temperature on the rate of oxidative weathering (Petsch et al., 2001; Soulet et al., 2021). Both represent important topics for future research.

Figure 10O2 dependency of the oxidative weathering rate of organic matter and pyrite sulfur. The gray line denotes the (pO2PAL)0.5 relationship assumed in previous biogeochemical models (Lasaga and Ohmoto, 2002; Daines et al., 2017). Solid and dashed black lines represent the empirical Monod-type relationships for oxidative weathering of organic matter (solid) and pyrite sulfur (dashed) based on the results obtained from a 1-D weathering model (Bolton et al., 2006; Daines et al., 2017), which are adopted in the standard model of the CANOPS-GRB model. PAL stands for present atmospheric level.

The present riverine flux of sulfur, JSr, is estimated at 2.6 Tmol S yr−1 (Raiswell and Canfield, 2012), representing the dominant source to the oceans. Riverine flux is written as the sum of the gypsum weathering flux Jgypw and the oxidative weathering of pyrite Jpyrw and depends directly or indirectly on the oxidation state of the atmosphere:

(98) J S r = J gyp w + J pyr w .

Based on previous studies (Berner, 2009; Wortmann and Paytan, 2012; Bergman et al., 2004; Markovic et al., 2015), a 3 : 1 ratio in modern rivers of SO42- from gypsum vs. pyrite weathering is assumed. Gypsum weathering flux is assumed to be proportional to its sedimentary reservoir size, GYP, and fR:

(99) J gyp w = f R GYP GYP J gyp w ,

where * represents the present value.

In the previous version of the CANOPS (Ozaki et al., 2019a), oxidative weathering of pyrite was divided into biogenic and abiotic weathering fluxes. In this study, we simplify this (Eq. 95). In addition, oxidative weathering of Fe(II)-bearing minerals is ignored in this study, which simplifies the framework of the global O2 budget (Sect. 2.3.5).

2.5.4 Volcanic degassing

Volcanic outgassing fluxes of carbon and sulfur are assumed to be proportional to their respective crustal reservoir sizes:


We set the reference value of the volcanic outgassing flux of organic carbon, Jorgm,, at 1.25 Tmol C yr−1 (Bergman et al., 2004). The estimates of modern volcanic fluxes of sulfur fall within the range of ∼0.3–3 Tmol S yr−1 (Kagoshima et al., 2015; Catling and Kasting, 2017; Raiswell and Canfield, 2012; Walker and Brimblecombe, 1985). We adopted a recent estimate of 0.8 Tmol S yr−1 (Kagoshima et al., 2015).

2.5.5 Sedimentary reservoirs

We extend the original model framework to the explicit calculation of the secular evolution of the sedimentary reservoirs, linking the biogeochemical cycles in the ocean–atmosphere system to the rock cycle. The mass balance equation for sedimentary organic carbon (ORG) can be written as follows:

(103) d ORG d t = J org b - J org w - J org m ,

where Jorgb denotes the sum of the burial rate of marine and terrigenous organic matter (Jorgb,ocn+Jorgb,lnd), the primary source of sedimentary organic carbon. Primary outputs are oxidative weathering, volcanic outgassing, and metamorphism. Previous estimates of the present reservoir size of ORG fall in the range of 1000–1300 Emol (1 E = 1018) (Berner, 1989; Garrels and Perry, 1974; Mackenzie et al., 1993). We assumed 1250 Emol for the reference value of ORG.

The sedimentary reservoir sizes of pyrite sulfur (PYR) and gypsum sulfur (GYP) are also written as the balance between the input (burial) and outputs (weathering and outgassing):


where Jpyrb represents the sum of pyrite precipitation rates in the water column and sediments, Jpyrb,wc+Jpyrb,sed. Previous estimates of present reservoir sizes of GYP and PYR fall in the range of 77–300 Emol and 155–300 Emol (Berner, 2006; Bottrell and Newton, 2006; Yaroshevsky, 2006; Kump, 1989; Lasaga, 1989; Holser et al., 1989; Sleep, 2005; Schlesinger and Bernhardt, 2013), respectively. We adopted 200 and 200 Emol for GYP and PYR.

2.6 Atmosphere model

2.6.1 Hydrogen escape

The rate of hydrogen escape is assumed to be diffusion limited as it is today. Thus, the total concentration of all H-bearing compounds in the lower stratosphere determines the rate of hydrogen escape (Walker, 1977). For Proterozoic-Phanerozoic atmospheres, CH4 appears to have been the dominant hydrogen-bearing species in the stratosphere, and the flux, JHesc (mol yr−1), is calculated as

(106) J Hesc = s M CH 4 atm ,

where MCH4atm denotes the abundance of CH4 in the atmosphere (mol) and s (=3.7×10-5 yr−1) is a proportional coefficient (Goldblatt et al., 2006).

2.6.2 Photochemistry

CANOPS-GRB includes parameterized O2–O3–CH4 photochemistry that allows quantification of the abundances of atmospheric O2 and CH4. The rate of oxidation of CH4 is calculated by the following empirical parameterization that was obtained from a 1-D photochemistry model (Claire et al., 2006):

(107) J CH 4 ox = k CH 4 ox M O 2 atm M CH 4 atm ,

where MO2atm and MCH4atm denote the abundance of O2 and CH4 in the atmosphere (mol). The reaction rate kCH4ox (mol−1 yr−1) is expressed as a polynomial function of the reservoir sizes of O2 and CH4 (Ozaki and Reinhard, 2021):

(108) log k CH 4 ox = α 0 j + α 1 j φ O 2 + α 2 j φ O 2 2 + α 3 j φ O 2 3 + α 4 j φ O 2 4 + α 5 j φ O 2 5 + α 6 j φ O 2 6 ,

where aij represents fitting coefficients for given atmospheric CH4 levels and ϕO2 is logpO2 (in bar) (supplementary Table 4 of Ozaki and Reinhard, 2021). The oxidation rate was evaluated using Fig. 3 of Claire et al. (2006), showing the oxidation rate as a function of pO2 and pCH4. We took the relationship between kCH4ox and pO2 for pCH4 of 10−6, 10−5, 10−4, 10−3, and 2×10-3 bar, and kCH4ox is calculated as a function of pO2 and pCH4 with a log-linear interpolation method. Note that the default photochemical parameterization presented above limits the applicability of CANOPS-GRB v1.0 to Earth-like planets around the Sun-like host stars. Current work is focused on elaborating parameterized photochemistry across a wider range of spectral energy distributions.

2.6.3 Mass balance

CANOPS-GRB accounts for the atmospheric concentrations of O2 and CH4. The atmospheric concentration of O2 is determined by the biogenic source (from the ocean and terrestrial ecosystems) and the consumption through the series of oxidation reaction (the continental weathering of kerogen and pyrite, volcanic outgassing, and photochemical oxidation of methane):

(109) d M O 2 atm d t = J O 2 air - sea + J O 2 air - lnd - J Hesc + 2 J CH 4 ox - J org w + J org m - 2 J pyr w + J pyr m ,

where MO2atm denotes the mass of O2 in the atmosphere (moles) and the first and second term on the right-hand side represents the biogenic flux of O2 from marine and terrestrial ecosystems. The third term denotes O2 consumption via photochemistry, and the fourth and fifth terms are the O2 consumption via organic C and pyrite S sub-cycles. Note that the hydrogen escape to space is represented as the O2 sink because the hydrogen escape via CH4 followed by the oxidation of carbon to CO2 is represented as

(110) CH 4 + O 2 + h ν 4 H + CO 2 .

On the other hand, the photochemical oxidation of CH4 can be written as follows:

(111) CH 4 + 2 O 2 + h ν 2 H 2 O + CO 2 .

Thus, the hydrogen escape to space represents the net gain of oxidizing power to the system (see Eq. 12).

The abundance of CH4 in the atmosphere, MCH4atm, is determined by input from the ecosystems and the consumption of CH4 via photolysis, as well as by the hydrogen escape:

(112) d M CH 4 atm d t = J CH 4 air - sea + J CH 4 air - lnd - J Hesc + J CH 4 ox .

No abiotic CH4 input via hydrothermal systems is included.

3 Validation against the modern global ocean

Here, a steady-state simulation mimicking the present-day condition was run to evaluate the overall performance of CANOPS-GRB. To do this, the ocean model was run until reaching steady state, assuming the present atmospheric O2 level and reference values of boundary fluxes (weathering and volcanic fluxes). The simulated circulation and biogeochemistry for the modern global ocean was compared with modern oceanographic observations from the Global Ocean Data Analysis Project (Key et al., 2015; Olsen et al., 2016).

3.1 Distribution of circulation tracers

Comparisons of model output with circulation tracers, such as potential temperature (θ) and radiocarbon (Δ14C), permit a test of the physical exchange scheme. Figure 11 depicts the simulated patterns of physical tracers with observational data. The physical circulation in the model generally agrees well with oceanic observations, although we note that model temperatures for low- to mid-latitudes above 1000 m water depth tend to be higher than observed because temperature distribution in the real ocean is strongly controlled by vertical structure and advective processes that are not captured in our simple circulation scheme. Despite this model limitation, the modeled temperature distribution generally reproduces the observed distribution. The Δ14C minimum in the model for the low- to mid-latitude region corresponds well with observations. The modeled background radiocarbon for young deep waters (about -150±25 ‰) is closer to the value for the Southern Ocean (approximately −150 ‰) than for North Atlantic deep waters (approximately −80 ‰), and old deep waters (-200±15 ‰) correspond to the South Pacific. We conclude that the simulated circulation tracers generally match well with ocean data.

Figure 11Simulated steady-state depth profiles of (a) potential temperature, θ, (b) radiocarbon, Δ14C, (c) DIP (dissolved inorganic phosphorus, PO43-), and (d) dissolved oxygen, O2. Solid and dashed white lines denote the simulated profiles for the LD and HD regions, respectively. Simulation results are compared with the dataset from the Global Ocean Data Analysis Project (GLODAP) database (GLODAPv2_2019; Key et al., 2015; Olsen et al., 2016). The color represents the density of observational points.

3.2 Ocean biogeochemistry

Having demonstrated that CANOPS-GRB's ocean circulation model does a reasonable job of representing water mass exchange, we next assess the performance of the oceanic biogeochemistry model by comparing its output to ocean biogeochemical data. Model-generated global fluxes and inventories of C, P, N and S cycles are summarized in Fig. 12. These compare well with independent observational estimates. Below, we provide a brief discussion of globally integrated biogeochemical flux estimates.

Figure 12Schematics of the simulated material flow in the ocean for the reference run. (a) Organic carbon (in Gt C yr−1), (b) phosphorus (in Tmol P yr−1), (c) nitrogen (in Tg N yr−1), and (d) sulfur (in Tmol S yr−1). NPPocn is oceanic net primary production. EX represents export production. MX is the mass of X in the ocean. τX is the residence time of X in the ocean. Pmol stands for 1015 mol, and Emol stands for 1018 mol.


3.2.1 Distribution of nutrients and oxygen

The simulated vertical profile of phosphate captures the characteristic features and values of observational data (Fig. 11c). More specifically, the distribution in the low- to mid-latitude region is more similar to that in the Pacific and Indian oceans, and distribution in the high- to mid-latitude region is similar to that in the Southern Ocean. This is a consequence of limiting high-latitude productivity (preformed DIP is 1.1 µM) which results in higher concentrations in the ocean interior. The model dissolved O2 profile for low- to mid-latitude areas shows a minimum of approximately 100 µM at water depth of 1000 m, corresponding to the oxygen minimum zone (Fig. 11d). In contrast, dissolved O2 for the high- to mid-latitude sector (HD) shows a monotonically decreasing trend. This is because of oxygen consumption via POM decomposition during downwelling.

3.2.2 Carbon cycling

The marine export and new production in our model is 9.1 Gt C yr−1 (8.36 Gt C yr−1 at L and 0.74 Gt C yr−1 at H). This is consistent with previously estimated global values of 8.5–12 Gt C yr−1 (Dunne et al., 2007; Laws et al., 2000; Sarmiento and Gruber, 2006; Heinze et al., 2009). In particular, our estimate is close to the mid-point of the previously estimated range of 9.6±3.6 Gt C yr−1 (Dunne et al., 2007). This is a marked improvement from earlier studies with box models that have underestimated marine new production by a factor of 2 or more (Archer et al., 2000; Shaffer et al., 2008). Simulated global oceanic NPP is 45.5 Gt C yr−1. This is also consistent with the previous estimated range of 44–65 Gt C yr−1 (Prentice et al., 2001; Woodward, 2007; Carr et al., 2006; Berelson et al., 2007).

Figure 13Simulated steady-state depth profiles of organic C and reactive P flux density for the LD region. In (a), gray dots denote observations of depositional/settling flux density, whereas black dots represent observations of burial flux density compiled from literature survey (Baturin, 2007; Betts and Holland, 1991; Colman and Holland, 2000; Lutz et al., 2002). Solid gray and black lines denote the simulated POC depositional and burial flux densities obtained from the reference run. (b) Gray dots denote the benthic P efflux density obtained from literature survey (Hartnett and Devol, 2003; Hensen et al., 1998; Ingall and Jahnke, 1994, 1997; McManus et al., 1997; Colman and Holland, 2000; Schenau and De Lange, 2001; Zabel et al., 1998), whereas solid gray and black lines represent the simulated benthic P efflux density and burial flux density of reactive P obtained from the reference run. The burial flux density of reactive P is not shown due to the sparseness of such observations.

The global marine POC flux depends largely on water depth. Model-generated fluxes compare well with independent estimates of deposition, burial, and regeneration. The gray line in Fig. 13a shows the simulated sinking flux density of POC in the water column for LD region, compared with observations (Archer et al., 2002; Betts and Holland, 1991; Lutz et al., 2002; Baturin, 2007). The preferential consumption of labile compounds (G1 and G2) during the settling process leads to a continuous decrease in reactivity and therefore remineralization rates from the surface ocean down to the deep. Our estimate lies well within the range of observations. The model tends to give lower fluxes than observed above 2000 m water depth and higher below 5000 m water depth. This is probably because of the assumption of homogeneous productivity in the surface ocean. In the real ocean, oceanic productivity is generally greater at the continental margins than in the pelagic gyre regions (Lutz et al., 2002). This is a model limitation, but the simulated biological pump is sufficient to describe the general characteristics of global ocean biogeochemistry.

Of total exported POC, 91 % (8.25 Gt C yr−1) is decomposed in the water column and the rest (0.85 Gt C yr−1) sinks to the sediment surface (Fig. 12a). The simulated global POC depositional flux is comparable not only with observational estimates of 0.93 Gt C yr−1 (Muller-Karger et al., 2005) and 0.67±0.48 Gt C yr−1 for off-shore regions (Dunne et al., 2007) but also with an estimate using EMIC (0.87 Gt C yr−1) (Ridgwell and Hargreaves, 2007). The depositional fluxes of Corg in marginal (<2000 m) and deep-sea sediments (>2000 m) are estimated at 0.58 and 0.27 Gt C yr−1, respectively. These estimates are slightly lower than previous estimates of 0.62–1.98 and 0.31–0.62 Gt C yr−1 (Bohlen et al., 2012; Dunne et al., 2007; Muller-Karger et al., 2005; Burdige, 2007).

In our standard run, benthic remineralization removes 7.9 % of the exported POC (0.72 Gt C yr−1), equivalent of 84 % of the global POC sedimentation rate. As a result, only 1.5 % (0.135 Gt C yr−1 or 11.3 Tmol C yr−1) of the global POC export production is ultimately buried in marine sediments. Our model demonstrates that much (91 %) of the total burial occurs on the continental margins (<2000 m water depth), where the settling flux and burial efficiency are relatively high. Previous studies (Dunne et al., 2007; Muller-Karger et al., 2005) estimated a Corg burial rate of 0.29±0.15 and >0.06±0.06 Gt C yr−1 at the margin. Our estimate of 0.123 Gt C yr−1 lies between these values, whereas our estimate for the deep sea, 0.012 Gt C yr−1, is on the lower end of previous estimates of 0.012±0.02 Gt C yr−1 (Dunne et al., 2007) and 0.017±0.005 Gt C yr−1 (Hayes et al., 2021). In addition, Sarmiento and Gruber (2006) and Hayes et al. (2021) estimated the burial rate below 1000 m as 0.02±0.006 Gt C yr−1; our estimate of 0.019 Gt C yr−1 is consistent with this. Combined with the prescribed burial rate of terrigenous Corg 0.036 Gt C yr−1 (3 Tmol C yr−1), the total burial rate is calculated to be 0.17 Gt C yr−1 (14.3 Tmol C yr−1). This is somewhat higher than previous estimates (Berner, 1982; Burdige, 2005; Muller-Karger et al., 2005), but given the rather large uncertainty we consider it defensible.

Figure 14 shows OPD as a function of water depth. Although the benthic data could be biased towards highly specific environments, such as sediments underlying upwelling areas and continental margins, our estimates capture the general features of modern observations.

Figure 14Oxygen penetration depth (OPD) as a function of water depth. Colored dots denote the observational data obtained from literature survey (Bradley et al., 2020; Donis et al., 2016; Nierop et al., 2017; Rowe et al., 2008; McManus et al., 2005; Martin and Sayles, 2014; Pfeifer et al., 2002; Hyacinthe et al., 2001; Hartnett et al., 1998; Hedges et al., 1999; Morford and Emerson, 1999; Devol and Christensen, 1993; Gundersen and Jorgensen, 1990; Sachs et al., 2009). The color represents the O2 concentration of bottom water, [O2]bw, with gray dots for the unknown [O2]bw value. The simulated OPD obtained from the reference run is shown as a black line.

3.2.3 Phosphorus cycling

The removal of phosphate from surface waters occurs through photosynthetic fixation by primary producers and subsequent export in the form of POP into deeper waters, where it is largely remineralized back into DIP. Through this process there is a vertical partitioning of DIP within the ocean with reduced surface concentrations. Phosphorus export production is 7.16 Tmol P yr−1, which is coupled with carbon according to the POM compositional ratio (C:P=106:1 for our standard model). The remineralization in the water column (6.49 Tmol P yr−1) and total sedimentation rate (0.672 Tmol P yr−1) are also proportional to those of POC. In contrast, the benthic DIP flux is decoupled from the carbon flux. Figure 13b shows modeled DIP benthic efflux and burial flux together with observed fluxes. Some observational data showing a relatively large abyssal (4–6 km) benthic flux are from upwelling regions in the South Atlantic (Hensen et al., 1998). The deviation is therefore not critical for our globally averaged model. Our model gives the total benthic efflux of DIP as 0.517 Tmol P yr−1, which is roughly three times the riverine reactive P input rate. This is within the range of previous estimates of 0.05–1.25 Tmol P yr−1 (Wallmann, 2003, 2010; Compton et al., 2000; Colman and Holland, 2000).

The preservation efficiency (here defined as burial flux divided by the export flux) of P is 2.1 %. This is higher than that of organic carbon (1.5 %), indicating that more P is trapped in marine sediments than might be expected from Redfield stoichiometry. In marine sediments overlain by oxic bottom waters, a fraction of the DIP released to pore waters from POM decomposition can be absorbed by iron oxyhydroxide or precipitated as authigenic fluorapatite (Fig. 9). Therefore, the global averaged Corg/Preac ratio of buried sediments is generally less than the Redfield of 106 (approximately 65±25 based on observations; Algeo and Ingall, 2007). The modeled global average Corg/Preac ratio of buried sediment is 73, which is consistent with this. The P burial fluxes of organic P, Fe-bound P, and authigenic P are estimated at 0.044, 0.032, and 0.079 Tmol P yr−1, respectively.

The Corg/Porg ratio of burying organic matter shows a nonlinear relationship with respect to sedimentation rate. The observed Corg/Porg ratios are generally greater than the Redfield value of 106, especially for sediments in oxygen minimum zones (OMZs), which are characterized by a high depositional flux of organic matter (Corg/Porg ratios up to 600 for the present open ocean) (Ingall and Cappellen, 1990). For example, the averaged Corg/Porg molar ratio at the Peru–Chile OMZ and Black Sea are 600, and the estimated burial ratio of sapropel S1 of Mediterranean Sea is in the range of 400–800 (Slomp et al., 2002). This reflects the preferential regeneration of P relative to C during microbial remineralization of marine organic matter and reflects the more labile nature of P biochemicals relative to most non-phosphorus-containing organic carbon compounds. Additional rationale for this observation is that P is preferentially targeted for remineralization to support subsequent biological productivity as an essential and potentially limiting nutrient. Our model demonstrates that we can reproduce the first-order relationship between Corg/Porg and sediment accumulation rate (Fig. 15).

The modeled marine DIP inventory is 2.75×1015 mol, consistent with the observational estimate of around 3×1015 mol (e.g., Delaney, 1998; Guidry et al., 2000). Given the riverine reactive P input flux of 0.155 Tmol P yr−1, the phosphorus residence time is estimated at 18 kyr, which is also consistent with previous estimates of 20 kyr or shorter (Benitez-Nelson, 2000; Ruttenberg, 2003).

Figure 15Corg/Porg ratios of buried sediments as a function of sedimentation rate. Black dots represent the observational data (Ingall and Cappellen, 1990). The simulated Corg/Porg ratios for the LD region obtained from our reference run is shown as a black line.

3.2.4 Nitrogen cycling

Nitrogen export production is 1603 Tg N yr−1, which is coupled with carbon according to the C:N stoichiometry of organic matter. Simulated N fixation required for the N balance in the ocean is 180 Tg N yr−1, which is higher than the range of previous estimates of 110–150 Tg N yr−1 (Luo et al., 2012; Gruber and Sarmiento, 1997; Galloway et al., 2004; Karl et al., 2002; Fowler et al., 2013; Duce et al., 2008; Deutsch et al., 2007), while a recent study (Großkopf et al., 2012) suggests a higher value of ∼180 Tg N yr−1. This discrepancy is partly because atmospheric deposition is ignored in the CANOPS-GRB model, which contributes 25.8 Tg N yr−1 (Wang et al., 2019). Gruber and Sarmiento (2002) estimated the preindustrial value of the total source of N as 188±44 Tg N yr−1. Our estimate of 196 Tg N yr−1 is within this range.

Nitrogen fluxes in an oxic water column are tightly coupled with the Corg fluxes, whereas decoupling appears in suboxic environments. Simulated denitrification in the water column is 102 Tg N yr−1, within the range of the observational estimates (50–150 Tg N yr−1) (DeVries et al., 2012, 2013; Brandes and Devol, 2002; Gruber, 2008; Gruber and Sarmiento, 2002; Oschlies et al., 2008; Wang et al., 2019). Modeled benthic denitrification is 62 Tg N yr−1, which is lower than the estimated range of 90–300 Tg N yr−1 (DeVries et al., 2012, 2013; Brandes and Devol, 2002; Eugster and Gruber, 2012; Devol, 2015; Wang et al., 2019) by a factor of 1.5–5, suggesting that further efforts are required to improve representation of this process. One possible explanation for this discrepancy is that our model is not sufficient to express benthic N cycling because we ignore localized upwelling regions (such as the eastern Tropical Pacific and the Arabian Sea) and coastal regions where benthic denitrification is significant POM decomposition pathway in favor of globally averaged parameterizations. The separate treatment of continental shelves and margin sediments from the pelagic ocean could improve this issue. We also ignore another denitrification mechanism: anaerobic ammonium oxidation (anammox), which will often play an important role in the loss of fixed nitrogen in marine sediments and pelagic anoxic zones (Karthäuser et al., 2021; Kuypers et al., 2005).

The modeled DIN inventory is 4.5×105 Tg N. Given the total source flux of 196 Tg N yr−1, the residence time of DIN is estimated at 2.3 kyr.

3.2.5 Sulfur cycling

MSR is a major early diagenetic pathway of carbon oxidation in organic-rich sediments deposited below oxygenated waters. For the standard run, aerobic oxidation is a dominant process in the water column, but MSR contributes 37 % of benthic degradation. CANOPS-GRB estimates a global rate of benthic sulfate reduction at 11 Tmol S yr−1. This is lower than the previously reported value of gross MSR (40–75 Tmol S yr−1; Canfield and Farquhar, 2009; Jørgensen and Kasten, 2006) but agrees better with net MSR (Bowles et al., 2014). Bowles et al. (2014) have estimated global net MSR at 6.2 and 11.3 Tmol S yr−1 for z>200 m depth and z>0 m depth, respectively. Our estimate is within this range. MSR is most pronounced on the shelf where high fluxes of organic matter to the seafloor lead to shallow OPD, high sulfide production, and consequently high pyrite precipitation (Fig. 16).

Figure 16MSR and pyrite burial flux density as a function of sedimentation rate. Gray and black dots depict observational data compilation of depth-integrated MSR flux density and pyrite burial flux density for normal (oxic) marine sediments (Berner and Canfield, 1989; Canfield, 1989; Raiswell and Canfield, 2012). The unit of sedimentation rate was converted from g cm−2 yr−1 to cm yr−1, assuming a dry bulk density of 2.5 g cm−3 and porosity of 0.9. Solid lines are the results obtained from the reference run of the CANOPS-GRB model.

In the CANOPS-GRB model, pyrite burial efficiency epyr (Sect. 2.4.4) for sediments underlying oxic bottom waters is set such that simulated seawater [SO42-] of the reference run is consistent with the modern value of 28.9 mM. The tuned value of 11.7 % agrees well with observations suggesting that pyrite precipitation rate is about 10 %–20 % of the rate of MSR (Fig. 16). Our reference value is also consistent with other estimates of 11 %–20 % (Bottrell and Newton, 2006; Tarhan et al., 2015; Turchyn and Schrag, 2006).

The sulfate inventory of our reference state is 39×1018 mol. Given the total source flux of 3.4 Tmol S yr−1, the residence time of sulfate is 11.5 Myr.

3.3 Global oxygen cycling

The global O2 budget for our reference state is shown in Fig. 17. The simulated O2 inventory in the ocean-atmosphere system is 38×1018 mol (atmosphere =38×1018 mol, ocean =0.23×1018 mol). Organic carbon burial represents a major O2 source flux (marine = 11.3 Tmol O2 equiv. yr−1 and terrigenous = 3 Tmol O2 equiv. yr−1). Pyrite burial and hydrogen escape to space contribute 2.6 and 0.001 Tmol O2 equiv. yr−1, respectively. Given the total source–sink flux of 16.9 Tmol O2 yr−1, the residence time of O2 in the ocean–atmosphere system of our reference state is estimated as 2.26 Myr, which is consistent with previous estimates of 2–4 Myr (Berner, 1989; Berner, 2004a; Garrels and Perry, 1974).

Figure 17Schematics of global redox (O2) budget for the reference run. Arrows represent the O2 flux in terms of 1012 mol O2 equiv. yr−1. PAL stands for present atmospheric level. Pmol is 1015 mol. Emol is 1018 mol. ORG represents sedimentary organic carbon. PYR represents sedimentary pyrite sulfur. NPPocn is oceanic net primary production. EX stands for export production. MP is oceanic phosphate inventory. SO4 represents the concentration of sulfate.


4 Sensitivity experiment

Based on the results obtained above, we conclude that the CANOPS-GRB model is sufficient to describe basic biogeochemical characteristics in the modern ocean–atmosphere system. As a next step, we assess the dynamic response of the full model by performing sensitivity experiments with respect to P availability in surface environments.

4.1 Dynamic response to changes in P weathering

Here, we conduct a sensitivity experiment with respect to the P weathering rate in order to see how the atmospheric and oceanic O2 levels respond to changes in P availability in the exogenic system over a wide range of timescales. Specifically, we performed four simulations, varying the values of fP in Eq. (2) over 2 orders of magnitude. The reference state presented in the previous section is assumed for the initial condition, and the full model is allowed to evolve freely for 3 billion model years. These experiments demonstrate how P availability in surface environments affects global biogeochemical cycles and redox states of the atmosphere and oceans over a diverse range of timescales.

The simulated transient response is shown in Fig. 18. As expected, lower P availability leads to a lower oceanic P inventory (Fig. 18a), resulting in suppressed biological productivity in the ocean (Fig. 18b). Given the residence time of P in the ocean (18 kyr, see Sect. 3.2.3), these responses occur within 105 years. The suppressed biological productivity leads a decline of burial rate of organic matter in sediments (Fig. 18c). Specifically, 10 % and 1 % of fP give rise to the burial rate of marine Corg of 1 and 0.13 Tmol C yr−1 at 105 years, respectively (cf. the reference value of 11.3 Tmol C yr−1).

On the timescales of 105–106 years, the system reaches a quasi-steady state, but there is still a large redox imbalance due mainly to the suppression of Corg burial (Fig. 18i). This gives rise to deoxygenation of the atmosphere on a timescale of millions of years (Fig. 18d). Note that once the ocean interior becomes anoxic, the enhanced P recycling and preservation of organic matter in anoxic marine sediments tend to buffer the atmospheric deoxygenation (Fig. 18a and c). However, these passive responses do not alter the fundamental behavior: lower P availability results in lower atmospheric O2 levels. After atmospheric deoxygenation (>4 Myr), the system again reaches its quasi-steady state. Specifically, fP values of 10 % and 1 % result in atmospheric O2 levels of 9 % PAL and 0.6 % PAL, respectively.

Figure 18Biogeochemical responses obtained from the CANOPS-GRB model with different values of P availability, fP. (a) Oceanic phosphate inventory, MP. (b) Oceanic net primary production (NPPocn). (c) Burial rate of organic carbon (Corg) in marine sediments. (d) Atmospheric partial pressure of O2. PAL stands for present atmospheric level. (e) Atmospheric CH4 mixing ratio. (f) Sulfate concentration in the surface ocean layer. (g) Sedimentary reservoir size of pyrite sulfur (PYR). (h) Sedimentary reservoir size of organic carbon (ORG). (i) Global redox budget (GRB). For the fP=1 % run (red line), the calculation stopped when the atmospheric O2 level decreased to 10-5 PAL due to the numerical instability.


The following change is driven by the response of oceanic S cycle, which is characterized by the long residence time of 11.5 Myr (see Sect. 3.2.5). Ocean anoxia promotes the MSR and subsequent precipitation of pyrite in the ocean interior. However, our model demonstrates that the decline of seawater SO42- on a timescale of tens of millions of years is small (Fig. 18f) because the rate of MSR depends not only on the oceanic redox state but on the availability of organic matter for the MSR. The significant reduction of seawater SO42- occurs on the longer timescales (>100 Myr) for extremely low fP scenarios (0.016 and 0.01), in which atmospheric O2 levels decrease to <1 % PAL. These scenarios also accompany a growth of sedimentary S from gypsum to pyrite (Fig. 18g).

On longer timescales, sedimentary reservoirs affect the redox state of the atmosphere and oceans. The present result demonstrates that fP of 1 % finally leads to the catastrophic decrease in atmospheric O2 level at around 0.9 billion years (Fig. 18d). The simulation was stopped at this point due to the numerical instability. For other scenarios, the system reaches a new steady state after roughly 3 billion model years.

Biogenic CH4 production tends to be enhanced in anoxic oceans. However, the present result demonstrates that CH4 degassing to the atmosphere is inhibited by both limited availability of organic matter for methanogenesis and the anaerobic oxidation of CH4 by SO42-. Once seawater [SO42-] decreases below 1 mM, CH4 can escape from oceans to the atmosphere, promoting the buildup of CH4 in the atmosphere. Nevertheless, because of limited biological activity, atmospheric CH4 levels remain comparable to the modern value (∼1 ppmv) (Fig. 18e).

4.2 O2 budget for the less-oxygenated scenario

Figure 19 shows the O2 budget of the less-oxygenated state (fP=1.6 % scenario). Because P availability exerts a primary control on biospheric O2 production, the strongly suppressed P delivery to the ocean leads to low oceanic P levels and commensurately low biological productivity (0.08 Pmol and 1.3 Gt C yr−1, respectively). As a consequence, the atmospheric O2 level is low (0.75 % PAL). In this scenario the ocean interior is globally anoxic, and the preservation of organic C in marine sediments is enhanced. However, the suppressed biological productivity results in a low overall burial rate of organic C (0.9 Tmol O2 equiv. yr−1; ∼9 % of the reference value). When combined with the burial rate of terrigenous organic C, total O2 production by the organic C sub-cycle is 0.97 Tmol O2 equiv. yr−1. This O2 source is balanced by the sum of oxidative weathering and metamorphism. The role of the pyrite S sub-cycle in the global redox budget is also shown in Fig. 19. Most of the SO42- entering the anoxic ocean is buried as pyrite, representing a major O2 source (2.64 Tmol O2 equiv. yr−1). This O2 source is balanced by oxidation of sedimentary pyrite S through weathering (1.33 Tmol O2 equiv. yr−1) and metamorphism (1.31 Tmol O2 equiv. yr−1). In other words, the O2 budget for the weakly oxygenated Earth system is largely affected by the crustal S sub-cycle. This is in marked contrast to the well-oxygenated Earth system, where the O2 budget is mainly controlled by organic C sub-cycle.

The present result demonstrates that low atmospheric O2 states (∼1 % PAL) can be achieved in scenarios where the availability of P is strongly limited. However, a slight decrease of fP to 1 % leads to the destabilization of global O2 budget, providing implications for the stability and evolution of atmospheric O2 levels during the Proterozoic. This point will be further systematically examined in future work.

Figure 19Schematics of global redox (O2) budget for the fP=1.6 % (=10-1.8) run. Arrows represent the O2 flux in terms of 1012 mol O2 equiv. yr−1. PAL stands for present atmospheric level. Pmol is 1015 mol. Emol is 1018 mol. Gt C is 1015 g C. ORG represents sedimentary organic C. PYR represents sedimentary pyrite S. NPPocn is oceanic net primary production. EX stands for export production. MP is oceanic phosphate inventory. SO4 represents the concentration of sulfate.


5 Discussion

The reference run under the present condition demonstrates generally good agreement with modern observations (Sect. 3). The water circulation scheme provides an adequate representation of general ocean circulation, resulting in robust and reliable tracer distributions that are comparable to the modern observations. This provides a mechanistic foundation for simulating generalized ocean biogeochemical cycles. The ocean biogeochemistry module includes a series of biogeochemical processes in oxic–anoxic–sulfidic environments. The reference run gives rise to the distributions of nutrients and dissolved O2 that capture fundamental properties observed in the modern ocean. Integrated biogeochemical fluxes of the global ocean, such as biological productivity, material flow in the water column, and burial into sediments are also consistent with observational data. In our reference run the estimated organic carbon burial and oxidative weathering fluxes are relatively high compared to some previous estimates, though there remains significant uncertainty in globally integrated organic carbon weathering and burial fluxes. Further work will also be needed to better quantify the biogeochemical cycling in the continental shelf, which is a major locus of organic matter burial. In addition, some future developments to the N cycle may be needed, especially with regard to denitrification (e.g., anammox, coastal benthic denitrification). Despite these remaining challenges, our biogeochemical model is adequate for representing the general property of the coupled C–N–P–O2–S cycles.

A new scheme for oxidative weathering of organic matter and pyrite sulfur, mass balance calculation of O2 in the atmosphere, and time evolution of sedimentary reservoirs are explicitly included in the CANOPS-GRB model. These are a significant improvement from the previous versions of CANOPS (Lenton, 2020). The simplified framework for the global O2 budget is also useful to understand the response of complex biogeochemical systems. The computational efficiency of our CANOPS-GRB model allows us to conduct simulations over billions of model years with reasonable wall times (on the order of weeks), providing a useful tool for exploring the wide range of topics about the oxygenation history of Earth's atmosphere.

Sensitivity experiments with respect to the terrestrial weathering rate of P were conducted in order to see how the redox state of the ocean–atmosphere system responds to varying P availability in the surface system (Sect. 4). The CANOPS-GRB model appears to adequately simulate the biogeochemical dynamics over a wide range of timescales and is applicable for quantitative assessment of the evolution and stability of Earth's O2 cycling. Perhaps even more importantly, our results encourage us to perform further systematic examinations with Earth system models that have different complexities. Such an “Earth system model intercomparison” would be a critical step towards better mechanistic understanding of the stability and dynamics of atmospheric O2 levels over Earth's history.

Due to a lack of explicit Fe cycling and anaerobic metabolisms (such as anoxygenic photosynthesis), the current version of the model cannot be applicable for the simulation under the Archean-like weakly oxygenated (pO2<10-5 PAL) conditions. These topics are left to future studies, but it would be an achievable goal (Ozaki et al., 2018, 2019b; Van De Velde et al., 2021). The model design presented here also ignores the interaction between the surface system and the mantle (e.g., subduction) except for the degassing of reducing gases from the mantle. We note, however, that the surface–mantle interaction would have exerted a primary control on the long-term redox budget of Earth's surface system through the Earth's history (Canfield, 2004; Eguchi et al., 2020; Hayes and Waldbauer, 2006) and may be important for the discussion about the distant future (Ozaki and Reinhard, 2021). The importance of mantle and solid Earth controls on surficial environments is a crucially important topic for future research.

The CANOPS-GRB model has the basic capability to simulate the time evolution of the abundance of atmospheric biosignature gases (O2 and CH4) on a wide range of timescales. While the biogeochemical model is based on process studies to the extent possible, many processes are derived from empirical calibrations to Earth-like planets around Sun-like stars. Clearly, some these parameterizations, such as the photochemical parameterization among O2–O3–CH4 (Eq. 107), must be modified when applying the model to a range of habitable Earth-like exoplanets.

6 Conclusions

A new Earth system box model was developed – CANOPS-GRB v1.0. The new code release provides an improved description of the coupled C–N–P–O2–S biogeochemical cycles in the ocean–atmosphere–crust system, which can be utilized to examine the dynamics and stability of Earth's O2 cycle over a wide range of timescales. The computational efficiency and simple model design of CANOPS-GRB make it relatively easy to modify existing processes or add entirely new processes and components. CANOPS-GRB is thus a new and uniquely flexible tool capable of providing a coherent mechanistic framework for quantifying the biogeochemical cycles regulating Earth's O2 cycle. CANOPS-GRB is also a useful tool for the development of more comprehensive, low- to intermediate-complexity Earth system box models of biogeochemistry.

CANOPS-GRB provides an important step forward when coupled to new and existing geochemical proxy data. The accumulating geological and geochemical records have led to a new hypothesis for the evolution of atmospheric O2 levels on Earth. CANOPS-GRB was designed to facilitate simulation of a wide range of past conditions so as to permit more explicit testing of hypothesis about the function of biogeochemical cycles and their effect on the redox budget throughout Earth's history. Through the model–data synergy, CANOPS-GRB has great potential to provide an integrated, quantitative, and statistically informative picture of biogeochemical states, opening new perspectives on a wide range of scientific questions in research seeking to understand Earth's chemical evolution and the cause-and-effect relationships with evolving biosphere in particular.

CANOPS-GRB also provides significant steps forward in our predictive understanding of the links between geology, biogeochemistry, and the evolution of Earth's biosphere. It will allow for a fundamentally new and more precise quantitative understanding of evolving atmospheric biosignatures (O2, O3, CH4) on Earth and will broaden the interpretive power of Earth system evolution in the search for life beyond our planet. Additional elaboration of the CANOPS-GRB code could represent an important avenue for developing a more robust tool for diagnosing atmospheric biosignatures for future analysis of extrasolar worlds. In summary, it is anticipated that CANOPS-GRB will have many applications for problems linking the coupled evolution of life and the atmosphere on Earth and habitable rocky exoplanets.

Code availability

The bulk of CANOPS-GRB is written in Fortran as a standalone model. The model code can be found at GitHub (; Ozaki, 2022). This model is still undergoing regular development, and it is recommended that potential users contact the corresponding author (Kazumi Ozaki; to obtain the latest version.

Data availability

The GLODAPv2.2019 merged and adjusted data product used in Fig. 11 can be found at NOAA NCEI under (Olsen et al., 2019). These data and ancillary information are also available via their web pages (last access: 18 April 2021) and (last access: 18 April 2021).

Author contributions

KO designed the study, wrote the code, and ran model simulations. DBC and CTR contributed to code debugging. KO wrote the manuscript, with input from DBC, CTR, and ET.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


This work was supported by JSPS KAKENHI. Christopher T. Reinhard acknowledges the NASA Exobiology Program and the NASA Interdisciplinary Consortia for Astrobiology Research (ICAR). We are grateful to two anonymous reviewers for their insightful comments.

Financial support

This research has been supported by the JSPS KAKENHI (grant nos. JP16K05618, 19K21055, JP20K04066, and 22K18738), the NASA Astrobiology Institute (grant nos. 80NSSC19K0461, 80NSSC21K0594).

Review statement

This paper was edited by Paul Halloran and reviewed by two anonymous referees.


Alcott, L. J., Mills, B. J. W., and Poulton, S. W.: Stepwise Earth oxygenation is an inherent property of global biogeochemical cycling, Science, 366, 1333–1337,, 2019. 

Algeo, T. J. and Ingall, E.: Sedimentary Corg:P ratios, paleocean ventilation, and Phanerozoic atmospheric pO2, Palaeogeogr. Palaeocl., 256, 130–155,, 2007. 

Anderson, L. D., Delaney, M. L., and Faul, K. L.: Carbon to phosphorus ratios in sediments: Implications for nutrient cycling, Global Biogeochem. Cycles, 15, 65–79,, 2001. 

Archer, D., Kheshgi, H., and Maier-Reimer, E.: Dynamics of fossil fuel CO2 neutralization by marine CaCO3, Global Biogeochem. Cycles, 12, 259–276,, 1998. 

Archer, D. E., Eshel, G., Winguth, A., Broecker, W., Pierrehumbert, R., Tobis, M., and Jacob, R.: Atmospheric pCO2 sensitivity to the biological pump in the ocean, Global Biogeochem. Cycles, 14, 1219–1230,, 2000. 

Archer, D. E., Morford, J. L., and Emerson, S. R.: A model of suboxic sedimentary diagenesis suitable for automatic tuning and gridded global domains, Global Biogeochem. Cycles, 16, 17-11–17-21,, 2002. 

Armstrong, R. A., Lee, C., Hedges, J. I., Honjo, S., and Wakeham, S. G.: A new, mechanistic model for organic carbon fluxes in the ocean based on the quantitative association of POC with ballast minerals, Deep-Sea Res. Pt. II, 49, 219–236,, 2001. 

Arndt, S., Regnier, P., Goddéris, Y., and Donnadieu, Y.: GEOCLIM reloaded (v 1.0): a new coupled earth system model for past climate change, Geosci. Model Dev., 4, 451–481,, 2011. 

Baturin, G. N.: Issue of the relationship between primary productivity of organic carbon in ocean and phosphate accumulation (Holocene-Late Jurassic), Lithol. Min. Resour., 42, 318–348,, 2007. 

Beal, E. J., Claire, M. W., and House, C. H.: High rates of anaerobic methanotrophy at low sulfate concentrations with implications for past and present methane levels, Geobiology, 9, 131–139,, 2011. 

Belcher, C. M. and McElwain, J. C.: Limits for combustion in low O2 redefine paleoatmospheric predictions for the Mesozoic, Science, 321, 1197–1200,, 2008. 

Bellefroid, E. J., Hood, A. v. S., Hoffman, P. F., Thomas, M. D., Reinhard, C. T., and Planavsky, N. J.: Constraints on Paleoproterozoic atmospheric oxygen levels, P. Natl Acad. Sci. USA, 115, 8104–8109,, 2018. 

Benitez-Nelson, C. R.: The biogeochemical cycling of phosphorus in marine systems, Earth-Sci. Rev., 51, 109–135,, 2000. 

Berelson, W. M.: Particle settling rates increase with depth in the ocean, Deep-Sea Res. Pt. II, 49, 237–251,, 2001a. 

Berelson, W. M.: The Flux of Particulate Organic Carbon Into the Ocean Interior: A Comparison of Four U.S. JGOFS Regional Studies, Oceanography, 14, 59–67, 2001b. 

Berelson, W. M., Balch, W. M., Najjar, R., Feely, R. A., Sabine, C., and Lee, K.: Relating estimates of CaCO3 production, export, and dissolution in the water column to measurements of CaCO3 rain into sediment traps and dissolution on the sea floor: A revised global carbonate budget, Global Biogeochem. Cycles, 21, GB1024,, 2007. 

Bergman, N. M., Lenton, T. M., and Watson, A. J.: COPSE: A new model of biogeochemical cycling over Phanerozoic time, Am. J. Sci., 304, 397–437,, 2004. 

Berner, R. A.: Early diagenesis: A theoretical approach, Princeton University Press, Princeton, 256 pp., ISBN 0-691-08258-8, 1980. 

Berner, R. A.: Burial of organic carbon and pyrite sulfur in the modern ocean; its geochemical and environmental significance, Am. J. Sci., 282, 451–473,, 1982. 

Berner, R. A.: Biogeochemical cycles of carbon and sulfur and their effect on atmospheric oxygen over phanerozoic time, Palaeogeogr. Palaeocl., 75, 97–122,, 1989. 

Berner, R. A.: The Phanerozoic Carbon Cycle: CO2 and O2, Oxford University Press, ISBN 0-19-517333-3, 2004a. 

Berner, R. A.: A model for calcium, magnesium and sulfate in seawater over Phanerozoic time, Am. J. Sci., 304, 438–453,, 2004b. 

Berner, R. A.: GEOCARBSULF: A combined model for Phanerozoic atmospheric O2 and CO2, Geochim. Cosmochim. Ac., 70, 5653–5664,, 2006. 

Berner, R. A.: Phanerozoic atmospheric oxygen: New results using the GEOCARBSULF model, Am. J. Sci., 309, 603–606,, 2009. 

Berner, R. A. and Canfield, D. E.: A new model for atmospheric oxygen over Phanerozoic time, Am. J. Sci., 289, 333–361,, 1989. 

Berner, R. A. and Westrich, J. T.: Bioturbation and the early diagenesis of carbon and sulfur, Am. J. Sci., 285, 193–206,, 1985. 

Betts, J. N. and Holland, H. D.: The oxygen content of ocean bottom waters, the burial efficiency of organic carbon, and the regulation of atmospheric oxygen, Palaeogeogr. Palaeocl., 97, 5–18,, 1991. 

Bohlen, L., Dale, A. W., and Wallmann, K.: Simple transfer functions for calculating benthic fixed nitrogen losses and C:N:P regeneration ratios in global biogeochemical models, Global Biogeochem. Cycles, 26, GB3029,, 2012. 

Bolton, E. W., Berner, R. A., and Petsch, S. T.: The Weathering of Sedimentary Organic Matter as a Control on Atmospheric O2: II. Theoretical Modeling, Am. J. Sci., 306, 575–615,, 2006. 

Bottrell, S. H. and Newton, R. J.: Reconstruction of changes in global sulfur cycling from marine sulfate isotopes, Earth-Sci. Rev., 75, 59–83,, 2006. 

Boudreau, B. P.: A method-of-lines code for carbon and nutrient diagenesis in aquatic sediments, Comput. Geosci., 22, 479–496,, 1996. 

Bowles, M. W., Mogollón, J. M., Kasten, S., Zabel, M., and Hinrichs, K.-U.: Global rates of marine sulfate reduction and implications for sub–sea-floor metabolic activities, Science, 344, 889–891,, 2014. 

Bradley, J. A., Arndt, S., Amend, J. P., Burwicz, E., Dale, A. W., Egger, M., and LaRowe, D. E.: Widespread energy limitation to life in global subseafloor sediments, Sci. Adv., 6, eaba0697,, 2020. 

Brandes, J. A. and Devol, A. H.: A global marine-fixed nitrogen isotopic budget: Implications for Holocene nitrogen cycling, Global Biogeochem. Cycles, 16, GB001856,, 2002. 

Broecker, W. S. and Peng, T.-H.: Tracers in the sea, Eldigio Pr, New York, 690 pp., ISBN 9993186724, 1982. 

Burdige, D. J.: Burial of terrestrial organic matter in marine sediments: A re-assessment, Global Biogeochem. Cycles, 19, GB4011,, 2005. 

Burdige, D. J.: Preservation of Organic Matter in Marine Sediments: Controls, Mechanisms, and an Imbalance in Sediment Organic Carbon Budgets?, Chem. Rev., 107, 467–485,, 2007. 

Canfield, D. E.: Sulfate reduction and oxic respiration in marine sediments: implications for organic carbon preservation in euxinic environments, Deep-Sea Res. Pt. A., 36, 121–138,, 1989. 

Canfield, D. E.: Sulfate reduction in deep-sea sediments, Am. J. Sci., 291, 177–188,, 1991. 

Canfield, D. E.: Organic Matter Oxidation in Marine Sediments, in: Interactions of C, N, P and S Biogeochemical Cycles and Global Change, edited by: Wollast, R., Mackenzie, F. T., and Chou, L., Springer Berlin Heidelberg, Berlin, 333–363, ISBN 978-3-642-76066-2, 1993. 

Canfield, D. E.: The evolution of the Earth surface sulfur reservoir, Am. J. Sci., 304, 839–861,, 2004. 

Canfield, D. E. and Farquhar, J.: Animal evolution, bioturbation, and the sulfate concentration of the oceans, P. Natl Acad. Sci. USA, 106, 8123–8127,, 2009. 

Canfield, D. E., Zhang, S., Frank, A. B., Wang, X., Wang, H., Su, J., Ye, Y., and Frei, R.: Highly fractionated chromium isotopes in Mesoproterozoic-aged shales and atmospheric oxygen, Nat. Commun., 9, 2871,, 2018. 

Carr, M.-E., Friedrichs, M. A. M., Schmeltz, M., Noguchi Aita, M., Antoine, D., Arrigo, K. R., Asanuma, I., Aumont, O., Barber, R., Behrenfeld, M., Bidigare, R., Buitenhuis, E. T., Campbell, J., Ciotti, A., Dierssen, H., Dowell, M., Dunne, J., Esaias, W., Gentili, B., Gregg, W., Groom, S., Hoepffner, N., Ishizaka, J., Kameda, T., Le Quéré, C., Lohrenz, S., Marra, J., Mélin, F., Moore, K., Morel, A., Reddy, T. E., Ryan, J., Scardi, M., Smyth, T., Turpie, K., Tilstone, G., Waters, K., and Yamanaka, Y.: A comparison of global estimates of marine primary production from ocean color, Deep-Sea Res. Pt. II, 53, 741–770,, 2006. 

Catling, D. C. and Kasting, J. F.: Atmospheric Evolution on Inhabited and Lifeless Worlds, Cambridge University Press, ISBN 978-0-521-84412-3, 2017. 

Catling, D. C. and Zahnle, K. J.: The Archean atmosphere, Sci. Adv., 6, eaax1420,, 2020. 

Cha, H. J., Lee, C. B., Kim, B. S., Choi, M. S., and Ruttenberg, K. C.: Early diagenetic redistribution and burial of phosphorus in the sediments of the southwestern East Sea (Japan Sea), Marine Geol., 216, 127–143,, 2005. 

Claire, M. W., Catling, D. C., and Zahnle, K. J.: Biogeochemical modelling of the rise in atmospheric oxygen, Geobiology, 4, 239–269,, 2006. 

Cole, D. B., Reinhard, C. T., Wang, X., Gueguen, B., Halverson, G. P., Gibson, T., Hodgskiss, M. S. W., McKenzie, N. R., Lyons, T. W., and Planavsky, N. J.: A shale-hosted Cr isotope record of low atmospheric oxygen during the Proterozoic, Geology, 44, 555–558,, 2016. 

Cole, D. B., Ozaki, K., and Reinhard, C. T.: Atmospheric Oxygen Abundance, Marine Nutrient Availability, and Organic Carbon Fluxes to the Seafloor, Global Biogeochem. Cycles, 36, e2021GB007052,, 2022. 

Colman, A. S. and Holland, H. D.: The global diagenetic flux of phosphorus from marine sediments to the oceans: redox sensitivity and the control of atmosphreic oxygen levels, in: Marine authigenesis: from global to microbial, edited by: Glenn, C. R., Prevot-Lucas, L., and Lucas, J., SEPM (Society for Sedimentary Geology), 53–75, ISBN 1-56576-064-6, 2000. 

Compton, J., Mallinson, D., Glenn, C. R., Filippelli, G., Follmi, K., Shields, G. A., and Zanin, Y.: Variations in the global phosphorus cycle, in: Marine authigenesis: from global to microbial, edited by: Glenn, C. R., Prevot-Lucas, L., and Lucas, J., SEPM (Society for Sedimentary Geology), 21–33, 2000. 

Crichton, K. A., Wilson, J. D., Ridgwell, A., and Pearson, P. N.: Calibration of temperature-dependent ocean microbial processes in the cGENIE.muffin (v0.9.13) Earth system model, Geosci. Model Dev., 14, 125–149,, 2021. 

Crockford, P. W., Hayles, J. A., Bao, H., Planavsky, N. J., Bekker, A., Fralick, P. W., Halverson, G. P., Bui, T. H., Peng, Y., and Wing, B. A.: Triple oxygen isotope evidence for limited mid-Proterozoic primary productivity, Nature, 559, 613–616,, 2018. 

Daines, S. J., Mills, B. J. W., and Lenton, T. M.: Atmospheric oxygen regulation at low Proterozoic levels by incomplete oxidative weathering of sedimentary organic carbon, Nat. Commun., 8, 14379,, 2017. 

Dale, A. W., Meyers, S. R., Aguilera, D. R., Arndt, S., and Wallmann, K.: Controls on organic carbon and molybdenum accumulation in Cretaceous marine sediments from the Cenomanian–Turonian interval including Oceanic Anoxic Event 2, Chem. Geol., 324–325, 28–45,, 2012. 

Delaney, M. L.: Phosphorus accumulation in marine sediments and the oceanic phosphorus cycle, Global Biogeochem. Cycles, 12, 563–572,, 1998. 

Dellwig, O., Leipe, T., März, C., Glockzin, M., Pollehne, F., Schnetger, B., Yakushev, E. V., Böttcher, M. E., and Brumsack, H.-J.: A new particulate Mn–Fe–P-shuttle at the redoxcline of anoxic basins, Geochim. Cosmochim. Ac., 74, 7100–7115,, 2010. 

Derry, L. A.: Causes and consequences of mid-Proterozoic anoxia, Geophys. Res. Lett., 42, 2015GL065333,, 2015. 

Des Marais, D. J., Harwit, M. O., Jucks, K. W., Kasting, J. F., Lin, D. N., Lunine, J. I., Schneider, J., Seager, S., Traub, W. A., and Woolf, N. J.: Remote Sensing of Planetary Properties and Biosignatures on Extrasolar Terrestrial Planets, Astrobiology, 2, 153–181,, 2002. 

Deutsch, C., Sarmiento, J. L., Sigman, D. M., Gruber, N., and Dunne, J. P.: Spatial coupling of nitrogen inputs and losses in the ocean, Nature, 445, 163,, 2007. 

Devol, A. and Christensen, J. P.: Benthic fluxes and nitrogen cycling in sediments of the continental margin of the eastern North Pacific, J. Marine Res., 51, 345–372, 1993. 

Devol, A. H.: Denitrification, Anammox, and N2 Production in Marine Sediments, Ann. Rev. Mar. Sci., 7, 403–423,, 2015. 

DeVries, T., Deutsch, C., Primeau, F., Chang, B., and Devol, A.: Global rates of water-column denitrification derived from nitrogen gas measurements, Nat. Geosci., 5, 547,, 2012. 

DeVries, T., Deutsch, C., Rafter, P. A., and Primeau, F.: Marine denitrification rates determined from a global 3-D inverse model, Biogeosciences, 10, 2481–2496,, 2013. 

Doney, S. C., Lindsay, K., Caldeira, K., Campin, J. M., Drange, H., Dutay, J. C., Follows, M., Gao, Y., Gnanadesikan, A., Gruber, N., Ishida, A., Joos, F., Madec, G., Maier-Reimer, E., Marshall, J. C., Matear, R. J., Monfray, P., Mouchet, A., Najjar, R., Orr, J. C., Plattner, G. K., Sarmiento, J., Schlitzer, R., Slater, R., Totterdell, I. J., Weirig, M. F., Yamanaka, Y., and Yool, A.: Evaluating global ocean carbon models: The importance of realistic physics, Global Biogeochem. Cycles, 18, GB3017,, 2004. 

Donis, D., McGinnis, D. F., Holtappels, M., Felden, J., and Wenzhofer, F.: Assessing benthic oxygen fluxes in oligotrophic deep sea sediments (HAUSGARTEN observatory), Deep-Sea Res. Pt. I, 111, 1–10,, 2016. 

Duce, R. A., LaRoche, J., Altieri, K., Arrigo, K. R., Baker, A. R., Capone, D. G., Cornell, S., Dentener, F., Galloway, J., Ganeshram, R. S., Geider, R. J., Jickells, T., Kuypers, M. M., Langlois, R., Liss, P. S., Liu, S. M., Middelburg, J. J., Moore, C. M., Nickovic, S., Oschlies, A., Pedersen, T., Prospero, J., Schlitzer, R., Seitzinger, S., Sorensen, L. L., Uematsu, M., Ulloa, O., Voss, M., Ward, B., and Zamora, L.: Impacts of Atmospheric Anthropogenic Nitrogen on the Open Ocean, Science, 320, 893–897,, 2008. 

Dunne, J. P., Sarmiento, J. L., and Gnanadesikan, A.: A synthesis of global particle export from the surface ocean and cycling through the ocean interior and on the seafloor, Global Biogeochem. Cycles, 21, GB4006,, 2007. 

Eguchi, J., Seales, J., and Dasgupta, R.: Great Oxidation and Lomagundi events linked by deep cycling and enhanced degassing of carbon, Nat. Geosci., 13, 71–76,, 2020. 

Etheridge, D. M., Steele, L. P., Francey, R. J., and Langenfelds, R. L.: Atmospheric methane between 1000 A.D. and present: Evidence of anthropogenic emissions and climatic variability, J. Geophys. Res., 103, 15979–15993,, 1998. 

Eugster, O. and Gruber, N.: A probabilistic estimate of global marine N-fixation and denitrification, Global Biogeochem. Cycles, 26, GB4013,, 2012. 

Fakhraee, M., Planavsky, N. J., and Reinhard, C. T.: The role of environmental factors in the long-term evolution of the marine biological pump, Nat. Geosci., 13, 812–816,, 2020. 

Fiebig, J., Woodland, A. B., D'Alessandro, W., and Püttmann, W.: Excess methane in continental hydrothermal emissions is abiogenic, Geology, 37, 495–498,, 2009. 

Filippelli, G. M.: Carbon and phosphorus cycling in anoxic sediments of the Saanich Inlet, British Columbia, Marine Geol., 174, 307–321,, 2001. 

Föllmi, K. B.: The phosphorus cycle, phosphogenesis and marine phosphate-rich deposits, Earth-Sci. Rev., 40, 55–124,, 1996. 

Fowler, D., Coyle, M., Skiba, U., Sutton, M. A., Cape, J. N., Reis, S., Sheppard, L. J., Jenkins, A., Grizzetti, B., Galloway, J. N., Vitousek, P., Leach, A., Bouwman, A. F., Butterbach-Bahl, K., Dentener, F., Stevenson, D., Amann, M., and Voss, M.: The global nitrogen cycle in the twenty-first century, Phil. Trans. R. Soc. B, 368, 20130164,, 2013. 

Francois, R., Honjo, S., Krishfield, R., and Manganini, S.: Factors controlling the flux of organic carbon to the bathypelagic zone of the ocean, Global Biogeochem. Cycles, 16, 1087,, 2002. 

Froelich, P. N., Klinkhammer, G. P., Bender, M. L., Luedtke, N. A., Heath, G. R., Cullen, D., Dauphin, P., Hammond, D., Hartman, B., and Maynard, V.: Early oxidation of organic matter in pelagic sediments of the eastern equatorial Atlantic: suboxic diagenesis, Geochim. Cosmochim. Ac., 43, 1075–1090,, 1979. 

Galbraith, E. D. and Martiny, A. C.: A simple nutrient-dependence mechanism for predicting the stoichiometry of marine ecosystems, P. Natl Acad. Sci. USA, 112, 8199–8204,, 2015. 

Galloway, J. N., Dentener, F. J., Capone, D. G., Boyer, E. W., Howarth, R. W., Seitzinger, S. P., Asner, G. P., Cleveland, C. C., Green, P. A., Holland, E. A., Karl, D. M., Michaels, A. F., Porter, J. H., Townsend, A. R., and Vöosmarty, C. J.: Nitrogen Cycles: Past, Present, and Future, Biogeochemistry, 70, 153–226,, 2004. 

Garcia, H. E. and Gordon, L. I.: Oxygen solubility in seawater: Better fitting equations, Limnol. Oceanogr., 37, 1307–1312,, 1992. 

Garrels, R. M. and Lerman, A.: Phanerozoic cycles of sedimentary carbon and sulfur, P. Natl Acad. Sci. USA, 78, 4652–4656, 1981. 

Garrels, R. M. and Perry, J. E. A.: Cycling of carbon, sulfur, and oxygen through geologic time, The Sea, Wiley-Interscience, New York, edited by: Goldberg, E. D., 303–336, ISBN 067401734X, 1974. 

Goldblatt, C., Lenton, T. M., and Watson, A. J.: Bistability of atmospheric oxygen and the Great Oxidation, Nature, 443, 683–686, 2006. 

Graham, W. F. and Duce, R. A.: Atmospheric pathways of the phosphorus cycle, Geochim. Cosmochim. Ac., 43, 1195–1208,, 1979. 

Großkopf, T., Mohr, W., Baustian, T., Schunck, H., Gill, D., Kuypers, M. M. M., Lavik, G., Schmitz, R. A., Wallace, D. W. R., and LaRoche, J.: Doubling of marine dinitrogen-fixation rates based on direct measurements, Nature, 488, 361,, 2012. 

Gruber, N.: Chapter 1 – The Marine Nitrogen Cycle: Overview and Challenges, in: Nitrogen in the Marine Environment, 2nd edn., Academic Press, San Diego, 1–50,, 2008. 

Gruber, N. and Sarmiento, J. L.: Global patterns of marine nitrogen fixation and denitrification, Global Biogeochem. Cycles, 11, 235–266,, 1997. 

Gruber, N. and Sarmiento, J. L.: Biogeochemical/physical interactions in elemental cycles, in: THE SEA: Biological-Physical Interactions in the Oceans, edited by: Robinson, A. R., McCarthy, J. J., and Rothschild, B. J., John Wiley and Sons, New York, 337–399, 2002. 

Guidry, M. W., Mackenzie, F. T., and Arvidson, R. S.: Role of tectonics in phosphorus distribution and cycling, in: Marine Authigenesis: From Global to Microbial, edited by: Glenn, C. R., Prevot-Lucas, L., and Lucas, J., SEPM, 35–51, 2000. 

Gundersen, J. K. and Jorgensen, B. B.: Microstructure of diffusive boundary layers and the oxygen uptake of the sea floor, Nature, 345, 604,, 1990. 

Halevy, I., Peters, S. E., and Fischer, W. W.: Sulfate Burial Constraints on the Phanerozoic Sulfur Cycle, Science, 337, 331–334,, 2012. 

Handoh, I. C. and Lenton, T. M.: Periodic mid-Cretaceous oceanic anoxic events linked by oscillations of the phosphorus and oxygen biogeochemical cycles, Global Biogeochem. Cycles, 17, 1092,, 2003. 

Hartnett, H. E. and Devol, A. H.: Role of a strong oxygen-deficient zone in the preservation and degradation of organic matter: a carbon budget for the continental margins of northwest Mexico and Washington State, Geochim. Cosmochim. Ac., 67, 247–264,, 2003. 

Hartnett, H. E., Keil, R. G., Hedges, J. I., and Devol, A. H.: Influence of oxygen exposure time on organic carbon preservation in continental margin sediments, Nature, 391, 572–575, 1998. 

Hayes, C. T., Costa, K. M., Anderson, R. F., Calvo, E., Chase, Z., Demina, L. L., Dutay, J.-C., German, C. R., Heimbürger-Boavida, L.-E., Jaccard, S. L., Jacobel, A., Kohfeld, K. E., Kravchishina, M. D., Lippold, J., Mekik, F., Missiaen, L., Pavia, F. J., Paytan, A., Pedrosa-Pamies, R., Petrova, M. V., Rahman, S., Robinson, L. F., Roy-Barman, M., Sanchez-Vidal, A., Shiller, A., Tagliabue, A., Tessin, A. C., van Hulten, M., and Zhang, J.: Global Ocean Sediment Composition and Burial Flux in the Deep Sea, Global Biogeochem. Cycles, 35, e2020GB006769,, 2021. 

Hayes, J. M. and Waldbauer, J. R.: The carbon cycle and associated redox processes through time, Phil. Trans. R. Soc. B, 361, 931–950,, 2006. 

Hedges, J. I., Hu, F. S., Devol, A. H., Hartnett, H. E., Tsamakis, E., and Keil, R. G.: Sedimentary organic matter preservation; a test for selective degradation under oxic conditions, Am. J. Sci., 299, 529–555,, 1999. 

Heinze, C., Kriest, I., and Maier-Reimer, E.: Age offsets among different biogenic and lithogenic components of sediment cores revealed by numerical modeling, Paleoceanography, 24, PA4214,, 2009. 

Henrichs, S. M. and Reeburgh, W. S.: Anaerobic mineralization of marine sediment organic matter: Rates and the role of anaerobic processes in the oceanic carbon economy, Geomicrobiol. J., 5, 191–237,, 1987. 

Hensen, C., Landenberger, H., Zabel, M., and Schulz, H. D.: Quantification of diffusive benthic fluxes of nitrate, phosphate, and silicate in the southern Atlantic Ocean, Global Biogeochem. Cycles, 12, 193–210,, 1998. 

Hitchcock, D. R. and Lovelock, J. E.: Life detection by atmospheric analysis, Icarus, 7, 149–159,, 1967. 

Holland, H. D.: The Chemistry of the Atmosphere and Oceans, John Wiley & Sons, New York, ISBN 0471035092, 1978. 

Holser, W. T., Maynard, J. B., and Cruikshank, K. M.: Modelling the natural cycle of sulphur through Phanerozoic time, in: Evolution of the Global Biogeochemical Sulphur Cycle, edited by: Brimblecombe, P., and Lein, A. Y., John Wiley & Sons Ltd, New York, 21–56, 1989. 

Honjo, S.: Material fluxes and modes of sedimentation in the mesopelagic and bathypelagic zones, J. Marine Res., 38, 53–97, 1980. 

Honjo, S. and Manganini, S. J.: Annual biogenic particle fluxes to the interior of the North Atlantic Ocean; studied at 34 N 21 W and 48 N 21 W, Deep-Sea Res. Pt. II, 40, 587–607,, 1993. 

Hotinski, R. M., Kump, L. R., and Najjar, R. G.: Opening Pandora's Box: The impact of open system modeling on interpretations of anoxia, Paleoceanography, 15, 267–279,, 2000. 

Hyacinthe, C., Anschutz, P., Carbonel, P., Jouanneau, J. M., and Jorissen, F. J.: Early diagenetic processes in the muddy sediments of the Bay of Biscay, Marine Geol., 177, 111–128,, 2001. 

Ingall, E. and Jahnke, R.: Evidence for enhanced phosphorus regeneration from marine sediments overlain by oxygen depleted waters, Geochim. Cosmochim. Ac., 58, 2571–2575,, 1994. 

Ingall, E. and Jahnke, R.: Influence of water-column anoxia on the elemental fractionation of carbon and phosphorus during sediment diagenesis, Marine Geol., 139, 219–229,, 1997. 

Ingall, E. D. and Cappellen, P. V.: Relation between sedimentation rate and burial of organic phosphorus and organic carbon in marine sediments, Geochim. Cosmochim. Ac., 54, 373–386,, 1990. 

Ingall, E. D., Bustin, R. M., and Van Cappellen, P.: Influence of water column anoxia on the burial and preservation of carbon and phosphorus in marine shales, Geochim. Cosmochim. Ac., 57, 303–316,, 1993. 

Ittekkot, V.: The abiotically driven biological pump in the ocean and short-term fluctuations in atmospheric CO2 contents, Global Planet. Change, 8, 17–25,, 1993. 

Jahnke, R. A.: The global ocean flux of particulate organic carbon: Areal distribution and magnitude, Global Biogeochem. Cycles, 10, 71–88,, 1996. 

Joos, F., Sarmiento, J. L., and Siegenthaler, U.: Estimates of the effect of Southern Ocean iron fertilization on atmospheric CO2 concentrations, Nature, 349, 772–775,, 1991. 

Jørgensen, B. B.: Mineralization of organic matter in the sea bed—the role of sulphate reduction, Nature, 296, 643,, 1982. 

Jørgensen, B. B. and Kasten, S.: Sulfur cycling and methane oxidation, in: Marine Geochemistry, edited by: Schulz, H. D. and Zabel, M., Springer Berlin Heidelberg, 271–309,, 2006. 

Kagoshima, T., Sano, Y., Takahata, N., Maruoka, T., Fischer, T. P., and Hattori, K.: Sulphur geodynamic cycle, Sci. Rep.-UK, 5, 8330,, 2015. 

Karl, D., Michaels, A., Bergman, B., Capone, D., Carpenter, E., Letelier, R., Lipschultz, F., Paerl, H., Sigman, D., and Stal, L.: Dinitrogen fixation in the world's oceans, in: The Nitrogen Cycle at Regional to Global Scales, edited by: Boyer, E. W., and Howarth, R. W., Springer Netherlands, Dordrecht, 47–98,, 2002. 

Karl, D. M., Beversdorf, L., Björkman, K. M., Church, M. J., Martinez, A., and Delong, E. F.: Aerobic production of methane in the sea, Nat. Geosci., 1, 473–478,, 2008. 

Karthäuser, C., Ahmerkamp, S., Marchant, H. K., Bristow, L. A., Hauss, H., Iversen, M. H., Kiko, R., Maerz, J., Lavik, G., and Kuypers, M. M. M.: Small sinking particles control anammox rates in the Peruvian oxygen minimum zone, Nat. Commun., 12, 3235,, 2021. 

Kashiyama, Y., Ozaki, K., and Tajika, E.: Impact of the Evolution of Carbonate Ballasts on Marine Biogeochemistry in the Mesozoic and Associated Changes in Energy Delivery to Subsurface Waters, Paleontol. Res., 15, 89–99,, 2011. 

Katsev, S. and Crowe, S. A.: Organic carbon burial efficiencies in sediments: The power law of mineralization revisited, Geology, 43, 607–610,, 2015. 

Key, R. M., Olsen, A., van Heuven, S., Lauvset, S. K., Velo, A., Lin, X., Schirnick, C., Kozyr, A., Tanhua, T., Hoppema, M., Jutterström, S., Steinfeldt, R., Jeansson, E., Ishii, M., Perez, F. F., and Suzuki, T.: Global Ocean Data Analysis Project, Version 2 (GLODAPv2), NDP093_GLODAPv2, 2015. 

Kharecha, P., Kasting, J., and Siefert, J.: A coupled atmosphere–ecosystem model of the early Archean Earth, Geobiology, 3, 53–76,, 2005. 

Klaas, C. and Archer, D. E.: Association of sinking organic matter with various types of mineral ballast in the deep sea: Implications for the rain ratio, Global Biogeochem. Cycles, 16, 63-61–63-14,, 2002. 

Knox, F. and McElroy, M. B.: Changes in atmospheric CO2: Influence of the marine biota at high latitude, J. Geophys. Res., 89, 4629–4637,, 1984. 

Krissansen-Totton, J., Garland, R., Irwin, P., and Catling, D. C.: Detectability of Biosignatures in Anoxic Atmospheres with the James Webb Space Telescope: A TRAPPIST-1e Case Study, The Astronom. J., 156, 114,, 2018. 

Kump, L. R.: Chemical stability of the atmosphere and ocean, Palaeogeogr. Palaeocl., 75, 123–136,, 1989. 

Kump, L. R.: The rise of atmospheric oxygen, Nature, 451, 277–278,, 2008. 

Kuypers, M. M. M., Lavik, G., Woebken, D., Schmid, M., Fuchs, B. M., Amann, R., Jørgensen, B. B., and Jetten, M. S. M.: Massive nitrogen loss from the Benguela upwelling system through anaerobic ammonium oxidation, P. Natl. Acad. Sci. USA, 102, 6478–6483,, 2005. 

Kuznetsov, I., Neumann, T., and Burchard, H.: Model study on the ecosystem impact of a variable C:N:P ratio for cyanobacteria in the Baltic Proper, Ecol. Model., 219, 107–114,, 2008. 

Laakso, T. A. and Schrag, D. P.: Regulation of atmospheric oxygen during the Proterozoic, Earth Planet. Sc. Lett., 388, 81–91,, 2014. 

Larsson, U., Hajdu, S., Walve, J., and Elmgren, R.: Baltic Sea nitrogen fixation estimated from the summer increase in upper mixed layer total nitrogen, Limnol. Oceanogr., 46, 811–820,, 2001. 

Lasaga, A. C.: A new approach to isotopic modeling of the variation of atmospheric oxygen through the Phanerozoic, Am. J. Sci., 289, 411–435,, 1989. 

Lasaga, A. C. and Ohmoto, H.: The oxygen geochemical cycle: dynamics and stability, Geochim. Cosmochim. Ac., 66, 361–381,, 2002. 

Laws, E. A., Falkowski, P. G., Smith, W. O., Ducklow, H., and McCarthy, J. J.: Temperature effects on export production in the open ocean, Global Biogeochem. Cycles, 14, 1231–1246,, 2000. 

Ledwell, J. R., Watson, A. J., and Law, C. S.: Mixing of a tracer in the pycnocline, J. Geophys. Res., 103, 21499–21529,, 1998. 

Lenton, T. M.: Fire Feedbacks on Atmospheric Oxygen, in: Fire Phenomena and the Earth System, edited by: Belcher, C. M., 289–308,, 2013. 

Lenton, T. M.: On the use of models in understanding the rise of complex life, Interface Focus, 10, 20200018,, 2020. 

Lenton, T. M. and Watson, A. J.: Redfield revisited: 1. Regulation of nitrate, phosphate, and oxygen in the ocean, Global Biogeochem. Cycles, 14, 225–248,, 2000a. 

Lenton, T. M. and Watson, A. J.: Redfield revisited: 2. What regulates the oxygen content of the atmosphere?, Global Biogeochem. Cycles, 14, 249–268,, 2000b. 

Lenton, T. M., Daines, S. J., and Mills, B. J. W.: COPSE reloaded: An improved model of biogeochemical cycling over Phanerozoic time, Earth-Sci. Rev., 178, 1–28,, 2018. 

Lenton, T. M., Dahl, T. W., Daines, S. J., Mills, B. J. W., Ozaki, K., Saltzman, M. R., and Porada, P.: Earliest land plants created modern levels of atmospheric oxygen, P. Natl Acad. Sci. USA, 113, 9704–9709,, 2016. 

Lin, S. and Morse, J. W.: Sulfate reduction and iron sulfide mineral formation in Gulf of Mexico anoxic sediments, Am. J. Sci., 291, 55–89,, 1991. 

Liss, P. S. and Slater, P. G.: Flux of Gases across the Air-Sea Interface, Nature, 247, 181–184, 1974. 

Lord, N. S., Ridgwell, A., Thorne, M. C., and Lunt, D. J.: An impulse response function for the “long tail” of excess atmospheric CO2 in an Earth system model, Global Biogeochem. Cycles, 30, 2–17,, 2016. 

Lovelock, J. E.: A Physical Basis for Life Detection Experiments, Nature, 207, 568–570,, 1965. 

Lovelock, J. E.: Gaia as seen through the atmosphere, Atmos. Environ., 6, 579–580,, 1972. 

Lovelock, J. E.: Thermodynamics and the recognition of alien biospheres, P. Roy. Soc. Lond. B, 189, 167–181,, 1975. 

Lumpkin, R. and Speer, K.: Global Ocean Meridional Overturning, J. Phys. Oceanogr., 37, 2550–2562,, 2007. 

Luo, Y.-W., Doney, S. C., Anderson, L. A., Benavides, M., Berman-Frank, I., Bode, A., Bonnet, S., Boström, K. H., Böttjer, D., Capone, D. G., Carpenter, E. J., Chen, Y. L., Church, M. J., Dore, J. E., Falcón, L. I., Fernández, A., Foster, R. A., Furuya, K., Gómez, F., Gundersen, K., Hynes, A. M., Karl, D. M., Kitajima, S., Langlois, R. J., LaRoche, J., Letelier, R. M., Marañón, E., McGillicuddy Jr., D. J., Moisander, P. H., Moore, C. M., Mouriño-Carballido, B., Mulholland, M. R., Needoba, J. A., Orcutt, K. M., Poulton, A. J., Rahav, E., Raimbault, P., Rees, A. P., Riemann, L., Shiozaki, T., Subramaniam, A., Tyrrell, T., Turk-Kubo, K. A., Varela, M., Villareal, T. A., Webb, E. A., White, A. E., Wu, J., and Zehr, J. P.: Database of diazotrophs in global ocean: abundance, biomass and nitrogen fixation rates, Earth Syst. Sci. Data, 4, 47–73,, 2012. 

Lutz, M., Dunbar, R., and Caldeira, K.: Regional variability in the vertical flux of particulate organic carbon in the ocean interior, Global Biogeochem. Cycles, 16, 11-11–11-18,, 2002. 

Lyons, T. W. and Gill, B. C.: Ancient Sulfur Cycling and Oxygenation of the Early Biosphere, Elements, 6, 93–99,, 2010. 

Lyons, T. W., Reinhard, C. T., and Planavsky, N. J.: The rise of oxygen in Earth's early ocean and atmosphere, Nature, 506, 307–315,, 2014. 

Mackenzie, F. T., Ver, L. M., Sabine, C., Lane, M., and Lerman, A.: C, N, P, S Global Biogeochemical Cycles and Modeling of Global Change, in: Interactions of C, N, P and S Biogeochemical Cycles and Global Change, edited by: Wollast, R., Mackenzie, F. T., and Chou, L., Springer Berlin Heidelberg, Berlin, Heidelberg, 1–61, 1993. 

Maier-Reimer, E.: Geochemical cycles in an ocean general circulation model. Preindustrial tracer distributions, Global Biogeochem. Cycles, 7, 645–677,, 1993. 

Markovic, S., Paytan, A., and Wortmann, U. G.: Pleistocene sediment offloading and the global sulfur cycle, Biogeosciences, 12, 3043–3060,, 2015. 

Martin, J. H., Knauer, G. A., Karl, D. M., and Broenkow, W. W.: VERTEX: carbon cycling in the northeast Pacific, Deep-Sea Res. Pt. A, 34, 267–285,, 1987. 

Martin, W. R. and Sayles, F. L.: The Recycling of Biogenic Material at the Sea Floor, in: Treatise on Geochemistry (Second Edition), edited by: Turekian, K. K., Elsevier, Oxford, 33–59,, 2014. 

Mayor, M. and Queloz, D.: A Jupiter-mass companion to a solar-type star, Nature, 378, 355–359,, 1995. 

McManus, J., Berelson, W. M., Coale, K. H., Johnson, K. S., and Kilgore, T. E.: Phosphorus regeneration in continental margin sediments, Geochim. Cosmochim. Ac., 61, 2891–2907,, 1997. 

McManus, J., Berelson, W. M., Klinkhammer, G. P., Hammond, D. E., and Holm, C.: Authigenic uranium: Relationship to oxygen penetration depth and organic carbon rain, Geochim. Cosmochim. Ac., 69, 95–108,, 2005. 

Meadows, V. S.: Reflections on O2 as a Biosignature in Exoplanetary Atmospheres, Astrobiology, 17, 1022–1052,, 2017. 

Meadows, V. S., Reinhard, C. T., Arney, G. N., Parenteau, M. N., Schwieterman, E. W., Domagal-Goldman, S. D., Lincowski, A. P., Stapelfeldt, K. R., Rauer, H., DasSarma, S., Hegde, S., Narita, N., Deitrick, R., Lustig-Yaeger, J., Lyons, T. W., Siegler, N., and Grenfell, J. L.: Exoplanet Biosignatures: Understanding Oxygen as a Biosignature in the Context of Its Environment, Astrobiology, 18, 630–662,, 2018. 

Middelburg, J. J., Soetaert, K., Herman, P. M. J., and Heip, C. H. R.: Denitrification in marine sediments: A model study, Global Biogeochem. Cycles, 10, 661–673,, 1996. 

Middelburg, J. J., Soetaert, K., and Herman, P. M. J.: Empirical relationships for use in global diagenetic models, Deep-Sea Res. Pt. I, 44, 327–344,, 1997. 

Millero, F. J.: The oxidation of H2S in Black Sea waters, Deep-Sea Res. Pt. A, 38, S1139–S1150,, 1991. 

Millero, F. J.: Chemical Oceanography, 3rd edn., Taylor & Francis Group CRC Press, Boca Raton, 496 pp., 2006. 

Millero, F. J., Plese, T., and Fernandez, M.: The dissociation of hydrogen-sulfide in seawater, Limnol. Oceanogr., 33, 269–274, 1988. 

Morford, J. L. and Emerson, S.: The geochemistry of redox sensitive trace metals in sediments, Geochim. Cosmochim. Ac., 63, 1735–1750,, 1999. 

Muller-Karger, F. E., Varela, R., Thunell, R., Luerssen, R., Hu, C., and Walsh, J. J.: The importance of continental margins in the global carbon cycle, Geophys. Res. Lett., 32, L01602,, 2005. 

National Academies of Sciences, E. and Medicine: An Astrobiology Strategy for the Search for Life in the Universe, The National Academies Press, Washington, D.C., 188 pp.,, 2019. 

Nierop, K. G. J., Reichart, G.-J., Veld, H., and Sinninghe Damsté, J. S.: The influence of oxygen exposure time on the composition of macromolecular organic matter as revealed by surface sediments on the Murray Ridge (Arabian Sea), Geochim. Cosmochim. Ac., 206, 40–56,, 2017. 

Oguz, T., Ducklow, H. W., and Malanotte-Rizzoli, P.: Modeling distinct vertical biogeochemical structure of the Black Sea: Dynamical coupling of the oxic, suboxic, and anoxic layers, Global Biogeochem. Cycles, 14, 1331–1352,, 2000. 

Oguz, T., Murray, J. W., and Callahan, A. E.: Modeling redox cycling across the suboxic–anoxic interface zone in the Black Sea, Deep-Sea Res. Pt. I, 48, 761–787,, 2001. 

Olsen, A., Key, R. M., van Heuven, S., Lauvset, S. K., Velo, A., Lin, X., Schirnick, C., Kozyr, A., Tanhua, T., Hoppema, M., Jutterström, S., Steinfeldt, R., Jeansson, E., Ishii, M., Pérez, F. F., and Suzuki, T.: The Global Ocean Data Analysis Project version 2 (GLODAPv2) – an internally consistent data product for the world ocean, Earth Syst. Sci. Data, 8, 297–323,, 2016. 

Olsen, A., Lange, N., Key, R. M., Tanhua, T., Álvarez, M., Becker, S., Bittig, H. C., Carter, B. R., Cotrim da Cunha, L., Feely, R. A., van Heuven, S. M. A. C., Hoppema, M., Ishii, M., Jeansson, E., Jones, S. D., Jutterström, S., Karlsen, M. K., Kozyr, A., Lauvset, S. K., Lo Monaco, C., Murata, A., Pérez, F. F., Pfeil, B., Schirnick, C., Steinfeldt, R., Suzuki, T., Telszewski, M., Tilbrook, B., Velo, A., and Wanninkhof, R.: Global Ocean Data Analysis Project version 2.2019 (GLODAPv2.2019) (NCEI Accession 0186803), version 2.2019, NOAA National Centers for Environmental Information [data set],, 2019. 

Olson, S. L., Reinhard, C. T., and Lyons, T. W.: Limited role for methane in the mid-Proterozoic greenhouse, P. Natl Acad. Sci. USA, 113, 11447–11452,, 2016. 

Oschlies, A., Schulz, K. G., Riebesell, U., and Schmittner, A.: Simulated 21st century's increase in oceanic suboxia by CO2-enhanced biotic carbon export, Global Biogeochem. Cycles, 22, GB4008,, 2008. 

Ozaki, K.: kazumi-ozaki/CANOPS-GRBv1: CANOPS-GRBv1 (Version v1), Zenodo [code],, 2022. 

Ozaki, K. and Reinhard, C. T.: The future lifespan of Earth's oxygenated atmosphere, Nat. Geosci., 14, 138–142,, 2021. 

Ozaki, K. and Tajika, E.: Biogeochemical effects of atmospheric oxygen concentration, phosphorus weathering, and sea-level stand on oceanic redox chemistry: Implications for greenhouse climates, Earth Planet. Sc. Lett., 373, 129–139,, 2013. 

Ozaki, K., Tajima, S., and Tajika, E.: Conditions required for oceanic anoxia/euxinia: Constraints from a one-dimensional ocean biogeochemical cycle model, Earth Planet. Sc. Lett., 304, 270–279,, 2011. 

Ozaki, K., Tajika, E., Hong, P. K., Nakagawa, Y., and Reinhard, C. T.: Effects of primitive photosynthesis on Earth's early climate system, Nat. Geosci., 11, 55–59,, 2018. 

Ozaki, K., Reinhard, C. T., and Tajika, E.: A sluggish mid-Proterozoic biosphere and its effect on Earth's redox balance, Geobiology, 17, 3–11,, 2019a. 

Ozaki, K., Thompson, K. J., Simister, R. L., Crowe, S. A., and Reinhard, C. T.: Anoxygenic photosynthesis and the delayed oxygenation of Earth's atmosphere, Nat. Commun., 10, 3026,, 2019b. 

Pallud, C. and Van Cappellen, P.: Kinetics of microbial sulfate reduction in estuarine sediments, Geochim. Cosmochim. Ac., 70, 1148–1162,, 2006. 

Papadomanolaki, N. M., Lenstra, W. K., Wolthers, M., and Slomp, C. P.: Enhanced phosphorus recycling during past oceanic anoxia amplified by low rates of apatite authigenesis, Sci, Adv,, 8, eabn2370,, 2022. 

Petsch, S. T. and Berner, R. A.: Coupling the geochemical cycles of C, P, Fe, and S; the effect on atmospheric O2 and the isotopic records of carbon and sulfur, Am. J. Sci., 298, 246–262,, 1998. 

Petsch, S. T., Eglinton, T. I., and Edwards, K. J.: 14C-Dead Living Biomass: Evidence for Microbial Assimilation of Ancient Organic Carbon During Shale Weathering, Science, 292, 1127–1131,, 2001. 

Pfeifer, K., Hensen, C., Adler, M., Wenzhfer, F., Weber, B., and Schulz, H. D.: Modeling of subsurface calcite dissolution, including the respiration and reoxidation processes of marine sediments in the region of equatorial upwelling off Gabon, Geochim. Cosmochim. Ac., 66, 4247–4259,, 2002. 

Planavsky, N. J., Cole, D. B., Reinhard, C. T., Diamond, C., Love, G. D., Luo, G., Zhang, S., Konhauser, K. O., and Lyons, T. W.: No evidence for high atmospheric oxygen levels 1,400 million years ago, P. Natl Acad. Sci. USA, 113, E2550–E2551,, 2016. 

Planavsky, N. J., Cole, D. B., Isson, T. T., Reinhard, C. T., Crockford, P. W., Sheldon, N. D., and Lyons, T. W.: A case for low atmospheric oxygen levels during Earth's middle history, Emerging Topics in Life Sciences, 2, 149–159,, 2018. 

Prentice, I. C., Farquhar, G. D., Fasham, M. J. R., Goulden, M. L., Heimann, M., Jaramillo, V. J., Kheshgi, H. S., Le Quere, C., Scholes, R. J., and Wallace, D. W. R.: The carbon cycle and atmospheric carbon dioxide, in: Climate Change 2001: the Scientific Basis, edited by: Houghton, J. T., Ding, Y., Griggs, D. J., Noguer, N., van der Linden, P. J., Xiaosu, D., Maskell, K., and Johnson, C. A., Cambridge University Press, New York, 2001. 

Quigg, A., Finkel, Z. V., Irwin, A. J., Rosenthal, Y., Ho, T.-Y., Reinfelder, J. R., Schofield, O., Morel, F. M. M., and Falkowski, P. G.: The evolutionary inheritance of elemental stoichiometry in marine phytoplankton, Nature, 425, 291–294,, 2003. 

Raiswell, R. and Canfield, D. E.: The Iron Biogeochemical Cycle Past and Present, Geochemical Perspectives, 1, 1–2, 2012. 

Raynaud, D., Jouzel, J., Barnola, J. M., Chappellaz, J., Delmas, R. J., and Lorius, C.: The Ice Record of Greenhouse Gases, Science, 259, 926–934,, 1993. 

Redfield, A. C., Ketchum, B. H., and Richards, F. A.: The influence of organisms on the composition of sea-water, in: The Sea, edited by: Hill, M. N., Interscience Publishers, New York, 26–77, 1963. 

Reimers, C. E., Jahnke, R. A., and McCorkle, D. C.: Carbon fluxes and burial rates over the continental slope and rise off central California with implications for the global carbon cycle, Global Biogeochem. Cycles, 6, 199–224,, 1992. 

Reinhard, C. T., Olson, S. L., Schwieterman, E. W., and Lyons, T. W.: False Negatives for Remote Life Detection on Ocean-Bearing Planets: Lessons from the Early Earth, Astrobiology, 17, 287–297,, 2017a. 

Reinhard, C. T., Planavsky, N. J., Gill, B. C., Ozaki, K., Robbins, L. J., Lyons, T. W., Fischer, W. W., Wang, C., Cole, D. B., and Konhauser, K. O.: Evolution of the global phosphorus cycle, Nature, 541, 386–389,, 2017b. 

Reinhard, C. T., Olson, S. L., Kirtland Turner, S., Pälike, C., Kanzaki, Y., and Ridgwell, A.: Oceanic and atmospheric methane cycling in the cGENIE Earth system model – release v0.9.14, Geosci. Model Dev., 13, 5687–5706,, 2020. 

Ridgwell, A. and Hargreaves, J. C.: Regulation of atmospheric CO2 by deep-sea sediments in an Earth system model, Global Biogeochem. Cycles, 21, GB2008,, 2007. 

Romaniello, S. J. and Derry, L. A.: An intermediate-complexity model for simulating marine biogeochemistry in deep time: Validation against the modern global ocean, Geochem. Geophys. Geosyst., 11, Q08001,, 2010. 

Rowe, G. T., Morse, J., Nunnally, C., and Boland, G. S.: Sediment community oxygen consumption in the deep Gulf of Mexico, Deep-Sea Res. Pt. II, 55, 2686–2691,, 2008. 

Ruttenberg, K. C.: Reassessment of the oceanic residence time of phosphorus, Chem. Geol., 107, 405–409,, 1993. 

Ruttenberg, K. C.: The Global Phosphorus Cycle, in: Treatise on Geochemistry, edited by: Turekian, K. K., Pergamon, Oxford, 585–643,, 2003. 

Sachs, O., Sauter, E. J., Schlüter, M., Rutgers van der Loeff, M. M., Jerosch, K., and Holby, O.: Benthic organic carbon flux and oxygen penetration reflect different plankton provinces in the Southern Ocean, Deep-Sea Res. Pt. I, 56, 1319–1335,, 2009. 

Sagan, C., Thompson, W. R., Carlson, R., Gurnett, D., and Hord, C.: A search for life on Earth from the Galileo spacecraft, Nature, 365, 715–721, 1993. 

Sarmiento, J. L. and Gruber, N.: Ocean biogeochemical dynamics, Princeton University Press, ISBN 0-691-01707-7, 2006. 

Sarmiento, J. L. and Toggweiler, J. R.: A new model for the role of the oceans in determining atmospheric PCO2, Nature, 308, 621–624,, 1984. 

Schenau, S. J. and De Lange, G. J.: Phosphorus regeneration vs. burial in sediments of the Arabian Sea, Marine Chem., 75, 201–217,, 2001. 

Schlesinger, W. H. and Bernhardt, E. S.: The Global Cycles of Sulfur and Mercury, in: Biogeochemistry, 3rd edn., Academic Press, Boston, 469–486,, 2013. 

Schwieterman, E. W., Kiang, N. Y., Parenteau, M. N., Harman, C. E., DasSarma, S., Fischer, T. M., Arney, G. N., Hartnett, H. E., Reinhard, C. T., Olson, S. L., Meadows, V. S., Cockell, C. S., Walker, S. I., Grenfell, J. L., Hegde, S., Rugheimer, S., Hu, R., and Lyons, T. W.: Exoplanet Biosignatures: A Review of Remotely Detectable Signs of Life, Astrobiology, 18, 663–708,, 2018. 

Shaffer, G.: Phosphate pumps and shuttles in the Black Sea, Nature, 321, 515–517,, 1986. 

Shaffer, G. and Sarmiento, J. L.: Biogeochemical cycling in the global ocean: 1. A new, analytical model with continuous vertical resolution and high-latitude dynamics, J. Geophys. Res., 100, 2659–2672,, 1995. 

Shaffer, G., Malskær Olsen, S., and Pepke Pedersen, J. O.: Presentation, calibration and validation of the low-order, DCESS Earth System Model (Version 1), Geosci. Model Dev., 1, 17–51,, 2008. 

Sharoni, S. and Halevy, I.: Geologic controls on phytoplankton elemental composition, P. Natl Acad. Sci. USA, 119, e2113263118,, 2022. 

Siegenthaler, U. and Wenk, T.: Rapid atmospheric CO2 variations and ocean circulation, Nature, 308, 624–626,, 1984. 

Sleep, N. H.: Dioxygen over geological time, in: Metal ions in biological systems, edited by: Sigel, A., Sigel, H., and Sigel, R. K. O., Taylor & Francis Group, Boca Raton, 49–73, 2005. 

Slomp, C. P. and Van Cappellen, P.: The global marine phosphorus cycle: sensitivity to oceanic circulation, Biogeosciences, 4, 155–171,, 2007. 

Slomp, C. P., Thomson, J., and de Lange, G. J.: Enhanced regeneration of phosphorus during formation of the most recent eastern Mediterranean sapropel (S1), Geochim. Cosmochim. Ac., 66, 1171–1184,, 2002. 

Sloyan, B. M.: Spatial variability of mixing in the Southern Ocean, Geophys. Res. Lett., 32, L18603,, 2005. 

Soulet, G., Hilton, R. G., Garnett, M. H., Roylands, T., Klotz, S., Croissant, T., Dellinger, M., and Le Bouteiller, C.: Temperature control on CO2 emissions from the weathering of sedimentary rocks, Nat. Geosci., 14, 665–671,, 2021. 

Southam, J. R., Peterson, W. H., and Brass, G. W.: Dynamics of anoxia, Palaeogeogr. Palaeocl., 40, 183–198,, 1982. 

Steefel, C. I. and MacQuarrie, K. T. B.: Approaches to modeling of reactive transport in porous media, Rev. Miner. Geochem., 34, 85–129, 1996. 

Suess, E.: Particulate organic carbon flux in the oceans – surface productivity and oxygen utilization, Nature, 288, 260–263,, 1980. 

Tang, D., Shi, X., Wang, X., and Jiang, G.: Extremely low oxygen concentration in mid-Proterozoic shallow seawaters, Precambrian Res., 276, 145–157,, 2016. 

Tarhan, L. G., Droser, M. L., Planavsky, N. J., and Johnston, D. T.: Protracted development of bioturbation through the early Palaeozoic Era, Nat. Geosci., 8, 865,, 2015. 

Tarpgaard, I. H., Røy, H., and Jørgensen, B. B.: Concurrent low- and high-affinity sulfate reduction kinetics in marine sediment, Geochim. Cosmochim. Ac., 75, 2997–3010,, 2011. 

The LUVOIR Team: Mission Concept Study Final Report, in: arXiv e-prints,, 2019. 

Tostevin, R., Turchyn, A. V., Farquhar, J., Johnston, D. T., Eldridge, D. L., Bishop, J. K. B., and McIlvin, M.: Multiple sulfur isotope constraints on the modern sulfur cycle, Earth Planet. Sc. Lett., 396, 14–21,, 2014. 

Tromp, T. K., Van Cappellen, P., and Key, R. M.: A global model for the early diagenesis of organic carbon and organic phosphorus in marine sediments, Geochim. Cosmochim. Ac., 59, 1259–1284,, 1995. 

Tsunogai, S. and Noriki, S.: Particulate fluxes of carbonate and organic carbon in the ocean. Is the marine biological activity working as a sink of the atmospheric carbon?, Tellus B, 43, 265–266,, 1991. 

Turchyn, A. V. and Schrag, D. P.: Oxygen Isotope Constraints on the Sulfur Cycle over the Past 10 Million Years, Science, 303, 2004–2007,, 2004. 

Turchyn, A. V. and Schrag, D. P.: Cenozoic evolution of the sulfur cycle: Insight from oxygen isotopes in marine sulfate, Earth Planet. Sc. Lett., 241, 763–779,, 2006. 

Turnewitsch, R. and Pohl, C.: An estimate of the efficiency of the iron- and manganese-driven dissolved inorganic phosphorus trap at an oxic/euxinic water column redoxcline, Global Biogeochem. Cycles, 24, GB4025,, 2010. 

Tyrrell, T.: The relative influences of nitrogen and phosphorus on oceanic primary production, Nature, 400, 525–531,, 1999. 

Van Cappellen, P. and Ingall, E. D.: Benthic phosphorus regeneration, net primary production, and ocean anoxia: A model of the coupled marine biogeochemical cycles of carbon and phosphorus, Paleoceanography, 9, 677–692,, 1994. 

Van Cappellen, P. and Ingall, E. D.: Redox Stabilization of the Atmosphere and Oceans by Phosphorus-Limited Marine Productivity, Science, 271, 493–496,, 1996. 

Van Cappellen, P. and Wang, Y.: Cycling of iron and manganese in surface sediments; a general theory for the coupled transport and reaction of carbon, oxygen, nitrogen, sulfur, iron, and manganese, Am. J. Sci., 296, 197–243,, 1996. 

van de Velde, S. J., Hülse, D., Reinhard, C. T., and Ridgwell, A.: Iron and sulfur cycling in the cGENIE.muffin Earth system model (v0.9.21), Geosci. Model Dev., 14, 2713–2745,, 2021. 

Volk, T. and Hoffert, M. I.: Ocean carbon pumps: Analysis of relative strengths and efficiencies in ocean-driven atmospheric CO2 changes, in: The Carbon Cycle and Atmospheric CO2: Natural Variations Archean to Present, edited by: Sundquist, E. T. and Broecker, W. S., 99–110,, 1985. 

Walker, J. C. G.: Evolution of the atmosphere, Macmillan, New York, 318 pp., ISBN 0-02-854390-4, 1977. 

Walker, J. C. G. and Brimblecombe, P.: Iron and sulfur in the pre-biologic ocean, Precambrian Res., 28, 205–222,, 1985. 

Wallmann, K.: Feedbacks between oceanic redox states and marine productivity: A model perspective focused on benthic phosphorus cycling, Global Biogeochem. Cycles, 17, 1084,, 2003. 

Wallmann, K.: Phosphorus imbalance in the global ocean?, Global Biogeochem. Cycles, 24, GB4030,, 2010. 

Wang, W.-L., Moore, J. K., Martiny, A. C., and Primeau, F. W.: Convergent estimates of marine nitrogen fixation, Nature, 566, 205–211,, 2019. 

WebBook, N. C.: NIST Chemistry WebBook,, 2022. 

Westrich, J. T. and Berner, R. A.: The role of sedimentary organic matter in bacterial sulfate reduction: The G model tested, Limnol. Oceanogr., 29, 236–249,, 1984. 

Wheat, C. G., Feely, R. A., and Mottl, M. J.: Phosphate removal by oceanic hydrothermal processes: An update of the phosphorus budget in the oceans, Geochim. Cosmochim. Ac., 60, 3593–3608,, 1996. 

Wheat, C. G., McManus, J., Mottl, M. J., and Giambalvo, E.: Oceanic phosphorus imbalance: Magnitude of the mid-ocean ridge flank hydrothermal sink, Geophys. Res. Lett., 30, 1895,, 2003. 

Woodward, F. I.: Global primary production, Current Biology, 17, R269–R273,, 2007. 

Wortmann, U. G. and Paytan, A.: Rapid Variability of Seawater Chemistry Over the Past 130 Million Years, Science, 337, 334–336,, 2012. 

Yakushev, E. V. and Neretin, L. N.: One-dimensional modeling of nitrogen and sulfur cycles in the aphotic zones of the Black and Arabian Seas, Global Biogeochem. Cycles, 11, 401–414,, 1997. 

Yakushev, E. V., Pollehne, F., Jost, G., Kuznetsov, I., Schneider, B., and Umlauf, L.: Analysis of the water column oxic/anoxic interface in the Black and Baltic seas with a numerical model, Marine Chem., 107, 388–410,, 2007. 

Yamanaka, Y. and Tajika, E.: The role of the vertical fluxes of particulate organic matter and calcite in the oceanic carbon cycle: Studies using an ocean biogeochemical general circulation model, Global Biogeochem. Cycles, 10, 361–382,, 1996. 

Yao, W. and Millero, F.: The chemistry of the anoxic waters in the Framvaren Fjord, Norway, Aquatic Geochemistry, 1, 53–88,, 1995. 

Yaroshevsky, A. A.: Abundances of chemical elements in the Earth's crust, Geochem. Int., 44, 48–55,, 2006. 

Zabel, M., Dahmke, A., and Schulz, H. D.: Regional distribution of diffusive phosphate and silicate fluxes through the sediment–water interface: the eastern South Atlantic, Deep-Sea Res. Pt. I, 45, 277–300,, 1998. 

Zhang, S., Wang, X., Wang, H., Bjerrum, C. J., Hammarlund, E. U., Dahl, T. W., and Canfield, D. E.: Reply to Planavsky et al.: Strong evidence for high atmospheric oxygen levels 1,400 million years ago, P. Natl Acad. Sci. USA, 113, E2552–E2553,, 2016.  

Zhao, M., Zhang, S., Tarhan, L. G., Reinhard, C. T., and Planavsky, N.: The role of calcium in regulating marine phosphorus burial and atmospheric oxygenation, Nat. Commun., 11, 2232,, 2020. 

Short summary
A new biogeochemical model (CANOPS-GRB v1.0) for assessing the redox stability and dynamics of the ocean–atmosphere system on geologic timescales has been developed. In this paper, we present a full description of the model and its performance. CANOPS-GRB is a useful tool for understanding the factors regulating atmospheric O2 level and has the potential to greatly refine our current understanding of Earth's oxygenation history.