Ocean carbon and nitrogen isotopes in CSIRO Mk3L-COAL version 1.0: a tool for palaeoceanographic research

The isotopes of carbon (δ¹³C) and nitrogen (δ¹⁵N) are commonly used proxies for understanding the ocean. When used in tandem, they provide powerful insight into physical and biogeochemical processes. Here, we detail the implementation of δ¹³C and δ¹⁵N in the ocean component of an Earth system model. We evaluate our simulated δ¹³C and δ¹⁵N against contemporary measurements, place the model's performance alongside other isotope-enabled models and document the response of δ¹³C and δ¹⁵N to changes in ecosystem functioning. The model combines the Commonwealth Scientific and Industrial Research Organisation Mark 3L (CSIRO Mk3L) climate system model with the Carbon of the Ocean, Atmosphere and Land (COAL) biogeochemical model. The oceanic component of CSIRO Mk3L-COAL has a resolution of 1.6^∘ latitude × 2.8^∘ longitude and resolves multimillennial timescales, running at a rate of ∼400 years per day. We show that this coarse-resolution, computationally efficient model adequately reproduces water column and core-top δ¹³C and δ¹⁵N measurements, making it a useful tool for palaeoceanographic research. Changes to ecosystem function involve varying phytoplankton stoichiometry, varying CaCO₃ production based on calcite saturation state and varying N₂ fixation via iron limitation. We find that large changes in CaCO₃ production have little effect on δ¹³C and δ¹⁵N, while changes in N₂ fixation and phytoplankton stoichiometry have substantial and complex effects. Interpretations of palaeoceanographic records are therefore open to multiple lines of interpretation where multiple processes imprint on the isotopic signature, such as in the tropics, where denitrification, N₂ fixation and nutrient utilisation influence δ¹⁵N. Hence, there is significant scope for isotope-enabled models to provide more robust interpretations of the proxy records.

1 Introduction

Elements that are involved in reactions of interest, such as exchanges of carbon and nutrients, experience isotopic fractionation. Typically, the heavier isotope will be enriched in the reactant during kinetic fractionation, in more oxidised compounds during equilibrium fractionation and in the denser form during phase state fractionation (i.e. evaporation). Because fractionation against one isotope relative to the other is minuscule, the isotopic content of a sample is conventionally expressed as a δ value (δ^hE), where the ratio of the heavy to light element in solution (^hE:^lE) is compared to a standard ratio (^hE_std:^lE_std) in units of per mille (‰).

$$\begin{array}{}\text{(1)}& {\mathit{\delta }}^{\mathrm{h}}E=({\displaystyle \frac{{}^{\mathrm{h}}E:{}^{\mathrm{l}}E}{{}^{\mathrm{h}}{E}_{\mathrm{std}}:{}^{\mathrm{l}}{E}_{\mathrm{std}}}}-\mathrm{1})\cdot \mathrm{1000}\end{array}$$

The strength of fractionation against the heavier isotope during a given reaction, ϵ, is also expressed in per mille notation. Fractionation with an ϵ equal to 10 ‰, for example, will involve 990 units of ^hE for every 1000 units of ^lE at a hypothetical standard ratio (^hE_std:^lE_std) of 1:1. At more realistic standard ratios ⋘1:1, say 0.0112372:1 for a δ¹³C value of 0 ‰, a fractionation at 10 ‰ would involve $\sim \mathrm{0.0111123}(\mathrm{0.010}\cdot \frac{\mathrm{0.0112372}}{\mathrm{1.0112372}})$ units of ¹³C per unit of ¹²C. Slightly greater preference of one isotope over another in this case involves a preference for the lighter carbon isotope (¹²C) over the heavier (¹³C), which enriches the remaining dissolved inorganic carbon (DIC) in ¹³C and depletes the product. Certain isotopic preferences, or strengths of fractionation, therefore allow certain reactions to be detected in the environment.

The measurement of the stable isotopes of carbon (δ¹³C) and nitrogen (δ¹⁵N) have been fundamental for understanding how these important elements cycle within the ocean (e.g. Schmittner and Somes, 2016; Menviel et al., 2017a; Rafter et al., 2017; Muglia et al., 2018). We will now briefly introduce each isotope in turn.

The distribution of δ¹³C is dependent on air–sea gas exchange, ocean circulation and organic matter cycling. These contributions make the δ¹³C signature difficult to interpret, and several modelling studies have attempted to elucidate their roles (Tagliabue and Bopp, 2008; Schmittner et al., 2013). These studies have shown that preferential uptake of ¹²C over ¹³C by biology in surface waters enforces strong horizontal and vertical gradients in δ¹³C of DIC (δ¹³C_DIC), greatly enriching surface waters, particularly in subtropical gyres where vertical exchange with deeper waters is restricted (Tagliabue and Bopp, 2008; Schmittner et al., 2013). Meanwhile, air–sea gas exchange and carbon speciation control the δ¹³C_DIC reservoir over longer timescales (Schmittner et al., 2013). Because air–sea and speciation fractionation are temperature dependent, such that cooler conditions tend to elevate the δ¹³C_DIC of surface waters, they also tend to smooth the gradients produced by biology by working antagonistically to them. Despite this smoothing, biological fractionation drives strong gradients at the surface, which imparts unique δ¹³C signatures to the water masses that are carried into the interior. These insights have provided clear evidence of reduced ventilation rates in the deep ocean during glacial climates (Tagliabue et al., 2009; Menviel et al., 2017a; Muglia et al., 2018).

δ¹⁵N is determined by biological processes that add or remove fixed forms of nitrogen. It therefore records the relative rates of sources and sinks within the marine nitrogen cycle (Brandes and Devol, 2002). Dinitrogen (N₂) fixation is the largest source of fixed nitrogen to the ocean, the bulk of which occurs in warm, sunlit surface waters and introduces nitrogen with a δ¹⁵N of approximately −1 ‰ (Sigman et al., 2009). Denitrification is the largest sink of fixed nitrogen and occurs in deoxygenated water columns and sediments. Denitrification fractionates strongly against ¹⁵N at ∼25 ‰ (Cline and Kaplan, 1975). Fractionation during denitrification is most strongly expressed in the water column where ample nitrate (NO₃) is available, making water column denitrification responsible for elevating global mean δ¹⁵N above the −1 ‰ of N₂ fixers (Brandes and Devol, 2002). Meanwhile, denitrification occurring in the sediments only weakly fractionates against ¹⁵N (Sigman et al., 2009), providing only a slight enrichment of δ¹⁵N above that introduced by N₂ fixation. Variations in δ¹⁵N can therefore tell us about global changes in the ratio of sedimentary to water column denitrification, with increases in δ¹⁵N associated with increases in the proportion of denitrification occurring in the water column (Galbraith et al., 2013), but it can also reflect regional changes in N₂ fixation and denitrification (Ganeshram et al., 1995; Ren et al., 2009; Straub et al., 2013).

However, nitrogen isotopes are also subject to the effect of utilisation, which makes the interpretation of δ¹⁵N more complicated. Basically, when nitrogen is abundant, the preference for ¹⁴N over ¹⁵N increases but when nitrogen is limited this preference disappears (Altabet and Francois, 2001). Complete utilisation of nitrogen therefore reduces fractionation to 0 ‰. While this adds complexity, it also imbues δ¹⁵N as a proxy of nutrient utilisation by phytoplankton. As nitrogen supply to phytoplankton is controlled by physical delivery from below, changes in δ¹⁵N can be interpreted as changes in the physical supply (Studer et al., 2018). Phytoplankton fractionate against ¹⁵N at ∼5 ‰ (Wada, 1980) when bioavailable nitrogen is abundant. If nitrogen is utilised to completion, which occurs in much of the low to midlatitude ocean, then no fractionation will occur and the δ¹⁵N of organic matter will reflect the δ¹⁵N of the nitrogen that was supplied (Sigman et al., 2009). However, in the case where nitrogen is not consumed towards completion, which occurs in zones of strong upwelling/mixing near coastlines, the Equator and high latitudes, the bioavailable nitrogen pool will be enriched in ¹⁵N as phytoplankton preferentially consume ¹⁴N. As the remaining bioavailable N is continually enriched in ¹⁵N, the organic matter that settles into sediments beneath a zone of incomplete nutrient utilisation will bear this enriched δ¹⁵N signal. In combination with modelling (Schmittner and Somes, 2016), the δ¹⁵N record is able to provide evidence for a more efficient utilisation of bioavailable nitrogen during glacial times (Martinez-Garcia et al., 2014; Kemeny et al., 2018) and a less efficient one during the Holocene (Studer et al., 2018).

Complimentary measurements of δ¹³C and δ¹⁵N provide powerful, multi-focal insights into oceanographic processes. δ¹³C is largely a reflection of how water masses mix away the strong vertical and horizontal gradients enforced by biology, while δ¹⁵N simultaneously reflects changes in the major sources and sinks of the marine nitrogen cycle and how effectively nutrients are consumed at the surface. However, the interpretation of these isotopes is often difficult. They are subject to considerable uncertainty because there are multiple processes that imprint on the measured values. Our goal is to equip version 1.0 of the Commonwealth Scientific and Industrial Research Organisation Mark 3L (CSIRO Mk3L) climate system model with the Carbon of the Ocean, Atmosphere and Land (COAL) Earth system model with oceanic δ¹³C and δ¹⁵N such that this model can be used for interpreting palaeoceanographic records. First, we introduce CSIRO Mk3L-COAL. Second, we detail the equations that govern the implementation of carbon and nitrogen isotopes. Third, we assess our simulated isotopes against contemporary measurements from both the water column and sediments and compare the model performance against other isotope-enabled models. Finally, as a first test of the model, we take the opportunity to document how changes in ecosystem functioning affect δ¹³C and δ¹⁵N.

2 CSIRO Mk3L-COAL v1.0

The CSIRO Mk3L-COAL couples a computationally efficient climate system model (Phipps et al., 2013) with biogeochemical cycles in the ocean, atmosphere and land. The model is therefore based on the CSIRO Mk3L climate system model, where the “L” denotes that it is a low-resolution version of the CSIRO Mk3 model that contributed towards the third phase of Coupled Model Intercomparison Project (Meehl et al., 2007) and the Fourth Assessment Report of the Intergovernmental Panel on Climate Change (Solomon et al., 2007). See Smith (2007) for a complete discussion of the CSIRO family of climate models. The land biogeochemical component represents carbon, nitrogen and phosphorus cycles in the Community Atmosphere Biosphere Land Exchange (CABLE) (Mao et al., 2011). The ocean component currently represents carbon, alkalinity, oxygen, nitrogen, phosphorus and iron cycles. The atmospheric component conserves carbon and alters its radiative properties according to changes in its carbon content. For this paper, we focus on the ocean biogeochemical model (OBGCM).

Previous versions of the OBGCM have explored changes in oceanic properties under past (Buchanan et al., 2016), present (Buchanan et al., 2018) and future scenarios (Matear and Lenton, 2014, 2018). These studies demonstrate that the model can reproduce observed features of the global carbon cycle, nutrient cycling and organic matter cycling in the ocean. The OBGCM offers highly efficient simulations of these processes at computational speeds of ∼400 years per day when the ocean general circulation model (OGCM) is run offline (compared to ∼10 years per day in fully coupled mode). The ocean is made up of grid cells of 1.6^∘ in latitude by 2.8^∘ in longitude, with 21 vertical depth levels spaced by 25 m at the surface and 450 m in the deep ocean (Table 1). The OGCM time step is 1 h, while the OBGCM time step is 1 d. The ability of the OBGCM to reproduce large-scale dynamical and biogeochemical properties of the ocean coupled with its fast computational speed makes the OBGCM useful as a tool for palaeoceanographic research.

2.1 Ocean biogeochemical model (OBGCM)

The OBGCM is equipped with 13 prognostic tracers (Fig. 1). These can be grouped into carbon chemistry fields, oxygen fields, nutrient fields, age tracers and nitrous oxide (N₂O). Carbon chemistry fields include DIC, alkalinity (ALK), DI¹³C and radiocarbon (¹⁴C). Radiocarbon is simulated according to Toggweiler et al. (1989). Oxygen fields include dissolved oxygen (O₂) and abiotic dissolved oxygen (${\mathrm{O}}_{\mathrm{2}}^{\mathrm{abio}}$), a purely physical tracer from which true oxygen utilisation (TOU) can be calculated (Duteil et al., 2013). Nutrient fields include phosphate (PO₄), dissolved bioavailable iron (Fe), nitrate (NO₃) and ¹⁵NO₃. Although we define the phosphorus and nitrogen tracers as their dominant species, being PO₄ and NO₃, these tracers can also be thought of as total dissolved inorganic phosphorus and nitrogen pools. Remineralisation, for instance, implicitly accounts for the process of nitrification from ammonium (NH₄) to NO₃ (Paulmier et al., 2009) and therefore implicitly includes NH₄ and nitrite (NO₂) within the NO₃ tracer. Age tracer fields include years since subduction from the surface (Age_gbl) and years since entering a suboxic zone where O₂ concentrations are less than 10 mmol m⁻³ (Age_omz). Finally, N₂O in µmol m⁻³ is produced via nitrification and denitrification according to the temperature-dependent equations of Freing et al. (2012). All air–sea gas exchanges (CO₂, ¹³CO₂, O₂ and N₂O) and carbon speciation reactions are computed according to the Ocean Modelling Intercomparison Project phase 6 protocol (Orr et al., 2017).

Figure 1 A conceptual representation of the ocean biogeochemical model (OBGCM). The bottom panel shows organic matter cycling involving the isotopes of carbon and nitrogen. (1) Carbon chemistry reactions. (2) Air–sea gas exchange. (3) Biological uptake of nutrients and production of organic and inorganic matter. Particulate organic carbon (POC) is produced by the general phytoplankton group and N₂ fixers (diazotrophs), while particulate inorganic carbon (PIC) as calcium carbonate (CaCO₃) is produced by calcifiers. Export of POC by the general (G) phytoplankton group and N₂ fixers (D; diazotrophs) is herein referred to as C${}_{\mathrm{org}}^{\mathrm{G}}$ and C${}_{\mathrm{org}}^{\mathrm{D}}$ (see Sect. A1 in the Appendix), respectively. (4) Remineralisation of sinking organic matter under oxic and suboxic conditions. (5) Sedimentary oxic and suboxic remineralisation. (6) Nitrous oxide production and consumption.

Because the isotopes of carbon and nitrogen are influenced by biological processes and there is as yet no accepted standard for ecosystem model parameterisation in the community (see Hülse et al., 2017, for a more detailed discussion), we provide a thorough description of the ecosystem component of the OBGCM in Sect. A in the Appendix. Default parameters for the OBGCM are further provided in Sect. B in the Appendix. Briefly, the ecosystem model simulates the production, remineralisation and stoichiometry (elemental composition) of three types of primary producers: a general phytoplankton group, diazotrophs (N₂ fixers) and calcifiers.

3 Carbon and nitrogen isotope equations

3.1  δ¹³C

The OBGCM explicitly simulates the fractionation of ¹³C from the total DIC pool, where for simplicity we make the assumption that the total DIC pool represents the light isotope of carbon and is therefore DI¹²C. Fractionation occurs during air–sea gas exchange, equilibrium reactions and biological consumption in the euphotic zone.

The air–sea gas exchange of ¹³CO₂ is calculated as the exchange of CO₂ with additional fractionation factors applied to the sea–air and air–sea components (Zhang et al., 1995; Orr et al., 2017). The flux of ¹³CO₂ across the air–sea interface, F(¹³CO₂), therefore takes the form of CO₂ with additional terms that convert to units of ¹³C in both environments. Without any isotopic fractionation, the equation requires the gas piston velocity of carbon dioxide in m s⁻¹ ($k_{{\mathrm{CO}}_{\mathrm{2}}}$), the concentration of aqueous CO₂ in both mediums at the air–sea interface in mmol m⁻³ (${\mathrm{CO}}_{\mathrm{2}}^{\mathrm{air}}$ and ${\mathrm{CO}}_{\mathrm{2}}^{\mathrm{sea}}$) and the ratios of ¹³C:¹²C in both mediums (R_atm and R_sea):

$$\begin{array}{ll}\text{(2)}& {\displaystyle }& {\displaystyle }F\left({}^{\mathrm{13}}{\mathrm{CO}}_{\mathrm{2}}\right)=k_{{\mathrm{CO}}_{\mathrm{2}}}\cdot ({\mathrm{CO}}_{\mathrm{2}}^{\mathrm{air}}\cdot R_{\mathrm{atm}}-{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{sea}}\cdot R_{\mathrm{DIC}}),{\displaystyle }& {\displaystyle }\text{where}\\ {\displaystyle }& {\displaystyle }{R}_{\mathrm{DIC}}={\displaystyle \frac{{\mathrm{DI}}^{\mathrm{13}}\mathrm{C}}{{\mathrm{DI}}^{\mathrm{13}}\mathrm{C}+{\mathrm{DI}}^{\mathrm{12}}\mathrm{C}}}\\ {\displaystyle }& {\displaystyle }{R}_{\mathrm{atm}}={\displaystyle \frac{{}^{\mathrm{13}}{\mathrm{CO}}_{\mathrm{2}}}{{}^{\mathrm{13}}{\mathrm{CO}}_{\mathrm{2}}+{}^{\mathrm{12}}{\mathrm{CO}}_{\mathrm{2}}}}=\mathrm{0.011164381}.\end{array}$$

A transfer of ¹³C into the ocean is therefore positive and an outgassing is negative. The R_atm is set to a preindustrial atmospheric δ¹³C of −6.48 ‰ (Friedli et al., 1986).

The fractionation of carbon isotopes during air–sea exchange involves three components. These are (${\mathit{ϵ}}_{\mathrm{k}}^{{}^{\mathrm{13}}\mathrm{C}}$), a kinetic fractionation that occurs during transfer of gaseous CO₂ into or out of the ocean; (${\mathit{ϵ}}_{\mathrm{aq}\leftarrow \mathrm{g}}^{{}^{\mathrm{13}}\mathrm{C}}$), a fractionation that occurs as gaseous CO₂ becomes aqueous CO₂ (is dissolved in solution); and (${\mathit{ϵ}}_{\mathrm{DIC}\leftarrow \mathrm{g}}^{{}^{\mathrm{13}}\mathrm{C}}$), an equilibrium isotopic fractionation as carbon speciates into DIC constituents (H₂CO₂ ⇔ ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ ⇔ ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$). The kinetic fractionation during transfer, ${\mathit{ϵ}}_{\mathrm{k}}^{{}^{\mathrm{13}}\mathrm{C}}$, is constant at 0.99912, thus reducing the δ¹³C of carbon entering the ocean by 0.88 ‰. Conversely, carbon outgassing increases the δ¹³C of the ocean. The fractionations during dissolution (${\mathit{ϵ}}_{\mathrm{aq}\leftarrow \mathrm{g}}^{{}^{\mathrm{13}}\mathrm{C}}$) and speciation (${\mathit{ϵ}}_{\mathrm{DIC}\leftarrow \mathrm{g}}^{{}^{\mathrm{13}}\mathrm{C}}$) are both dependent on temperature. Fractionation during speciation is also dependent on the fraction of ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ relative to total DIC ($f_{{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}}$). These fractionation factors are parameterised as

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }{\mathit{ϵ}}_{\mathrm{aq}\leftarrow \mathrm{g}}^{{}^{\mathrm{13}}\mathrm{C}}={\displaystyle \frac{\mathrm{0.0049}\cdot T-\mathrm{1.31}}{\mathrm{1000}}}+\mathrm{1}\text{(4)}& {\displaystyle }& {\displaystyle }{\mathit{ϵ}}_{\mathrm{DIC}\leftarrow \mathrm{g}}^{{}^{\mathrm{13}}\mathrm{C}}={\displaystyle \frac{\mathrm{0.0144}\cdot T\cdot f_{{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}}-\mathrm{0.107}\cdot T+\mathrm{10.53}}{\mathrm{1000}}}+\mathrm{1}.\end{array}$$

Dissolution of CO₂ into seawater (${\mathit{ϵ}}_{\mathrm{aq}\leftarrow \mathrm{g}}^{{}^{\mathrm{13}}\mathrm{C}}$) therefore preferences the lighter isotope and lowers δ¹³C by between 1.32 ‰ and 1.14 ‰, while the speciation of gaseous CO₂ into DIC instead prefers the heavier isotope and raises δ¹³C by between 10.7 ‰ and 6.8 ‰ for temperatures between −2 and 35 ^∘C.

These fractionation factors are applied to the gaseous exchange of CO₂ (Eq. 2) to calculate carbon isotopic fractionation.

$$\begin{array}{ll}{\displaystyle }F\left({}^{\mathrm{13}}{\mathrm{CO}}_{\mathrm{2}}\right)=& {\displaystyle }k\cdot {\mathit{ϵ}}_{\mathrm{k}}^{{}^{\mathrm{13}}\mathrm{C}}\cdot {\mathit{ϵ}}_{\mathrm{aq}\leftarrow \mathrm{g}}^{{}^{\mathrm{13}}\mathrm{C}}\\ \text{(5)}& {\displaystyle }& {\displaystyle }\cdot ({\mathrm{CO}}_{\mathrm{2}}^{\mathrm{air}}\cdot R_{\mathrm{atm}}-{\displaystyle \frac{{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{sea}}\cdot R_{\mathrm{DIC}}}{{\mathit{ϵ}}_{\mathrm{DIC}\leftarrow \mathrm{g}}^{{}^{\mathrm{13}}\mathrm{C}}}})\end{array}$$

Because fractionation to aqueous CO₂ from DIC (${\mathit{ϵ}}_{\mathrm{aq}\leftarrow \mathrm{DIC}}^{{}^{\mathrm{13}}\mathrm{C}}$) is equal to $\frac{{\mathit{ϵ}}_{\mathrm{aq}\leftarrow \mathrm{g}}^{{}^{\mathrm{13}}\mathrm{C}}}{{\mathit{ϵ}}_{\mathrm{DIC}\leftarrow \mathrm{g}}^{{}^{\mathrm{13}}\mathrm{C}}}$, a strong preference to hold the heavy isotope in solution exists (${\mathit{ϵ}}_{\mathrm{aq}\leftarrow \mathrm{DIC}}^{{}^{\mathrm{13}}\mathrm{C}}=-\mathrm{11.9}$ ‰ to −7.9 ‰ between −2 and 35 ^∘C). Aqueous carbon that is transferred to the atmosphere is hence depleted in ¹³C. It is therefore the equilibrium fractionation associated with carbon speciation that is largely responsible for bolstering the oceanic δ¹³C signature above the atmospheric signature, as it tends to shift ¹³C towards the oxidised species (${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$), a tendency that strengthens under cooler conditions.

In the default version of CSIRO Mk3L-COAL v1.0, the fractionation of carbon during biological uptake (${\mathit{ϵ}}_{\mathrm{bio}}^{{}^{\mathrm{13}}\mathrm{C}}$) is set at 21 ‰ for general phytoplankton, 12 ‰ for diazotrophs (e.g. Carpenter et al., 1997) and at 2 ‰ for the formation of calcite (Ortiz et al., 1996). However, a variable fractionation rate for the general phytoplankton group may be activated and depends on the aqueous CO₂ concentration (CO₂(aq) in mmol m⁻³) and the growth rate (µ in d⁻¹, as a function of temperature and limiting resources) following Tagliabue and Bopp (2008):

$$\begin{array}{}\text{(6)}& {\mathit{ϵ}}_{\mathrm{bio}}^{{}^{\mathrm{13}}\mathrm{C}}=(\mathrm{0.371}-{\displaystyle \frac{\mathit{\mu }}{{\mathrm{CO}}_{\mathrm{2}\left(\mathrm{aq}\right)}}})/\mathrm{0.015}.\end{array}$$

An upper bound of 25 ‰ exists within Eq. (6) when $\frac{\mathit{\mu }}{{\mathrm{CO}}_{\mathrm{2}\left(\mathrm{aq}\right)}}$ approaches zero but a lower bound does not. We chose to limit ${\mathit{ϵ}}_{\mathrm{bio}}^{{}^{\mathrm{13}}\mathrm{C}}$ to a minimum of 15 ‰ given the reported variations of ${\mathit{ϵ}}_{\mathrm{bio}}^{{}^{\mathrm{13}}\mathrm{C}}$ from culture studies (e.g. Laws et al., 1995).

$$\begin{array}{}\text{(7)}& {\mathit{ϵ}}_{\mathrm{bio}}^{{}^{\mathrm{13}}\mathrm{C}}=max(\mathrm{15},{\mathit{ϵ}}_{\mathrm{bio}}^{{}^{\mathrm{13}}\mathrm{C}})\end{array}$$

Biological fractionation of ¹³C is then applied to the uptake and release of organic carbon.

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }{\mathrm{DI}}^{\mathrm{13}}\mathrm{C}=R_{\mathrm{DIC}}\cdot {\mathrm{C}}_{\mathrm{org}}\cdot (\mathrm{1}-{\displaystyle \frac{{\mathit{ϵ}}_{\mathrm{bio}}^{{}^{\mathrm{13}}\mathrm{C}}}{\mathrm{1000}}})\end{array}$$

Because biological fractionation is strong for the general phytoplankton group, which dominates export production throughout most of the ocean, this imparts a negative δ¹³C signature to the deep ocean. Subsequent remineralisation releases DIC with no fractionation. Finally, the concentration of DI¹³C is converted into a δ¹³C via

$$\begin{array}{}\text{(9)}& {\mathit{\delta }}^{\mathrm{13}}\mathrm{C}=({\displaystyle \frac{{\mathrm{DI}}^{\mathrm{13}}\mathrm{C}}{\mathrm{DIC}}}\cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{0.0112372}}}-\mathrm{1})\cdot \mathrm{1000},\end{array}$$

where 0.0112372 is the Pee Dee Belemnite standard (Craig, 1957).

3.2  δ¹⁵N

The OBGCM explicitly simulates the fractionation of ¹⁵N from the pool of bioavailable nitrogen. For simplicity, we treat this bioavailable pool as nitrate (NO₃), where total NO₃ is the sum of ¹⁵NO₃ and ¹⁴NO₃. We therefore chose to ignore fractionation during reactions involving ammonium, nitrite and dissolved organic nitrogen, which can vary in their isotopic composition independent of NO₃ but represent a small fraction of the bioavailable pool of nitrogen.

The isotopic signatures of N₂ fixation and atmospheric deposition, and the fractionation during water column denitrification (${\mathit{ϵ}}_{\mathrm{wc}}^{{}^{\mathrm{15}}\mathrm{N}}$) and sedimentary denitrification (${\mathit{ϵ}}_{\mathrm{sed}}^{{}^{\mathrm{15}}\mathrm{N}}$) determine the global δ¹⁵N of NO₃ (Brandes and Devol, 2002). Biological assimilation (${\mathit{ϵ}}_{\mathrm{bio}}^{{}^{\mathrm{15}}\mathrm{N}}$) and remineralisation are internal exchanges of the oceanic nitrogen cycle and affect the distribution of δ¹⁵NO₃. N₂ fixation and atmospheric deposition introduce ¹⁵NO₃ to the ocean with δ¹⁵N values of −1 ‰ and −2 ‰, respectively, while biological assimilation, water column denitrification and sedimentary denitrification fractionate against ¹⁵NO₃ at 5 ‰, 20 ‰ and 3 ‰, respectively (Sigman et al., 2009, Fig. 1).

The accepted standard ¹⁵N:¹⁴N ratio used to measure variations in nature is the average atmospheric ¹⁵N:¹⁴N ratio of 0.0036765. To minimise numerical errors caused by the OGCM, we set the atmospheric standard to 1. This scales up the ¹⁵NO₃ such that a δ¹⁵N value of 0 ‰ was equivalent to a ¹⁵N:¹⁴N ratio of 1:1.

Because we simulate NO₃ and ¹⁵NO₃ as tracers, our calculations require solving for an implicit pool of ¹⁴NO₃ during each reaction involving ¹⁵NO₃. The introduction of NO₃ at a fixed ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ of −1 ‰ due to remineralisation of N₂ fixer biomass provides a simple example with which we can begin to describe our equations. Setting the isotopic value of newly fixed NO₃ to −1 ‰ is simple because it removes any complications associated with fractionation. We note, however, that in reality the nitrogenase enzyme does fractionate during its conversion of aqueous N₂ (+0.7 ‰) to ammonium and the biomass that is subsequently produced can vary substantially depending of the type of nitrogenase enzyme used (vanadium versus molybdenum based) (McRose et al., 2019). However, we choose to implicitly account for these transformations and considerably simplify them by setting the δ¹⁵N of N₂ fixer biomass equal to −1 ‰, which reflects the biomass of N₂ fixers associated with the more common Mo-nitrogenase.

A ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ of −1 ‰ is equivalent to a ¹⁵N:¹⁴N ratio of 0.999 in our approach, where 0 ‰ equals a 1:1 ratio of ¹⁵N:¹⁴N. If the amount of NO₃ being added is known alongside its ¹⁵N:¹⁴N ratio, in this case 0.999 for N₂ fixation, we are able to calculate how much ¹⁵NO₃ is added. We begin with two equations that describe the system.

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{\mathrm{NO}}_{\mathrm{3}}={}^{\mathrm{15}}{\mathrm{NO}}_{\mathrm{3}}+{}^{\mathrm{14}}{\mathrm{NO}}_{\mathrm{3}}\text{(11)}& {\displaystyle }& {\displaystyle }{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}=({\displaystyle \frac{{}^{\mathrm{15}}{\mathrm{NO}}_{\mathrm{3}}{/}^{\mathrm{14}}{\mathrm{NO}}_{\mathrm{3}}}{{}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{std}}/{}^{\mathrm{14}}{\mathrm{N}}_{\mathrm{std}}}}-\mathrm{1})\cdot \mathrm{1000}\end{array}$$

Ultimately, we need to solve for the change in ¹⁵NO₃ associated with an introduction of NO₃ by N₂ fixation. Our known variables are the change in NO₃, the δ¹⁵N of that NO₃ and the ¹⁵N_std∕¹⁴N_std. Our two unknowns are ¹⁵NO₃ and ¹⁴NO₃. We must solve for ¹⁴NO₃ implicitly by describing it according to ¹⁵NO₃ by rearranging Eq. (11).

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }{}^{\mathrm{14}}{\mathrm{NO}}_{\mathrm{3}}={}^{\mathrm{15}}{\mathrm{NO}}_{\mathrm{3}}/\left(\right({\displaystyle \frac{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}}{\mathrm{1000}}}+\mathrm{1})\cdot {}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{std}}{/}^{\mathrm{14}}{\mathrm{N}}_{\mathrm{std}})\end{array}$$

This allows us to replace the ¹⁴NO₃ term in Eq. (10), such that

$$\begin{array}{}\text{(13)}& {\displaystyle }{\mathrm{NO}}_{\mathrm{3}}={}^{\mathrm{15}}{\mathrm{NO}}_{\mathrm{3}}+{}^{\mathrm{15}}{\mathrm{NO}}_{\mathrm{3}}/\left(\right({\displaystyle \frac{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}}{\mathrm{1000}}}+\mathrm{1})\cdot {}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{std}}{/}^{\mathrm{14}}{\mathrm{N}}_{\mathrm{std}}).\end{array}$$

In our example of N₂ fixation, we know the δ¹⁵N of the newly added NO₃ as being −1 ‰. We also know ¹⁵N_std∕¹⁴N_std as equal to 1:1, or 1. Our equation is simplified.

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle }{\mathrm{NO}}_{\mathrm{3}}={}^{\mathrm{15}}{\mathrm{NO}}_{\mathrm{3}}+{}^{\mathrm{15}}{\mathrm{NO}}_{\mathrm{3}}/\mathrm{0.999}\end{array}$$

We can now solve for ¹⁵NO₃ by rearranging the equation.

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }{}^{\mathrm{15}}{\mathrm{NO}}_{\mathrm{3}}={\displaystyle \frac{\mathrm{0.999}\cdot {\mathrm{NO}}_{\mathrm{3}}}{\mathrm{1}+\mathrm{0.999}}}\end{array}$$

The same calculation is applied to NO₃ addition via atmospheric deposition, except at a constant fraction of 0.998 (δ¹⁵N $=-\mathrm{2}$ ‰), and can be applied to any addition or subtraction of ¹⁵NO₃ relative to NO₃ where the isotopic signature is known.

Fractionating against ¹⁵NO₃ during biological assimilation (${\mathit{ϵ}}_{\mathrm{bio}}^{{}^{\mathrm{15}}\mathrm{N}}$), water column denitrification (${\mathit{ϵ}}_{\mathrm{wc}}^{{}^{\mathrm{15}}\mathrm{N}}$) and sedimentary denitrification (${\mathit{ϵ}}_{\mathrm{sed}}^{{}^{\mathrm{15}}\mathrm{N}}$) involves more considerations because we must account for the preference of ¹⁴NO₃ over ¹⁵NO₃. We begin with an ϵ of 5 ‰ for biological assimilation. This is equivalent to a ¹⁵NO₃:¹⁴NO₃ ratio of 0.995 when our atmospheric standard is equal to 1:1 using the following equation.

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle }\mathit{ϵ}=({\displaystyle \frac{{}^{\mathrm{15}}\mathrm{N}{/}^{\mathrm{14}}\mathrm{N}}{{}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{std}}{/}^{\mathrm{14}}{\mathrm{N}}_{\mathrm{std}}}}-\mathrm{1})\cdot \mathrm{1000}\end{array}$$

Note that a positive ϵ value returns a ¹⁵NO₃:¹⁴NO₃ ratio <1, while a negative δ¹⁵N in the previous example with N₂ fixation also returned a ¹⁵NO₃:¹⁴NO₃ ratio <1. This works because the reactions are in opposite directions. N₂ fixation adds NO₃, while assimilation removes NO₃. This means that 0.995 units of ¹⁵NO₃ are assimilated per unit of ¹⁴NO₃. As we have seen, a more useful way to quantify this is per unit of NO₃ assimilated into organic matter. Using Eq. (15), we find that ∼0.4987 units of ¹⁵NO₃ and ∼0.5013 units of ¹⁴NO₃ are assimilated per unit (1.0) of NO₃ when ϵ equals 5 ‰. Biological assimilation therefore leaves slightly more ¹⁵N in the unused NO₃ pool relative to ¹⁴N, which increases the δ¹⁵N of NO₃ while creating more ¹⁵N-deplete organic matter (δ¹⁵N_org).

However, we must also account for the effect that NO₃ availability has on fractionation. The preference of ¹⁴NO₃ over ¹⁵NO₃ strongly depends on the availability of NO₃, such that when NO₃ is abundant the preference for the lighter isotope will be strongest. This preference (fractionation) becomes weaker as NO₃ is depleted because cells will absorb any NO₃ that is available irrespective of its isotopic composition (Mariotti et al., 1981). Thus, as NO₃ is utilised, u, towards 100 % of its availability (u=1), the fractionation against ¹⁵NO₃ decreases to an ϵ of 0 ‰. This means that when u is equal to 1, no fractionation occurs and equal parts ¹⁵N and ¹⁴N (0.5:0.5 per unit NO₃) are assimilated. As we are interested in long timescales, we chose the accumulated product equations (Altabet and Francois, 2001) to approximate this process, where

$$\begin{array}{}\text{(17)}& {\displaystyle }& {\displaystyle }u=min(\mathrm{0.999},max(\mathrm{0.001},{\displaystyle \frac{{\mathrm{N}}_{\mathrm{org}}}{{\mathrm{NO}}_{\mathrm{3}}}}\left)\right)\text{(18)}& {\displaystyle }& {\displaystyle }{\mathit{ϵ}}_{\mathrm{u}}=\mathit{ϵ}\cdot {\displaystyle \frac{\mathrm{1}-u}{u}}\cdot \ln (\mathrm{1}-u).\end{array}$$

For numerical reasons, we limited the domain of u to (0.001,0.999) rather than (0,1), such that the utilisation-affected ϵ_u has a range of −4.997 to −0.035 ‰ for an ϵ of 5 ‰. ϵ_u is then converted into ratio units by dividing by 1000 and added to the ambient ¹⁵N:¹⁴N of NO₃ in the reactant pool to determine the ¹⁵N:¹⁴N of the product. In this case, it is the ¹⁵N:¹⁴N of newly created organic matter but could also be unused NO₃ effluxed from denitrifying cells in the case for denitrification.

$$\begin{array}{}\text{(19)}& {\displaystyle }{\displaystyle }{}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{org}}{:}^{\mathrm{14}}{\mathrm{N}}_{\mathrm{org}}={}^{\mathrm{15}}{\mathrm{NO}}_{\mathrm{3}}{:}^{\mathrm{14}}{\mathrm{NO}}_{\mathrm{3}}+{\mathit{ϵ}}_{\mathrm{u}}\end{array}$$

We then solve for how much ¹⁵NO₃ is assimilated into organic matter using Eq. (15) because we now know the change in NO₃ (ΔNO₃) and the ¹⁵N:¹⁴N of the product, which is ¹⁵N_org∕¹⁴N_org in our example of biological assimilation.

$$\begin{array}{}\text{(20)}& {\displaystyle }{\displaystyle }{\mathrm{\Delta }}^{\mathrm{15}}{\mathrm{NO}}_{\mathrm{3}}={\displaystyle \frac{{}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{org}}{/}^{\mathrm{14}}{\mathrm{N}}_{\mathrm{org}}\cdot \mathrm{\Delta }{\mathrm{NO}}_{\mathrm{3}}}{\mathrm{1}+{}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{org}}{/}^{\mathrm{14}}{\mathrm{N}}_{\mathrm{org}}}}\end{array}$$

Here, the change in ¹⁵NO₃ is equivalent to that assimilated into organic matter. Following assimilation into organic matter, the release of ¹⁵NO₃ through the water column during remineralisation occurs with no fractionation, such that the same δ¹⁵N signature is released throughout the water column.

We apply these calculations to each reaction in the nitrogen cycle that involves fractionation (assimilation, water column denitrification and sedimentary denitrification). They could be applied to any form of fractionation process with knowledge of ϵ, the isotopic ratio of the reactant, the amount of reactant that is used and the total amount of reactant available.

4 Model performance

CSIRO Mk3L-COAL adequately reproduces the large-scale thermohaline properties and circulation of the ocean under preindustrial conditions in numerous prior studies (Phipps et al., 2013; Matear and Lenton, 2014; Buchanan et al., 2016, 2018). Rather than reproduce these studies, we concentrate here on how the biogeochemical model performs relative to measurements of δ¹³C and δ¹⁵N in the water column (Eide et al., 2017, ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ data courtesy of the Sigman Lab at Princeton University) and in the sediments (Tesdal et al., 2013; Schmittner et al., 2017). We make these model–data comparisons alongside other isotope-enabled ocean general circulation models (Table 1).

All analyses of model performance were undertaken using the default parameterisation of the biogeochemical model, which is summarised in the tables of Sect. B in the Appendix. Each experiment was run towards steady state under preindustrial atmospheric conditions over many thousands of years. All results presented in this paper therefore reflect tracers that have achieved an equilibrium solution. We present annual averages of the equilibrium state in the following analysis.

Table 1 Models assessed against isotope data. The University of Victoria – Model of Ocean Biogeochemistry and Isotopes (UVic-MOBI) fields taken from Schmittner and Somes (2016). Pelagic Interactions Scheme for Carbon and Ecosystem Studies (PISCES) fields provided by Laurent Bopp. LOch–Vecode-Ecbilt-CLio-agIsm Model (LOVECLIM) fields taken from Menviel et al. (2017a). The isotope-enabled Community Earth System Model (iCESM) fields for δ¹³C (low resolution) provided by Alexandra Jahn (Jahn et al., 2015) and those for δ¹⁵N (high resolution) provided by Simon Yang (Yang and Gruber, 2016). PISCES and CESM model resolutions have a range of longitude–latitude spacings to reflect regions of finer resolution, including the Equator and polar regions.

4.1 δ¹³C of dissolved inorganic carbon (δ¹³C_DIC)

The recent reconstruction of preindustrial δ¹³C_DIC by Eide et al. (2017) provides a large dataset for comparison. We chose this dataset over the compilation of point location water column data of Schmittner et al. (2017) because it offers a gridded product where short-term and small-scale variability are smoothed, making for more appropriate comparison with model output.

Predicted values of δ¹³C_DIC from CSIRO Mk3L-COAL broadly replicated the preindustrial distribution. The predicted global mean of 0.41 ‰ reflected that of the reconstructed mean of 0.42 ‰ (Table 2). Spatial agreement was acceptable with a global correlation of 0.80 (G marker in Fig. 2). Regionally, the Southern Ocean performed well with the lowest root mean square (rms) error of 0.42 ‰, while a greater degree of disagreement in the values of δ¹³C_DIC existed in the middle and lower latitudes of each major basin, particularly in the Atlantic where model–data agreement (correlation, rms error and normalised standard deviation) was poorest. Subsurface δ¹³C_DIC was too low in the tropics of the major basins by ∼0.2 ‰ and too high in the North Pacific and North Atlantic by 0.4 ‰ to 0.6 ‰ (Fig. 3).

Eide et al. (2017)

Table 2 Comparison of global and region mean δ¹³C_DIC between observations (Eide et al., 2017) and model simulations. Means are annual averages and do not include the Arctic or the upper 200 m of the water column. All data were regridded onto the CSIRO Mk3L-COAL grid space.

Figure 2 Global and regional fits between data (Eide et al., 2017) and simulated δ¹³C of dissolved inorganic carbon displayed as Taylor diagrams (Taylor, 2001). Shading of the markers represents normalised bias. G indicates global; S indicates Southern Ocean (90–40^∘ S); A indicates Atlantic (40^∘ S–70^∘ N); P indicates Pacific (40^∘ S–70^∘ N); I indicates Indian Ocean (40^∘ S–70^∘ N). Measures of fit do not include the Arctic or the upper 200 m of the water column. All data were regridded onto the CSIRO Mk3L-COAL grid space before comparison.

Figure 3 Zonal mean observed (a, b, c) and modelled (d, e, f) δ¹³C of DIC produced by CSIRO Mk3L-COAL for each major basin. The red dashed line marks the upper 175 m and is used for comparison between observed and modelled distributions. Replicate figures for the other models are available in the Supplement.

These inconsistencies were likely related to physical and biological limitations of CSIRO Mk3L-COAL. δ¹³C_DIC in subsurface tropical waters was too low because restricted horizontal mixing and high carbon export drove very negative δ¹³C_DIC values. The very negative δ¹³C values were associated with very large oxygen minimum zones and were thus a product of poorly represented, fine-scale equatorial dynamics. Coarse-resolution OGCMs are known to have weak equatorial undercurrents that lead to oxygen minimum zones that are too large (Matear and Holloway, 1995; Oschlies, 2000) and CSIRO Mk3L-COAL is no exception. Alternatively, the large oxygen minimum zones could be due to our conservative treatment of organic matter remineralisation (Sect. A in the Appendix), where remineralisation is prevented when O₂ and NO₃ are unavailable. Organic matter therefore falls deeper into the interior through oxygen-deficient zones, leading to their vertical expansion. Almost certainly, however, it was the poorly represented dynamics within the Pacific basin that were responsible for high δ¹³C_DIC in the subsurface North Pacific, which contains low O₂ and low δ¹³C_DIC water due to northward transport from the tropics.

Another inconsistency was a positive bias in the upper 200–500 m, with values exceeding 2 ‰ in many areas of the lower latitudes. However, values as high as 2 ‰ have been measured in the upper 500 m of the Indo-Pacific (Schmittner et al., 2017). Given the difficulties associated with accounting for the Suess effect (invasion of isotopically light fossil fuel CO₂), it is possible that the upper ocean values of Eide et al. (2017) underestimate the preindustrial δ¹³C_DIC surface field. It is also equally possible that a fixed biological fractionation (${\mathit{ϵ}}_{\mathrm{bio}}^{{}^{\mathrm{13}}\mathrm{C}}$) of 21 ‰ may have driven unrealistic enrichment in the simulated field. High growth rates are thought to lower the strength of fractionation during carbon fixation (Laws et al., 1995). To explore the possibility of model–data mismatch caused by our choice to fix ${\mathit{ϵ}}_{\mathrm{bio}}^{{}^{\mathrm{13}}\mathrm{C}}$ at 21 ‰, we implemented biological fractionation that is dependent on phytoplankton growth rate and aqueous CO₂ concentration (Eq. 6). We found the implementation of a variable ${\mathit{ϵ}}_{\mathrm{bio}}^{{}^{\mathrm{13}}\mathrm{C}}$ reduced high values in the upper part of the low-latitude ocean but that this reduction was small (Fig. 4). The overwhelming effect was an increase in δ¹³C_DIC throughout the interior, itself caused by weaker fractionation in the tropical ocean. Global mean δ¹³C_DIC subsequently increased by 0.25 ‰. Meanwhile, model skill was unaffected (see CSIRO Mk3L-COAL (vary-${\mathit{ϵ}}_{\mathrm{bio}}^{{}^{\mathrm{13}}\mathrm{C}}$) in Fig. 2). Neither fixed nor variable biological fractionation could reproduce the low upper ocean values of the data.

Figure 4 The introduction of variable carbon fractionation by phytoplankton (a) and the consequent change in δ¹³C_DIC represented as a zonal mean (b) relative to a case where ${\mathit{ϵ}}_{\mathrm{bio}}^{{}^{\mathrm{13}}\mathrm{C}}$ is fixed at 21 ‰.

It is helpful to place our predicted δ¹³C_DIC alongside those of other global ocean models (Fig. 2; Table 2), both for skill assessment and to further understand the cause of the positive bias in the upper ocean. We take annually averaged, preindustrial δ¹³C_DIC distributions from the UVic-MOBI, PISCES, LOVECLIM and iCESM-low biogeochemical models, most of which have been used in significant palaeoceanographic modelling studies (Menviel et al., 2017a; Tagliabue et al., 2009; Schmittner and Somes, 2016). Predicted δ¹³C_DIC performs adequately in CSIRO Mk3L-COAL relative to these state-of-the-art models. LOVECLIM showed good fit in terms of global and regional means (Table 2) but had lower correlations (Fig. 2), suggesting that its values were accurate but its distribution biased. UVic-MOBI had high correlations, but it consistently overestimated the preindustrial field by ∼0.2 ‰. Interestingly, the bias of UVic-MOBI, which treats biological fractionation as a function of growth rate and aqueous CO₂, is similar to CSIRO Mk3L-COAL when this form of fractionation was activated (vary-${\mathit{ϵ}}_{\mathrm{bio}}^{{}^{\mathrm{13}}\mathrm{C}}$). PISCES and iCESM-low were the best performing models, equally demonstrating high correlations, low biases, accurate regional and global means and the lowest rms errors. This is perhaps not surprising considering the significantly finer vertical resolutions of these OGCMs and their more complex horizontal grid structure that enables an improved representation of ocean dynamics (Table 1). However, all models performed most poorly in the Atlantic Ocean, with poor correlations, high variability and greater biases.

Returning to the consistent positive bias in the upper ocean, most models (except iCESM-low) predicted upper ocean δ¹³C_DIC ≥ 2.0 ‰ (Figs. S1, S2, S3 and S4 in the Supplement) similar to CSIRO Mk3L-COAL. As each model has a unique representation of the marine ecosystem and consequently a unique treatment of biological fractionation, the common prediction of high upper ocean δ¹³C once more suggests that the upper ocean values between 200 and 500 m of Eide et al. (2017) may be too low. The underestimation of δ¹³C_DIC may be due to a neglect of biology introducing anthropogenic, isotopically depleted carbon to surface and subsurface layers via remineralisation (the biological Suess effect). This would in turn suggest that a higher global mean of 0.73 ‰ generated from a global compilation of foraminiferal δ¹³C (Schmittner et al., 2017) is perhaps a more accurate representation of preindustrial δ¹³C values.

Overall, CSIRO Mk3L-COAL performed acceptably in terms of its mean values and correlations but had consistently greater rms errors in major basins outside of the Southern Ocean. This indicates that CSIRO Mk3L-COAL exaggerated regional minima and maxima as discussed. Despite the regional biases of CSIRO Mk3L-COAL, the comparison demonstrates that all models have strengths and weaknesses. Given its low resolution and computational efficiency, CSIRO Mk3L-COAL performs adequately among other biogeochemical models in its simulation of δ¹³C_DIC.

4.2 δ¹³C of Cibicides foraminifera (δ¹³C_Cib)

We extended our assessment of modelled δ¹³C_DIC by comparing it to a compilation of benthic δ¹³C measured within the calcite of foraminifera from the genus Cibicides (Schmittner et al., 2017), a genus on which much of the palaeoceanographic δ¹³C records are based. For this comparison, we adjusted our predicted δ¹³C_DIC to predicted δ¹³C_Cib using the linear dependence on carbonate ion concentration and depth suggested by Schmittner et al. (2017):

$$\begin{array}{ll}{\displaystyle }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{Cib}}=& {\displaystyle }\mathrm{0.45}+{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}-\mathrm{2.2}\times {\mathrm{10}}^{-\mathrm{3}}\\ \text{(21)}& {\displaystyle }& {\displaystyle }\cdot {\mathrm{CO}}_{\mathrm{3}}-\mathrm{6.6}\times {\mathrm{10}}^{-\mathrm{5}}\cdot z.\end{array}$$

This adjustment accounts for slight fractionation during incorporation of DIC into foraminiferal calcite and is found to be partly explained by the concentration of ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ ions and pressure. A one-to-one comparison between δ¹³C_DIC and δ¹³C_Cib hence introduces some degree of error since this fractionation is not accounted for. Because we are interested in applying simulated δ¹³C_DIC to a palaeoceanographic context, we must first be able to convert our simulated δ¹³C_DIC to δ¹³C_Cib in an effort to make better comparisons, particularly as the distribution of ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ is subject to change. By adjusting our three-dimensional δ¹³C_DIC output using Eq. (21), we attain predicted δ¹³C_Cib (see inset titled “calibration” in Fig. 5). For good measure, we also computed measures of statistical fit for a traditional one-to-one comparison between δ¹³C_DIC and δ¹³C_Cib to assess the benefit of the calibration.

Figure 5 Measured versus modelled δ¹³C_Cib (N=690) of CSIRO Mk3L-COAL coloured by latitude. Red shading about the 1:1 line is an estimate of the variability implicit in the relationship between δ¹³C_Cib and δ¹³C_DIC of Schmittner et al. (2017). The inset on the bottom right shows the effect of the calibration of Eq. (21).

Measured δ¹³C_Cib data from Schmittner et al. (2017) were binned into model grid boxes and averaged for the comparison. Those measurements that fell within the OGCM's land mask were excluded. Transfer and averaging onto the coarse-resolution OGCM grid reduced the number of points for comparison from 1763 to 690, lowered the mean of measured δ¹³C_Cib from 0.76 ‰ to 0.52 ‰ and reduced the absolute range from $-\mathrm{0.9}\to \mathrm{2.1}$ to $-\mathrm{0.7}\to \mathrm{2.1}$.

Adjusted δ¹³C_Cib using Eq. (21) showed good fit to measured δ¹³C_Cib given the sparsity of data, with a global correlation of 0.64, a mean of 0.57 ‰ and an rms error of 0.63 ‰ (Table 3). If a one-to-one relationship between δ¹³C_DIC and δ¹³C_Cib was used, the global correlation was not affected and only slightly worse skill was detected in mean, rms error and standard deviation. Accounting for the regional influence of carbonate ion concentration and depth was therefore beneficial, likely because very negative and positive values were slightly adjusted towards the mid-range (inset in Fig. 5), but this was not necessary for an adequate comparison. This conclusion was also reached by Schmittner et al. (2017). Likewise, implementing variable fractionation by phytoplankton (${\mathit{ϵ}}_{\mathrm{bio}}^{{}^{\mathrm{13}}\mathrm{C}}$) had little effect except to increase values and slightly improve measures of skill (Table 3). Of the 690 data points used in the comparison, 419 fell within the error around what could be considered a good fit (shaded red area in Fig. 5). The error was taken as 0.29 ‰ and represents the standard deviation associated with the relationship between δ¹³C_DIC and δ¹³C_Cib measurements (Schmittner et al., 2017).

Table 3 Statistical comparison of core-top δ¹³C_Cib with predicted values produced by the CSIRO Mk3L-COAL ocean model.

Some notable over and underestimation occurred in the adjusted δ¹³C_Cib output that more or less mirrored those inconsistencies previously discussed for δ¹³C_DIC. Values as low as −1.9 ‰, well below measured δ¹³C_Cib minima of −0.7 ‰, existed in the equatorial subsurface Pacific and Indian oceans (i.e. where the oxygen minimum zones existed). This can be seen in Fig. 5, where some values in the equatorial band are well below the shaded region of good fit. Meanwhile, very high values of δ¹³C_Cib were predicted in Arctic surface waters. The exaggeration of these local minima and maxima reflects those found in the modelled δ¹³C_DIC distribution. Despite these local inconsistencies, CSIRO Mk3L-COAL shows good potential for direct comparisons to palaeoceanographic data sets of foraminiferal δ¹³C with or without calibration.

4.3 δ¹⁵N of nitrate (${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$)

We produced univariate measures of fit by comparing measurements of ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ with equivalent values from CSIRO Mk3L-COAL at the nearest point (Fig. 6; Table 4). Measured ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ values were collected over a 30-year period using a variety of collection and measurement methods with a distinct bias towards the Atlantic Ocean. To try and remove some temporal and spatial bias, we binned and averaged measurements into equivalent model grids.

Table 4 Comparison of global and region mean ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ between observations and model simulations. Model means are annual averages. All data were regridded onto the CSIRO Mk3L-COAL grid space. The δ¹⁵N data (5330 measurements courtesy of the Sigman Lab, Princeton University) were binned into corresponding grid boxes and averaged for direct comparison, which reduced the data to 2532 points. More than one data point of δ¹⁵N may therefore contribute to each simulated value.

Figure 6 Global and regional fits between observations and simulated ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ displayed as Taylor diagrams (Taylor, 2001). Shading of the markers represent normalised bias. G indicates global; S indicates Southern Ocean (90–40^∘ S); A indicates Atlantic (40^∘ S–70^∘ N); P indicates Pacific (40^∘ S–70^∘ N); I indicates Indian Ocean (40^∘ S–70^∘ N). The δ¹⁵N data (5330 measurements courtesy of the Sigman Lab, Princeton University) were binned into corresponding grid boxes and averaged for direct comparison, which reduced the data to 2532 points. More than one data point of δ¹⁵N may therefore contribute to each simulated value.

CSIRO Mk3L-COAL adequately reproduced the global patterns of ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$. We found excellent agreement in the volume-weighted means of ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ (Table 4). Tight agreement in the means was a consequence of reproducing similar values where the majority of observed data existed. Most ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ measurements have been taken from the upper 1000 m in the North Atlantic where values cluster at just under 5 ‰ (see Fig. 7a, c, e). Closer inspection of the Atlantic using depth and zonally averaged sections (Figs. 8 and 9) revealed that the model adequately reproduced the low δ¹⁵N signature of ∼4 ‰ caused by N₂ fixation occurring in the tropical Atlantic (Marconi et al., 2017). A basin-wide rate of Atlantic N₂ fixation equal to ∼33 Tg N yr⁻¹ lowered Atlantic values below 5 ‰ and was fundamental for reproducing the observations. Outside the Atlantic, where data are more sparse, the model successfully reproduced the meridional gradients across the Antarctic, Subantarctic and subtropical zones, the subsurface ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ maxima in the tropics of all major basins and the tongues of high and low values in surface waters of the Pacific consistent with changes in nitrate utilisation (Figs. 8 and 9).

Figure 7 Observed (a, c, e) and modelled (b, d, f) δ¹⁵N of NO₃ data (N=5004) plotted against depth (a, b), latitude (c, d) and longitude (e, f). Colour shading represents the density of data, such that the darker a mass of data points is, the more data are represented there.

Figure 8 Depth-averaged sections of modelled (colour contours) and observed (overlaid markers) ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$.

Figure 9 Zonally averaged sections of modelled (colour contours) and observed (overlaid markers) ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$. The global zonal average encompasses all basins.

Some important regional inconsistencies between the simulated and measured values did exist (refer to Figs. 8 and 9) and degraded the correlation. Much like the high values of δ¹³C_DIC that were transported too deeply into the North Atlantic interior, a low ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ signature was transported too far into the deep North Atlantic. CSIRO Mk3L-COAL therefore underestimated deep ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ before mixing through to the South Atlantic restored values towards the measurements. Subsurface values in the North Pacific were also underestimated, which can be attributed to the inability of the coarse-resolution OGCM to transport low O₂, high ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ water northwards from the eastern tropical Pacific. Simulated values in the Indian Ocean, specifically near the Arabian Sea, significantly underestimated the data because the suboxic zone was misrepresented in the Bay of Bengal. Misrepresentation of the northern Indian Ocean was responsible for very poor model–data fit in the Indian Ocean (Fig. 6). Meanwhile, the deep (>1500 m) eastern tropical Pacific tended to overestimate the data, due to a large, deep, unimodal suboxic zone. These physically driven inconsistencies in the oxygen field are common to other coarse-resolution models (Oschlies et al., 2008; Schmittner et al., 2008) and, like the δ¹³C distribution, were the main cause of the misfit between simulated and observed ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$. The correlations reflected these regional under- and overestimations, particularly in the Indian Ocean.

Finally, we placed CSIRO Mk3L-COAL in the context of other isotope-enabled global models: UVic-MOBI, PISCES and iCESM-high (Table 1). This comparison demonstrated that the modelled distribution of ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ was adequately placed among the current generation of models. The global and regional means were more accurately reproduced by CSIRO Mk3L-COAL than for UVic-MOBI, PISCES and iCESM-high (Table 4; also see shading in Fig. 6). Atlantic ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ was best reproduced by CSIRO Mk3L-COAL. Meanwhile, correlations tended to be slightly lower for CSIRO Mk3L-COAL than UVic-MOBI and iCESM-high, and consistently lower than PISCES (Fig. 6). UVic-MOBI underestimated the data but produced high correlations in the Southern Ocean and globally. Regionally, PISCES was best correlated to the measurements of ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ of the three models, although it had a consistent positive bias. iCESM-high was acceptably correlated with the data in the global sense but was highest in rms errors, particularly in the Pacific. CSIRO Mk3L-COAL therefore showed an acceptable measure of fit to the noisy and sparse ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ data and reproduced most regional patterns, albeit with misrepresentation in the Indian Ocean and some exaggerations of local minima/maxima as discussed. Future model–data comparisons with CSIRO Mk3L-COAL should therefore take these limitations into account. Overall, however, we find that CSIRO Mk3L-COAL broadly reproduced the ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ data. Annual rates of N₂ fixation, water column denitrification and sedimentary denitrification at roughly 122, 52 and 78 Tg N yr⁻¹, respectively, produced this agreement.

An important caveat to the ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ routines of CSIRO Mk3L-COAL should be noted. CSIRO Mk3L-COAL underwent significant tuning of water column and sedimentary denitrification parameterisations in order to reproduce known values of ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ during development. One important parameter is the lower threshold of NO₃ concentration at which point water column denitrification is shut off (Sect. A2.3). In CSIRO Mk3L-COAL, this is set at 30 mmol m⁻³, which is an arbitrary limit that was implemented to prevent water column denitrification from reducing NO₃ to zero in the large suboxic zones. Hence, a caveat of the current model is an inability for water column and sedimentary denitrification to realistically adjust as suboxia changes. However, the parameterisation does allow for targeted experiments where the ratio of water column to sedimentary denitrification can be controlled if, for instance, it is unclear how water column and sedimentary denitrification respond to certain conditions. This is currently the case during the Last Glacial Maximum, where expansive suboxic zones in the Pacific (Hoogakker et al., 2018) were counter-intuitively associated with reduced rates of water column denitrification (Ganeshram et al., 1995). We have, in this version, chosen to keep this parameterisation and note that future developments will focus on dynamic responses to variations in suboxia.

4.4 δ¹⁵N of organic matter (δ¹⁵N_org)

CSIRO Mk3L-COAL tracks the δ¹⁵N signature of organic matter (δ¹⁵N_org) that is deposited in the sediments. We compared the simulated δ¹⁵N_org to the core-top compilation of Tesdal et al. (2013) with 2176 records of δ¹⁵N_org. These records were binned and averaged onto the CSIRO Mk3L-COAL ocean grid, such that the 2176 records became 592. When comparing sediment core-top measurements of δ¹⁵N to that of the model, it is necessary to consider how δ¹⁵N_org is altered by early burial. As records in the compilation of Tesdal et al. (2013) are from bulk nitrogen, we can assume that the “diagenetic offset” as described by Robinson et al. (2012) is active. The diagenetic offset involves an increase in the δ¹⁵N of sedimentary nitrogen of between 0.5 ‰ and 4.1 ‰ relative to that of sinking particulate organic matter and appears to be related to pressure (Robinson et al., 2012), although the reasoning behind this relationship remains to be defined.

In light of the diagenetic offset, we make three comparisons with the compilation of Tesdal et al. (2013). A raw comparison is made, alongside an attempt to account for the diagenetic offset using two depth-dependent corrections (Table 5 and Fig. 10):

$$\begin{array}{}\text{(22)}& {\displaystyle }& {\displaystyle }{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{org}}^{\mathrm{cor}:\mathrm{1}}=\left\{\begin{array}{ll}{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{org}},& \text{if}z\left(\mathrm{km}\right)<\mathrm{1}\mathrm{km}\\ {\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{org}}+(\mathrm{1}\cdot z(\mathrm{km})+\mathrm{1}),& \text{if}z\left(\mathrm{km}\right)\ge \mathrm{1}\mathrm{km}\end{array}\right.\text{(23)}& {\displaystyle }& {\displaystyle }{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{org}}^{\mathrm{cor}:\mathrm{2}}={\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{org}}+\mathrm{0.9}\cdot z\left(\mathrm{km}\right).\end{array}$$

The first correction (${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{org}}^{\mathrm{cor}:\mathrm{1}}$) is taken from Robinson et al. (2012), while the second (${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{org}}^{\mathrm{cor}:\mathrm{2}}$) originates from how Schmittner and Somes (2016) treated sedimentary nitrogen isotope data in their study of the Last Glacial Maximum. Both are based on the observation that the diagenetic offset increases with pressure, in this case represented by depth (z) in kilometres (km).

Table 5 Statistical comparison of core-top δ¹⁵N_org with predicted values of the CSIRO Mk3L-COAL ocean model. The corrected vales (${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{org}}^{\mathrm{cor}:\mathrm{1}}$ and ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{org}}^{\mathrm{cor}:\mathrm{2}}$) account for alteration during early diagenesis following burial.

Following binning and averaging onto the model grid, the raw comparison immediately showed a consistent underestimation of the core-top data, with a predicted mean of 2.7 ‰ well below the observed mean of 4.7 ‰. Our correlation was 0.27, which indicates a limited ability to replicate regional patterns. This underestimation and low correlation is easily seen when predicted values are compared directly to the core-top data in Fig. 10. Like the nitrogen isotope model of Somes et al. (2010), we find that the offset between simulated and observed core-top bulk δ¹⁵N_org is roughly equivalent to the observed average diagenetic offset of $\sim \mathrm{2.3}\pm \mathrm{1.8}$ ‰. This indicates that diagenetic alteration of δ¹⁵N_org is active during early burial in the core-top data.

Including a diagenetic offset therefore improved agreement between our predicted δ¹⁵N_org and the core-top data considerably (Table 5 and Fig. 10). Both corrections accounted for the enrichment of δ¹⁵N in deeper regions and the minor diagenetic alteration in areas of high sedimentation that typically occurs in shallower sediments. The average δ¹⁵N_org increased to 4.5 ‰ for ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{org}}^{\mathrm{cor}:\mathrm{1}}$ and 5.2 ‰ for ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{org}}^{\mathrm{cor}:\mathrm{2}}$. Correlations increased from 0.27 to 0.47 and 0.53, respectively. The improvement was clearly observed in the Southern Ocean, where both the magnitude and spatial pattern of δ¹⁵N_org were well replicated by the model. Changes in the Southern Ocean over glacial–interglacial cycles reflect shifts in the global marine nitrogen cycle and nutrient utilisation (Martinez-Garcia et al., 2014; Studer et al., 2018), and the ability of CSIRO Mk3L-COAL to account for these patterns in the core-top data is encouraging for future study. We suggest that future palaeoceanographic model–data comparisons of δ¹⁵N_org use the depth correction of Schmittner and Somes (2016) as it provided the best correlations and reproduced Southern Ocean δ¹⁵N_org at 0.5 ‰ greater than the global mean (see Table 5).

Figure 10 Direct comparison of observed versus modelled δ¹⁵N_org incident on the sediments. Panels (a), (c) and (e) show spatial distribution of simulated δ¹⁵N_org overlain by core-top data from the compilation of Tesdal et al. (2013). Panels (b), (d) and (f) compare all core-top data against simulated δ¹⁵N_org. Panels (a) and (b) depict raw output of the model, while panels (c)–(f) depict the predicted values of the model following two depth-dependent offsets (Eqs. 22 and 23) that account for diagenetic alteration.

5 Ecosystem effects

As a first test of the isotope-enabled ocean model, we undertook simple ecosystem experiments to assess the effect on δ¹³C and δ¹⁵N. For reference, the assessment of model performance described above used model output with variable stoichiometry activated, a fixed 8 % rain ratio of CaCO₃ to organic carbon and a strong iron limitation of N₂ fixers that enforced a low degree of spatial coupling between N₂ fixers and denitrification zones. A summary of the biogeochemical effects of the different experiments is provided in Table 6.

Table 6 Summary of the biogeochemical effects of the different treatments of the ecosystem in CSIRO Mk3L-COAL. C_org is the total organic carbon exported from the euphotic zone composed of both general and diazotrophic phytoplankton groups (${\mathrm{C}}_{\mathrm{org}}^{\mathrm{G}}$ + ${\mathrm{C}}_{\mathrm{org}}^{\mathrm{D}}$; see Sect. A1 in the Appendix), while ${\mathrm{C}}_{{\mathrm{CaCO}}_{\mathrm{3}}}$ is the total export of CaCO₃ out of the euphotic zone. The sum of C_org and ${\mathrm{C}}_{{\mathrm{CaCO}}_{\mathrm{3}}}$ are equal to the global rate of carbon export referred to in the text. Sed : WC refers to the sedimentary to water column denitrification ratio. Note that the global mean δ¹³C_DIC is higher than reported in Table 2 because it includes the upper 200 m and the Arctic.

5.1 Variable versus Redfieldian stoichiometry

Enabling variable stoichiometry (see Sect. A3 in the Appendix) of the general phytoplankton group (${\mathrm{P}}_{\mathrm{org}}^{\mathrm{G}}$) over a Redfieldian ratio ($\mathrm{C}:\mathrm{N}:\mathrm{P}:{\mathrm{O}}_{\mathrm{2}}^{\mathrm{rem}}:{\mathrm{NO}}_{\mathrm{3}}^{\mathrm{rem}}$ $=\mathrm{106}:\mathrm{16}:\mathrm{1}:-\mathrm{138}:-\mathrm{94.4}$) altered the rate and distribution of organic matter export. Organic matter had more carbon and nitrogen per unit phosphorus in regions with low PO₄, such as the Atlantic Ocean (Fig. 11a), which elevated O₂ and NO₃ demand during oxic and suboxic remineralisation (denitrification), respectively. Lower ratios were produced in eutrophic regions such as the subarctic Pacific, Southern Ocean and tropical zones of upwelling. Overall, global mean C:P increased from the Redfieldian 106:1 to 117:1 and caused an increase in carbon export from 7.6 to 8.0 Pg C yr⁻¹. Approximately 0.1 Pg C yr⁻¹, or 25 % of the increase, was attributed purely to organic carbon export from N₂ fixation, which increased from 107 to 122 Tg N yr⁻¹ as higher N:P ratios in the tropics broadened their competitive niche. The total contribution of N₂ fixation to the increase in carbon export was likely greater than 25 %, as NO₃ also became more available to NO₃-limited ecosystems in the lower latitudes (Moore et al., 2013). The increase in carbon export under variable stoichiometry as compared to a Redfieldian ocean was therefore felt largely in the lower latitudes between 40^∘ S and 40^∘ N (Fig. 11b). Export production decreased poleward of 40^∘, particularly in the Southern Ocean, because C:P ratios were lower than the 106:1 Redfield ratio (Fig. 11a).

Distributions of both isotopes were affected by the change in carbon export and the marine nitrogen cycle. Global mean δ¹³C_DIC increased from 0.47 ‰ to 0.51 ‰ and ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ increased from 5.1 ‰ to 5.6 ‰. These are not great changes on the global scale and they had little influence on model–data measures of fit. However, the spatial distribution of these isotopes was significantly altered. Intermediate waters leaving the Southern Ocean were depleted in δ¹³C_DIC by up to 0.1 ‰ and ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ by up to 1 ‰, while the deep ocean, particularly the Pacific, was enriched in both isotopes to a similar degree (Fig. 12). Depletion of both isotopes in waters subducted between 40 and 60^∘ S reflected the local loss in export production as a result of lower C:P and N:P ratios, such that biological fractionation was unable to enrich DIC and NO₃ in the heavier isotope to the same degree as surface waters travelled north. Enrichment of δ¹³C in the deep ocean was the result of reduced carbon export in the Antarctic zone due to low C:P ratios, while enrichment of δ¹⁵N in the deep ocean was the result of increased tropical production that increased water column denitrification (${\mathit{ϵ}}_{\mathrm{wc}}^{{}^{\mathrm{15}}\mathrm{N}}$ =20 ‰). Lower C:P and N:P ratios in both the Antarctic and Subantarctic zones therefore elicited divergent isotope effects in deep and intermediate waters leaving the Southern Ocean.

Meanwhile, each isotope showed a different response in the suboxic zones of the tropics where variable stoichiometry increased the volume of suboxia (O₂<10 mmol m⁻³) by 0.5 %. The increase in water column denitrification caused by the expansion of suboxia increased ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$, while the local increase in carbon export that drove the increase in water column denitrification reduced δ¹³C_DIC in the same waters (Fig. 12). Overall, the increase in low-latitude carbon export caused an expansion of water column suboxia and elicited diverging behaviours in the isotopes, whereby ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ increased and δ¹³C_DIC decreased.

Figure 11 Simulated difference in the C:P ratio of exported organic matter due to variable stoichiometry as compared to Redfield stoichiometry (a) and the resulting change in carbon export out of the euphotic zone (b).

Figure 12 Differences in δ¹³C_DIC (a) and ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ (b) as a result of variable stoichiometry as compared to Redfield stoichiometry. Values are zonal means.

5.2 Calcifier dependence on calcite saturation state

The rate of calcification of planktonic foraminifera and coccolithophores is dependent on the calcite saturation state (Zondervan et al., 2001). In previous experiments, the production of CaCO₃ was fixed at a rate of 8 % per unit of organic carbon produced in accordance with the modelling study of Yamanaka and Tajika (1996) and produced 0.54 Pg CaCO₃ yr⁻¹. Now we investigate how spatial variations in the CaCO₃:C_org ratio ($R_{{\mathrm{CaCO}}_{\mathrm{3}}}$ in Eq. A17) affected δ¹³C_DIC and δ¹³C_Cib (see Sect. A1.3 in the Appendix). We applied three different values of η to Eq. (A18) to alter the quantity of CaCO₃ produced per unit of organic carbon (${\mathrm{C}}_{\mathrm{org}}^{\mathrm{G}}$) given the calcite saturation state (Ω_ca). The η coefficients were 0.53, 0.81 and 1.09. These numbers are equivalent to those in the experiments of Zhang and Cao (2016).

Figure 13 Global distribution of CaCO₃ export as a percentage of organic carbon (C_org) export (a) and the change in the CaCO₃ production field as a result of making CaCO₃ production dependent on calcite saturation state (η=1.09) compared to when it was a fixed 8 % of C_org (b). Areas where export production does not occur due to severely nutrient limited conditions are masked out.

Figure 14 Changes in the distribution of carbon isotopes (δ¹³C_DIC and δ¹³C_Cib; a) and carbon chemistry (dissolved inorganic carbon and alkalinity; b) as a result of increasing CaCO₃ production in surface waters between 40^∘ S and 40^∘ N.

Mean $R_{{\mathrm{CaCO}}_{\mathrm{3}}}$ was 4.5 %, 6.6 % and 9.5 %, and annual CaCO₃ production was 0.32, 0.47 and 0.68 Pg CaCO₃ yr⁻¹ in the three experiments. Although different in total CaCO₃ production, the three experiments shared the same spatial patterns. Low-latitude waters were high in $R_{{\mathrm{CaCO}}_{\mathrm{3}}}$, particularly the oligotrophic subtropical gyres, while high latitudes were low, particularly the Antarctic zone where mixing of old waters into the surface depressed the calcite saturation state (Fig. 13). These regional patterns in $R_{{\mathrm{CaCO}}_{\mathrm{3}}}$ therefore had the largest effect in areas of high export production. Productive, high-latitude areas like the Southern Ocean, subpolar Pacific and North Atlantic waters all produced less CaCO₃ when compared to an enforced 8 % rain ratio, while CaCO₃ production between 40^∘ S and 40^∘ N relative to a fixed $R_{{\mathrm{CaCO}}_{\mathrm{3}}}$ of 8 % was dependent on η. The highest η coefficient of 1.09 achieved greater export of CaCO₃ in the mid- to lower-latitude regions of high export production (Fig. 13). The consequence of increasing CaCO₃ production in the middle–lower latitudes was a loss of upper ocean alkalinity, subsequent outgassing of CO₂ and losses in the DIC inventory. Losses in global DIC were 95 and 130 Pg C as $R_{{\mathrm{CaCO}}_{\mathrm{3}}}$ increased from $\mathrm{4.6}\to \mathrm{6.6}\to \mathrm{9.5}$ % (Table 6), equivalent to one-fifth of the glacial increase in oceanic carbon (Ciais et al., 2011).

Despite the significant changes associated with the implementation of Ω_ca-dependent CaCO₃ production, effects were negligible on both δ¹³C_DIC and δ¹³C_Cib. Global mean δ¹³C_DIC was 0.51 ‰, when $R_{{\mathrm{CaCO}}_{\mathrm{3}}}$ was fixed at 8 %, and this changed to 0.52 ‰, 0.50 ‰ and 0.48 ‰ under η coefficients of 0.53, 0.81 and 1.09 (Table 6). Likewise, global mean δ¹³C_Cib was 0.59 ‰, when $R_{{\mathrm{CaCO}}_{\mathrm{3}}}$ was fixed at 8 %, and this changed to 0.60 ‰, 0.58 ‰ and 0.55 ‰. Minimal change in δ¹³C_Cib indicated minimal change in the ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ concentration (see Eq. 21), which varied by ≤2 mmol m⁻³ between experiments. Visual inspection of the change in δ¹³C_DIC and δ¹³C_Cib distributions showed an enrichment of these isotopes in the upper ocean north of 40^∘ S. Subsequent increases in η, which increased low-latitude CaCO₃ production, magnified the enrichment. Enrichment of δ¹³C_DIC and δ¹³C_Cib was caused by outgassing of CO₂ as surface alkalinity decreased in response to greater CaCO₃ production (Fig. 14). The change, however, was at most 0.1 ‰, which lies well within 1 standard deviation of variability known in the proxy data (Schmittner et al., 2017). We therefore find little scope for recognising even large variations in global CaCO₃ production (0.32 to 0.68 Pg CaCO₃ yr⁻¹) in the signature of carbon isotopes despite considerable effects on the oceanic inventory of DIC.

However, we stress that version 1.0 of CSIRO Mk3L-COAL does not include CaCO₃ burial or dissolution from the sediments according the calcite saturation state of overlying water (Boudreau, 2013). To neglect ocean–sediment CaCO₃ cycling is to neglect of an important aspect of the global carbon cycle active on millennial timescales (Sigman et al., 2010). Changes in CaCO₃ burial and dissolution could have a non-negligible effect on δ¹³C through altering whole ocean alkalinity and thereby air–sea gas exchange of CO₂, which would in turn affect surface δ¹³C as we have seen. While we do not address these effects here, we aim to do so in upcoming versions of the model equipped with carbon compensation dynamics.

Figure 15 Changes in the distribution of marine N₂ fixation caused by altering how limiting iron is to the growth of N₂ fixers via the coefficient $K^{{\mathrm{D}}_{\mathrm{Fe}}}$ in Eq. (A12). Iron limitation is sequentially relaxed from top to bottom.

5.3 Strength of coupling between N₂ fixation and denitrification

The degree to which N₂ fixers are spatially coupled to the tropical denitrification zones is controlled by altering the degree to which N₂ fixers are limited by iron ($K^{{\mathrm{D}}_{\mathrm{Fe}}}$) in Eq. (A12) (see Sect. A1.2 in the Appendix). Decreasing $K^{{\mathrm{D}}_{\mathrm{Fe}}}$ ensures that N₂ fixation becomes less dependent on iron supply and as such is released from regions of high aeolian deposition, such as the North Atlantic, to inhabit areas of low NO₃:PO₄ ratios. Areas of low NO₃:PO₄ exist in the tropics proximal to water column denitrification zones. Releasing N₂ fixers from Fe limitation therefore increases the spatial coupling between N₂ fixation and water column denitrification and increases the global rate of N₂ fixation.

Figure 16 Change in ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ caused by a stronger coupling between N₂ fixation and tropical regions of low NO₃:PO₄ concentrations (i.e. tropical upwelling zones with active water column denitrification). Panel (a) shows the global zonal mean change, while panel (b) shows the average change in the euphotic zone, here defined as the top 100 m. Areas with very low NO₃ (<0.1 mmol m⁻³) are masked out.

We steadily decreased iron limitation ($K^{{\mathrm{D}}_{\mathrm{Fe}}}$) to increase the strength of spatial coupling between N₂ fixers and the tropical denitrification zones (Fig. 15). As N₂ fixers coupled more strongly to regions of low NO₃:PO₄, the rate of N₂ fixation increased from 122 to 144 to 154 Tg N yr⁻¹ (Table 6). An expansion of suboxia from 2.1 % to 2.5 % to 2.7 % of global ocean volume in the tropics accompanied the increase in N₂ fixation, as did a decrease in global mean δ¹³C_DIC of 0.06 ‰ and 0.1 ‰, since greater rates of N₂ fixation stimulated tropical export production. Due to the expansion of the already large suboxic zones, which occurred in both horizontal and vertical directions, the amount of organic carbon that reached tropical sediments (20^∘ S to 20^∘ N) increased from 0.35 to 0.46 to 0.51 Pg C yr⁻¹.

The overarching consequence for ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ due to an expansion of the suboxic zones was an increase in the sedimentary to water column denitrification ratio from 1.5 to 1.9 to 2.2, which decreased mean ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ from 5.6 ‰ to 5.2 ‰ to 5.0 ‰ (Table 6). The increase in N₂ fixation (${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{org}}=-\mathrm{1}\mathrm{‰}$) and sedimentary denitrification (${\mathit{ϵ}}_{\mathrm{sed}}^{{}^{\mathrm{15}}\mathrm{N}}=\mathrm{3}\mathrm{‰}$) in the tropics was felt globally for ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ (Fig. 16). Lower ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ by 0.5 ‰ and 0.9 ‰ permeated water columns in the Southern Ocean and tropics, respectively. Meanwhile, ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ was up to 10 ‰ lower in surface waters of the tropical and subtropical Pacific, which is where the greatest increase in N₂ fixation and sedimentary denitrification occurred. The dramatic reduction in surface ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ was subsequently conveyed to the sediments as δ¹⁵N_org±1 ‰–2 ‰.

These simple experiments demonstrate that the insights garnered from sedimentary records of δ¹⁵N are open to multiple lines of interpretation. An expansion of the suboxic zones, normally associated with an increase in ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ (Galbraith et al., 2013), could instead cause a decrease in ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ if more organic matter reached the sediments to stimulate sedimentary denitrification. There is good evidence that the suboxic zones might have undergone a vertical expansion (Hoogakker et al., 2018) and that more organic matter reached the tropical sediments under glacial conditions (Cartapanis et al., 2016). The glacial decrease in bulk δ¹⁵N_org recorded in the eastern tropical Pacific (Ganeshram et al., 1995; Liu et al., 2008) therefore does not necessarily mean a decrease in suboxia. Rather, our experiments show that lower δ¹⁵N_org might also be caused by an increase in local N₂ fixation and sedimentary denitrification. The decrease in ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ associated with more sedimentary denitrification and local N₂ fixation demonstrates the complexity of interpreting sedimentary δ¹⁵N_org records in the lower latitudes.

6 Conclusions

The stable isotopes of carbon (δ¹³C) and nitrogen (δ¹⁵N) are proxies that have been fundamental for understanding the ocean. We have included both isotopes into the ocean component of an Earth system model, CSIRO Mk3L-COAL, to enable future studies with the capability for direct model–proxy data comparisons. We detailed how these isotopes are simulated, how to conduct model–data comparisons using both water column and sedimentary data and some basic assessment of changes caused by altered ecosystem functioning. We made three overall findings. First, CSIRO Mk3L-COAL performs well alongside a number of isotope-enabled global ocean GCMs. Second, alteration of δ¹³C during formation of foraminiferal calcite does not jeopardise simple one-to-one comparisons with simulated δ¹³C of DIC, while diagenetic alteration of bulk organic δ¹⁵N during early burial must be accounted for in model–data comparisons. Third, changes in how marine ecosystems function can have significant and complex effects on δ¹³C and δ¹⁵N. Our idealised experiments hence showed that the interpretation of palaeoceanographic records may suffer from multiple lines of interpretation, particularly records from the lower latitudes where multiple processes imprint on the isotopic signatures laid down in sediments. Future work will involve palaeoceanographic simulations of CSIRO Mk3L-COAL that seek to understand how the oceanic carbon and nitrogen cycles respond to and influence important climate transitions.

Appendix A: Ecosystem component of the OBGCM

A1 Export production

A1.1 General phytoplankton group (^G)

The production of organic matter by the general phytoplankton group (${\mathrm{P}}_{\mathrm{org}}^{\mathrm{G}}$) is measured in units of mmol phosphorus (P) m⁻³ d⁻¹ and is dependent on temperature (T), nutrients (PO₄, NO₃ and Fe) and irradiance (I):

$$\begin{array}{ll}\text{(A1)}& {\displaystyle }& {\displaystyle }{\mathrm{P}}_{\mathrm{org}}^{\mathrm{G}}={{\mathrm{S}}^{\mathrm{G}}}_{\mathrm{E}:\mathrm{P}}\cdot \mathit{\mu }(T{)}^{\mathrm{G}}\cdot min({\mathrm{P}}_{\lim }^{\mathrm{G}},{\mathrm{N}}_{\lim }^{\mathrm{G}},{\mathrm{Fe}}_{\lim }^{\mathrm{G}},F\left(I\right)),{\displaystyle }& {\displaystyle }\text{where}\\ {\displaystyle }& {\displaystyle }{{\mathrm{S}}^{\mathrm{G}}}_{\mathrm{E}:\mathrm{P}}=\mathrm{0.005}\mathrm{mmol}{\mathrm{PO}}_{\mathrm{4}}{\mathrm{m}}^{\mathrm{3}}\\ \text{(A2)}& {\displaystyle }& {\displaystyle }\mathit{\mu }(T{)}^{\mathrm{G}}=\mathrm{0.59}\cdot {\mathrm{1.0635}}^T\text{(A3)}& {\displaystyle }& {\displaystyle }F\left(I\right)=\mathrm{1}-e^{L\left(I\right)}\text{(A4)}& {\displaystyle }& {\displaystyle }L\left(I\right)={\displaystyle \frac{I\cdot \mathit{\alpha }\cdot \mathrm{PAR}}{\mathit{\mu }\left(T\right)}}.\end{array}$$

In the example above, S_E:P converts growth rates in units of d⁻¹ to mmol PO₄ m⁻³ d⁻¹. S_E:P conceptually represents the export to production ratio, and for simplicity we assume it does not change. μ(T) is the temperature-dependent maximum daily growth rate of phytoplankton (doubling per day), as defined by Eppley (1972). The light limitation term (F(I)) is the productivity versus irradiance equation used to describe phytoplankton growth defined by Clementson et al. (1998) and is dependent on I, the daily averaged shortwave incident radiation (W m⁻²), α, the initial slope of the productivity versus radiance curve (d⁻¹ (W m⁻²)⁻¹) and PAR, the fraction of shortwave radiation that is photosynthetically active.

The nutrient limitation terms (${\mathrm{P}}_{\lim }^{\mathrm{G}}$, ${\mathrm{N}}_{\lim }^{\mathrm{G}}$ and ${\mathrm{Fe}}_{\lim }^{\mathrm{G}}$) may be calculated in two ways.

If the option for static nutrient limitation is true, then Michaelis–Menten kinetics (Dugdale, 1967) is used:

$$\begin{array}{}\text{(A5)}& {\displaystyle }& {\displaystyle }{\mathrm{P}}_{\lim }^{\mathrm{G}}={\displaystyle \frac{{\mathrm{PO}}_{\mathrm{4}}}{{\mathrm{PO}}_{\mathrm{4}}+K{}_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{G}}}}\text{(A6)}& {\displaystyle }& {\displaystyle }{\mathrm{N}}_{\lim }^{\mathrm{G}}={\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}}{{\mathrm{NO}}_{\mathrm{3}}+K{}_{{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{G}}}}\text{(A7)}& {\displaystyle }& {\displaystyle }{\mathrm{Fe}}_{\lim }^{\mathrm{G}}={\displaystyle \frac{\mathrm{Fe}}{\mathrm{Fe}+K{}_{\mathrm{Fe}}^{\mathrm{G}}}}.\end{array}$$

Half-saturation coefficients ($K^{{\mathrm{G}}_{\mathrm{nutrient}}}$) show a large range across phytoplankton species (e.g. Timmermans et al., 2004), and so for simplicity, we set $K^{{\mathrm{G}}_{{\mathrm{PO}}_{\mathrm{4}}}}=\mathrm{0.1}$ mmol PO₄ m⁻³ (Smith, 1982), $K^{{\mathrm{G}}_{{\mathrm{NO}}_{\mathrm{3}}}}=\mathrm{0.75}$ mmol NO₃ m⁻³ (Eppley et al., 1969; Carpenter and Guillard, 1971) and $K^{{\mathrm{G}}_{\mathrm{Fe}}}=\mathrm{0.1}$ µmol Fe m⁻³ (Timmermans et al., 2001).

If the option for variable nutrient limitation is true (default), then optimal uptake kinetics (Smith et al., 2009) is used:

$$\begin{array}{ll}\text{(A8)}& {\displaystyle }& {\displaystyle }{\mathrm{P}}_{\lim }^{\mathrm{G}}={\mathrm{PO}}_{\mathrm{4}}/({\displaystyle \frac{{\mathrm{PO}}_{\mathrm{4}}}{\mathrm{1}-f_A}}+{\displaystyle \frac{V/A}{f_A\cdot \mathrm{N}:\mathrm{P}}})\text{(A9)}& {\displaystyle }& {\displaystyle }{\mathrm{N}}_{\lim }^{\mathrm{G}}={\mathrm{NO}}_{\mathrm{3}}/({\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}}{\mathrm{1}-f_A}}+{\displaystyle \frac{V/A}{f_A}})\text{(A10)}& {\displaystyle }& {\displaystyle }{\mathrm{Fe}}_{\lim }^{\mathrm{G}}={\displaystyle \frac{\mathrm{Fe}}{\mathrm{Fe}+K_{\mathrm{Fe}}}},{\displaystyle }& {\displaystyle }\text{where}\\ {\displaystyle }& {\displaystyle }{f}_A=max\left[\right(\mathrm{1}+\sqrt{{\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{V/A}}}{)}^{-\mathrm{1}},\\ \text{(A11)}& {\displaystyle }& {\displaystyle }(\mathrm{1}+\sqrt{{\displaystyle \frac{\left[{\mathrm{PO}}_{\mathrm{4}}\right]\cdot \mathrm{N}:\mathrm{P}}{V/A}}}{)}^{-\mathrm{1}}].\end{array}$$

Optimal uptake kinetics varies the two terms in the denominator of the Michaelis–Menten form according to the availability of nutrients. It therefore accounts for different phytoplankton communities with different abilities for nutrient uptake and does so using the f_A term. The V∕A term represents the maximum potential nutrient uptake, V, over the cellular affinity for that nutrient, A, and is set at 0.1.

A1.2 Diazotrophs (^D; N₂ fixers)

Organic matter produced by diazotrophs (${\mathrm{P}}_{\mathrm{org}}^{\mathrm{D}}$) is also measured in units of mmol phosphorus (P) m⁻³ d⁻¹ and is calculated in the same form of Eq. (A1), but using the maximum growth rate μ(T)^D of Kriest and Oschlies (2015), notable changes in the limitation terms and minimum thresholds that ensure the nitrogen fixation occurs everywhere in the ocean, except under sea ice. ${\mathrm{P}}_{\mathrm{org}}^{\mathrm{D}}$ is calculated via

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{P}}_{\mathrm{org}}^{\mathrm{D}}={{\mathrm{S}}^{\mathrm{D}}}_{\mathrm{E}:\mathrm{P}}\cdot \mathit{\mu }(T{)}^{\mathrm{D}}\cdot max(\mathrm{0.01},min({\mathrm{N}}_{\lim }^{\mathrm{D}},{\mathrm{P}}_{\lim },{\mathrm{Fe}}_{\lim }^{\mathrm{D}}))\\ \text{(A12)}& {\displaystyle }& {\displaystyle }\cdot (\mathrm{1}-\mathrm{ico}),{\displaystyle }& {\displaystyle }\text{where}\\ \text{(A13)}& {\displaystyle }& {\displaystyle }\mathit{\mu }(T{)}^{\mathrm{D}}=max(\mathrm{0.01},-\mathrm{0.0042}{T}^{\mathrm{2}}+\mathrm{0.2253}T-\mathrm{2.7819})\text{(A14)}& {\displaystyle }& {\displaystyle }{\mathrm{N}}_{\lim }^{\mathrm{D}}=e^{-{\mathrm{NO}}_{\mathrm{3}}}\text{(A15)}& {\displaystyle }& {\displaystyle }{\mathrm{P}}_{\lim }^{\mathrm{D}}={\displaystyle \frac{{\mathrm{PO}}_{\mathrm{4}}}{{\mathrm{PO}}_{\mathrm{4}}+K^{{\mathrm{D}}_{{\mathrm{PO}}_{\mathrm{4}}}}}}\text{(A16)}& {\displaystyle }& {\displaystyle }{\mathrm{Fe}}_{\lim }^{\mathrm{D}}=max(\mathrm{0.0},\tanh (\mathrm{2}\mathrm{Fe}-K^{{\mathrm{D}}_{\mathrm{Fe}}}\left)\right).\end{array}$$

The half-saturation values for PO₄ and Fe limitation are set at 0.1 mmol m⁻¹ and 0.5 µmol m⁻¹, respectively, in the default parameterisation. The motivation for making N₂ fixers strongly limited by Fe was the high cellular requirements of Fe for diazotrophy (see Sohm et al., 2011, and references therein). A dependency on light is omitted from the limitation term when ${\mathrm{P}}_{\mathrm{org}}^{\mathrm{D}}$ is produced. The omission of light is justified by its strong correlation with sea surface temperature (Luo et al., 2014) and its negligible effect on nitrogen fixation in the Atlantic Ocean (McGillicuddy, 2014). Finally, the fractional area coverage of sea ice (ico) is included to ensure that cold water N₂ fixation (Sipler et al., 2017) does not occur under ice, since a light dependency is omitted.

A1.3 Calcifiers

The calcifying group produces calcium carbonate (CaCO₃) in units of mmol carbon (C) m⁻³ d⁻¹. The production of CaCO₃ is always a proportion of the organic carbon export of the general phytoplankton group (${\mathrm{C}}_{\mathrm{org}}^{\mathrm{G}}$), according to

$$\begin{array}{}\text{(A17)}& {\displaystyle }{\displaystyle }{\mathrm{CaCO}}_{\mathrm{3}}={\mathrm{C}}_{\mathrm{org}}^{\mathrm{G}}\cdot R_{{\mathrm{CaCO}}_{\mathrm{3}}}.\end{array}$$

The ratio of CaCO₃ to ${\mathrm{C}}_{\mathrm{org}}^{\mathrm{G}}$ ($R_{{\mathrm{CaCO}}_{\mathrm{3}}}$) can be calculated in two ways.

If the option for fixed $R_{{\mathrm{CaCO}}_{\mathrm{3}}}$ is true (default), then $R_{{\mathrm{CaCO}}_{\mathrm{3}}}$ is set to 0.08 as informed by the experiments of Yamanaka and Tajika (1996). The production of CaCO₃ is thus 8 % of ${\mathrm{C}}_{\mathrm{org}}^{\mathrm{G}}$ everywhere.

If the option for variable $R_{{\mathrm{CaCO}}_{\mathrm{3}}}$ is true, then $R_{{\mathrm{CaCO}}_{\mathrm{3}}}$ varies as a function of the saturation state of calcite (Ω_ca) according to Ridgwell et al. (2007), where

$$\begin{array}{}\text{(A18)}& {\displaystyle }{\displaystyle }{R}_{{\mathrm{CaCO}}_{\mathrm{3}}}=\mathrm{0.022}\cdot ({\mathrm{\Omega }}_{\mathrm{ca}}-\mathrm{1}{)}^{\mathit{\eta }}.\end{array}$$

The exponent (η) is easily modified consistent with the parameterisations of Zhang and Cao (2016) and controls the rate of CaCO₃ production at a given value of Ω_ca.

A2 Remineralisation

A2.1 General phytoplankton group (^G)

Organic matter produced by the general phytoplankton group (in units of phosphorus: ${\mathrm{P}}_{\mathrm{org}}^{\mathrm{G}}$) at the surface is instantaneously remineralised each time step at depth levels beneath the euphotic zone using a power law scaled to depth (Martin et al., 1987). This power law defines the concentration of organic matter remaining at a given depth (${\mathrm{P}}_{\mathrm{org}}^{\mathrm{G},z}$) as a function of organic matter at the surface (${\mathrm{P}}_{\mathrm{org}}^{\mathrm{G},\mathrm{0}}$) and depth itself (z). Its form is as follows:

$$\begin{array}{}\text{(A19)}& {\displaystyle }{\displaystyle }{\mathrm{P}}_{\mathrm{org}}^{\mathrm{G},z}={\mathrm{P}}_{\mathrm{org}}^{\mathrm{G},\mathrm{0}}\cdot ({\displaystyle \frac{z}{z_{\mathrm{rem}}}}{)}^b,\end{array}$$

where z_rem in the denominator represents the depth at which remineralisation begins and is set to be 100 m everywhere. The OBGCM therefore does not consider sinking speeds or an interaction between organic matter and physical mixing. However, variations in the b exponent affect the steepness of the curve, thereby emulating sinking speeds and affecting the transfer and release of nutrients from the surface to the deep ocean.

Remineralisation of ${\mathrm{P}}_{\mathrm{org}}^{\mathrm{G}}$ through the water column is therefore dependent on the exponent b value in Eq. (A19). The b exponent is calculated in two ways.

If the option for static remineralisation is true, then b is set to −0.858 according to Martin et al. (1987).

If the option for variable remineralisation is true (default), then b is dependent on the component fraction of picoplankton (F_pico) in the ecosystem. The F_pico shows a strong inverse relationship to the transfer efficiency (T_eff) of organic matter from beneath the euphotic zone to 1000 m depth (Weber et al., 2016). Because F_pico is not explicitly simulated in OBGCM, we estimate F_pico from the export production field in units of carbon (${\mathrm{C}}_{\mathrm{org}}^{\mathrm{G}}$), calculate T_eff using the parameterisation of Weber et al. (2016) and subsequently calculate the b exponent:

$$\begin{array}{}\text{(A20)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{pico}}=\mathrm{0.51}-\mathrm{0.26}\cdot {\displaystyle \frac{{\mathrm{C}}_{\mathrm{org}}^{\mathrm{G}}\left(\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}\right)}{{\mathrm{C}}_{\mathrm{org}}^{\mathrm{G},\max }\left(\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}\right)}}\text{(A21)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{eff}}=\mathrm{0.47}-\mathrm{0.81}\cdot F_{\mathrm{pico}}\text{(A22)}& {\displaystyle }& {\displaystyle }b={\displaystyle \frac{\log \left(T_{\mathrm{eff}}\right)}{\log \left(\frac{\mathrm{1000}}{\mathrm{100}}\right)}}=\log \left(T_{\mathrm{eff}}\right).\end{array}$$

A2.2 Diazotrophs (^D)

Remineralisation of diazotrophs (${\mathrm{P}}_{\mathrm{org}}^{\mathrm{D}}$) is calculated in the same way as the general phytoplankton group (${\mathrm{P}}_{\mathrm{org}}^{\mathrm{G}}$), with the exception that the depth at which remineralisation occurs is raised from 100 to 25 m in Eq. (A19). This alteration emulates the release of NO₃ from N₂ fixers well within the euphotic zone, which in some cases can exceed the physical supply from below (Capone et al., 2005). Release of their N- and C-rich organic matter (see stoichiometry Sect. A3.2) therefore occurs higher in the water column than the general phytoplankton group.

A2.3 Suboxic environments

The remineralisation of ${\mathrm{P}}_{\mathrm{org}}^{\mathrm{G}}$ and ${\mathrm{P}}_{\mathrm{org}}^{\mathrm{D}}$ will typically require O₂ to be removed, except for in regions where oxygen concentrations are less than a particular threshold (Den${}^{{{\mathrm{O}}_{\mathrm{2}}}_{\lim }}$), which is set to 7.5 mmol O₂ m⁻³ and represents the onset of suboxia. In these regions, the remineralisation of organic matter begins to consume NO₃ via the process of denitrification. We calculate the fraction of organic matter that is remineralised by denitrification (F_den) via

$$\begin{array}{}\text{(A23)}& {\displaystyle }{\displaystyle }{F}_{\mathrm{den}}=(\mathrm{1}-e^{-\mathrm{0.5}\cdot {\mathrm{Den}}_{\lim }^{{\mathrm{O}}_{\mathrm{2}}}}+e^{{\mathrm{O}}_{\mathrm{2}}-\mathrm{0.5}\cdot {\mathrm{Den}}_{\lim }^{{\mathrm{O}}_{\mathrm{2}}}}{)}^{-\mathrm{1}},\end{array}$$

such that F_den rises and plateaus at 100 % in a sigmoidal function as O₂ is depleted from 7.5 to 0 mmol m⁻³.

Following this, the strength of denitrification is reduced if the ambient concentration of NO₃ is deemed to be limiting. Denitrification within the modern oxygen minimum zones only depletes NO₃ towards concentrations between 15 and 40 mmol m⁻³ (Codispoti and Richards, 1976; Voss et al., 2001). Without an additional constraint that weakens denitrification as NO₃ is drawn down, here defined as r_den, NO₃ concentrations quickly go to zero in simulated suboxic zones (Schmittner et al., 2008). We weaken denitrification by prescribing a lower bound at which NO₃ can no longer be consumed via denitrification, ${\mathrm{Den}}_{\lim }^{{\mathrm{NO}}_{\mathrm{3}}}$, which is set at 30 mmol NO₃ m⁻³.

$$\begin{array}{}\text{(A24)}& {\displaystyle }& {\displaystyle }{r}_{\mathrm{den}}=\mathrm{0.5}+\mathrm{0.5}\cdot \tanh (\mathrm{0.25}\cdot {\mathrm{NO}}_{\mathrm{3}}-\mathrm{0.25}\cdot {\mathrm{Den}}^{{{\mathrm{NO}}_{\mathrm{3}}}_{\lim }}-\mathrm{2.5})\text{(A25)}& {\displaystyle }& {\displaystyle }\text{if}{F}_{\mathrm{den}}>r_{\mathrm{den}},\text{then}{F}_{\mathrm{den}}=r_{\mathrm{den}}\end{array}$$

F_den is therefore reduced if NO₃ is deemed to be limiting and subsequently applied against both ${\mathrm{P}}_{\mathrm{org}}^{\mathrm{G}}$ and ${\mathrm{P}}_{\mathrm{org}}^{\mathrm{D}}$ to get the proportion of organic matter to be remineralised by O₂ and NO₃.

If the availability of O₂ and NO₃ is insufficient to remineralise all the organic matter at a given depth level, z, then the unremineralised organic matter will pass into the next depth level. Unremineralised organic matter will continue to pass into lower depth levels until the final depth level is reached, at which point all organic matter is remineralised by either water column or sedimentary processes. This version of CSIRO Mk3L-COAL does not consider burial of organic matter.

A2.4 Calcifiers

The dissolution of CaCO₃ is calculated using an e-folding depth-dependent decay, where the amount of CaCO₃ at a given depth z is defined by

$$\begin{array}{}\text{(A26)}& {\displaystyle }{\displaystyle }{{\mathrm{CaCO}}_{\mathrm{3}}}^z={{\mathrm{CaCO}}_{\mathrm{3}}}^{\mathrm{0}}\cdot e^{\frac{-z}{z_{\mathrm{dis}}}},\end{array}$$

where z_dis represents the depth at which e⁻¹ of CaCO₃ (∼0.37) produced at the surface remains undissolved.

Calcifiers are not susceptible to oxygen-limited remineralisation or the concentration of carbonate ion because the dissolution of CaCO₃ depends solely on this depth-dependent decay. All CaCO₃ reaching the final depth level is remineralised without considering burial. Future work will include a full representation of carbonate compensation.

A3 Stoichiometry

The elemental constitution, or stoichiometry, of organic matter affects the biogeochemistry of the water column through uptake (production) and release (remineralisation). The general phytoplankton group and diazotrophs both affect carbon chemistry, O₂ and nutrients (PO₄, NO₃ and Fe), while the calcifiers only affect carbon chemistry tracers (DIC, DI¹³C and ALK).

Alkalinity ratios for both the general and nitrogen-fixing groups are the negative of the N:P ratio, such that for a loss of 1 mmol of NO₃, alkalinity will increase at 1 mmol eq m⁻³ (Wolf-Gladrow et al., 2007).

A3.1 General phytoplankton group (^G)

The stoichiometry of the general phytoplankton group is calculated in two ways.

If the option for static stoichiometry is true, then the $\mathrm{C}:\mathrm{N}:\mathrm{Fe}:\mathrm{P}$ ratio is set according to the Redfield ratio of $\mathrm{106}:\mathrm{16}:\mathrm{0.00032}:\mathrm{1}$ (Redfield et al., 1937).

If the option for variable stoichiometry is true (default), then the $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratio of ${\mathrm{P}}_{\mathrm{org}}^{\mathrm{G}}$ is made dependent on the ambient nutrient concentration according to Galbraith and Martiny (2015):

$$\begin{array}{}\text{(A27)}& {\displaystyle }& {\displaystyle }\mathrm{C}:\mathrm{P}=({\displaystyle \frac{\mathrm{6.9}\cdot \left[{\mathrm{PO}}_{\mathrm{4}}\right]+\mathrm{6}}{\mathrm{1000}}}{)}^{-\mathrm{1}}\text{(A28)}& {\displaystyle }& {\displaystyle }\mathrm{N}:\mathrm{C}=\mathrm{0.125}+{\displaystyle \frac{\mathrm{0.03}\cdot \left[{\mathrm{NO}}_{\mathrm{3}}\right]}{\mathrm{0.32}+\left[{\mathrm{NO}}_{\mathrm{3}}\right]}}\text{(A29)}& {\displaystyle }& {\displaystyle }\mathrm{N}:\mathrm{P}=\mathrm{C}:\mathrm{P}\cdot \mathrm{N}:\mathrm{C}.\end{array}$$

Thus, the stoichiometry of ${\mathrm{P}}_{\mathrm{org}}^{\mathrm{G}}$ varies across the ocean according to the nutrient concentration, and the uptake and release of carbon, nutrients and oxygen (see Sect. A3.4) are dependent on the concentration of surface PO₄ and NO₃. The ratio of iron to phosphorus (Fe:P) remains fixed at 0.00032, such that 0.32 µmol of Fe is consumed per mmol of PO₄. We chose to maintain a fixed Fe:P ratio because phytoplankton communities from subtropical to Antarctic waters appear to show similar iron content (Boyd et al., 2015), despite changes in $\mathrm{C}:\mathrm{N}:\mathrm{P}$. However, the ratio of $\mathrm{C}:\mathrm{N}:\mathrm{Fe}$ does change as a result of varying $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratios, with higher C:Fe in oligotrophic environments and lower C:Fe in eutrophic regions.

A3.2 Diazotrophs (^D)

The stoichiometry of diazotrophs is fixed at a $\mathrm{C}:\mathrm{N}:\mathrm{P}:\mathrm{Fe}$ ratio of $\mathrm{331}:\mathrm{50}:\mathrm{1}:\mathrm{0.00064}$, which represents values reported in the literature (Kustka et al., 2003; Karl and Letelier, 2008; Mills and Arrigo, 2010). Diazotrophs do not consume NO₃; rather, they consume N₂, which is assumed to be of unlimited supply, and release NO₃ during remineralisation.

A3.3 Calcifiers

Calcifying organisms produce CaCO₃, which includes DIC, DI¹³C and ALK, and these tracers are consumed and released at a ratio of $\mathrm{1}:\mathrm{0.998}:\mathrm{2}$, respectively, relative to organic carbon. Thus, the ratio of $\mathrm{C}:{\mathrm{DI}}^{\mathrm{13}}\mathrm{C}:\mathrm{ALK}$ relative to each unit of phosphorus consumed by the general phytoplankton group is equal to the rain ratio of CaCO₃ to organic phosphorus multiplied by $\mathrm{106}:\mathrm{105.8}:\mathrm{212}$. This group has no effect on nutrient tracers or oxygen values.

A3.4 Stoichiometry of remineralisation

The requirements for oxygen (${\mathrm{O}}_{\mathrm{2}}^{\mathrm{rem}}:\mathrm{P}$) and nitrate (${\mathrm{NO}}_{\mathrm{3}}^{\mathrm{rem}}:\mathrm{P}$) during oxic and suboxic remineralisation, respectively, are calculated from the $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratios of organic matter via the equations of Paulmier et al. (2009). Additional knowledge of the hydrogen and oxygen content of the organic matter is also required to calculate ${\mathrm{O}}_{\mathrm{2}}^{\mathrm{rem}}:\mathrm{P}$ and ${\mathrm{NO}}_{\mathrm{3}}^{\mathrm{rem}}:\mathrm{P}$. However, the hydrogen and oxygen content of phytoplankton depends strongly on the proportions of lipids, carbohydrates and proteins that constitute the cell. As there is no empirical model for predicting these physiological components based on environmental variables, we continue Redfield's legacy by assuming that all organic matter is a carbohydrate of the form CH₂O. Future work, however, should address this obvious bias.

To calculate ${\mathrm{O}}_{\mathrm{2}}^{\mathrm{rem}}:\mathrm{P}$ and ${\mathrm{NO}}_{\mathrm{3}}^{\mathrm{rem}}:\mathrm{P}$, we therefore need to first calculate the amount of hydrogen and oxygen in organic matter via

$$\begin{array}{}\text{(A30)}& {\displaystyle }& {\displaystyle }\mathrm{H}:\mathrm{P}=\mathrm{2}\mathrm{C}:\mathrm{P}+\mathrm{3}\mathrm{N}:\mathrm{P}+\mathrm{3}\text{(A31)}& {\displaystyle }& {\displaystyle }\mathrm{O}:\mathrm{P}=\mathrm{C}:\mathrm{P}+\mathrm{4}.\end{array}$$

Once a $\mathrm{C}:\mathrm{N}:\mathrm{P}:\mathrm{H}:\mathrm{O}$ ratio for organic matter is known, we calculate ${\mathrm{O}}_{\mathrm{2}}^{\mathrm{rem}}:\mathrm{P}$ and ${\mathrm{NO}}_{\mathrm{3}}^{\mathrm{rem}}:\mathrm{P}$ in units of mmol m⁻³ P⁻¹ using the equations of Paulmier et al. (2009):

$$\begin{array}{ll}{\displaystyle }{\mathrm{O}}_{\mathrm{2}}^{\mathrm{rem}}:\mathrm{P}=& {\displaystyle }-(\mathrm{C}:\mathrm{P}+\mathrm{0.25}\mathrm{H}:\mathrm{P}-\mathrm{0.5}\mathrm{O}:\mathrm{P}-\mathrm{0.75}\mathrm{N}:\mathrm{P}\\ \text{(A32)}& {\displaystyle }& {\displaystyle }+\mathrm{1.25})-\mathrm{2}\mathrm{N}:\mathrm{P}{\displaystyle }{\mathrm{NO}}_{\mathrm{3}}^{\mathrm{rem}}:\mathrm{P}=& {\displaystyle }-(\mathrm{0.8}\mathrm{C}:\mathrm{P}+\mathrm{0.25}\mathrm{H}:\mathrm{P}-\mathrm{0.5}\mathrm{O}:\mathrm{P}-\mathrm{0.75}\mathrm{N}:\mathrm{P}\\ \text{(A33)}& {\displaystyle }& {\displaystyle }+\mathrm{1.25})+\mathrm{0.6}\mathrm{N}:\mathrm{P}.\end{array}$$

The calculation of ${\mathrm{O}}_{\mathrm{2}}^{\mathrm{rem}}:\mathrm{P}$ accounts for the oxygen that is also needed to oxidise ammonium to nitrate.

From these calculations, we find the following requirements of oxic and suboxic remineralisation, assuming the static stoichiometry option for the general phytoplankton group:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{O}}_{\mathrm{2}}^{\mathrm{rem}}:{\mathrm{P}}_{\mathrm{org}}^{\mathrm{G}}=\mathrm{138}\\ {\displaystyle }& {\displaystyle }{\mathrm{NO}}_{\mathrm{3}}^{\mathrm{rem}}:{\mathrm{P}}_{\mathrm{org}}^{\mathrm{G}}=\mathrm{94.4}\\ {\displaystyle }& {\displaystyle }{\mathrm{O}}_{\mathrm{2}}\mathrm{rem}:{\mathrm{P}}_{\mathrm{org}}^{\mathrm{D}}=\mathrm{431}\\ {\displaystyle }& {\displaystyle }{\mathrm{NO}}_{\mathrm{3}}^{\mathrm{rem}}:{\mathrm{P}}_{\mathrm{org}}^{\mathrm{D}}=\mathrm{294.8}.\end{array}$$

These numbers change dynamically alongside $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratios when the stoichiometry of organic matter is allowed to vary.

A4 Sedimentary processes

The remineralisation of organic matter within the sediments is provided as an option in the OBGCM. Sedimentary denitrification, and its slight preference for the light isotopes of fixed nitrogen (${\mathit{ϵ}}_{\mathrm{sed}}^{{}^{\mathrm{15}}\mathrm{N}}$ =3 ‰), is an important component of the marine nitrogen cycle and its isotopes. It acts as an additional sink of NO₃ and reduces the δ¹⁵N value of the global ocean by offsetting the strong fractionation of water column denitrification (${\mathit{ϵ}}_{\mathrm{wc}}^{{}^{\mathrm{15}}\mathrm{N}}$ =20 ‰).

If sedimentary processes are active, the empirical model of Bohlen et al. (2012) is used to estimate the rate of sedimentary denitrification, where the removal of NO₃ is dependent on the rate of particulate organic carbon (${\mathrm{C}}_{\mathrm{org}}^{\mathrm{G}}$ + ${\mathrm{C}}_{\mathrm{org}}^{\mathrm{D}}$) arriving at the sediments and the ambient concentrations of oxygen and nitrate. In the following, we assume that the concentrations of NO₃ and O₂ that are available in the sediments are two-thirds of the concentration in the overlying water column based on observations of transport across the diffusive boundary layer (Gundersen and Jorgensen, 1990).

$$\begin{array}{}\text{(A34)}& {\displaystyle }\mathrm{\Delta }{\mathrm{NO}}_{\mathrm{3}}\left(\mathrm{sed}\right)=& {\displaystyle }(\mathit{\alpha }+\mathit{\beta }\cdot {\mathrm{0.98}}^{({\mathrm{O}}_{\mathrm{2}}-{\mathrm{NO}}_{\mathrm{3}})})\cdot ({\mathrm{C}}_{\mathrm{org}}^{\mathrm{G}}+{\mathrm{C}}_{\mathrm{org}}^{\mathrm{D}})\text{(A35)}& {\displaystyle }& {\displaystyle }\text{where }\mathit{\alpha }=\mathrm{0.04}\text{and}\mathit{\beta }=\mathrm{0.1}\end{array}$$

In the example above, both the α and β values were halved from the values of Bohlen et al. (2012) to raise global mean NO₃ concentrations and lower the sedimentary to water column denitrification ratio to between 1 and 2. If NO₃ is not available, the remaining organic matter is remineralised using oxygen if the environment is sufficiently oxygenated. An additional limitation is set for sediments underlying hypoxic waters (O₂<40 mmol m⁻³), where oxic remineralisation is weakened towards zero according to a hyperbolic tangent function $(\mathrm{0.5}+\mathrm{0.5}\cdot \tanh (\mathrm{0.2}\cdot {\mathrm{O}}_{\mathrm{2}}-\mathrm{5}\left)\right)$. If oxygen is also limiting, the remaining organic matter is remineralised via sulfate reduction. As sulfate is not explicitly simulated, we assumed that sulfate is always available to account for the remaining organic matter.

Thus, sedimentary denitrification is heavily dependent on the rate of organic matter arriving at the sediments. However, a large amount of sedimentary remineralisation is not captured using only these parameterisations because the coarse resolution of the OGCM enables it to resolve only the largest continental shelves, such as the shallow Indonesian seas. Many small areas of raised bathymetry in pelagic environments are also unresolved by the OGCM. To address this insufficiency and increase the global rate of sedimentation and sedimentary denitrification, we coupled a subgrid-scale bathymetry to the coarse-resolution OGCM following the methodology of Somes et al. (2013) using the Earth topography 5 min grid (ETOPO5) 1/12^∘ dataset. For each latitude-by-longitude grid point, we calculated the fraction of area that would be represented by shallower levels in the OGCM if this finer-resolution bathymetry were used. At each depth level above the deepest level, the fractional area represented by sediments on the subgrid-scale bathymetry can be used to remineralise all forms of exported matter (${\mathrm{C}}_{\mathrm{org}}^{\mathrm{G}}$, ${\mathrm{C}}_{\mathrm{org}}^{\mathrm{D}}$ and CaCO₃) via sedimentary processes.

Also following the methodology of Somes et al. (2013), we included an option to amplify sedimentary denitrification in the upper 250 m to account for narrow continental shelves that are not resolved by the OGCM. Narrow shelves experience strong rates of upwelling and productivity, and hence high rates of sedimentary denitrification (Gruber and Sarmiento, 1997). To amplify shallow rates of sedimentary denitrification, we included an optional acceleration factor (Γ_sed), set to 3.0 in the default parameterisation, dependent on the total fraction of shallower depths not covered by the subgrid-scale bathymetry:

$$\begin{array}{}\text{(A36)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }{\mathrm{NO}}_{\mathrm{3}}\left(\mathrm{sed}\right)=\mathrm{\Delta }{\mathrm{NO}}_{\mathrm{3}}\left(\mathrm{sed}\right)\cdot \left(\right(\mathrm{1}-F_{\mathrm{sgb}})\cdot {\mathrm{\Gamma }}_{\mathrm{sed}}+\mathrm{1}).\end{array}$$

For those grids with a low fraction covered by the subgrid-scale bathymetry (F_sgb), the amplification of sedimentary denitrification is therefore greatest.

Appendix B: Parameterisation of the OBGCM ecosystem component

Default parameters for the marine ecosystem component of CSIRO Mk3L-COAL are outlined in Tables B1, B2 and B3. The values presented in these tables are required as input when running the ocean model.

Table B1 Parameter values controlling export production in the ecosystem component of the CSIRO Mk3L-COAL ocean model. Default settings: Michaelis–Menten == False, Optimal Uptake == True, fix == True and Vary CaCO3 == False.

Table B2 Parameter values controlling remineralisation in the ecosystem component of the CSIRO Mk3L-COAL ocean model. Default settings: ReminPico == True, fix == True, den == True and sedfluxes == True.

Table B3 Parameter values controlling stoichiometry in the ecosystem component of the CSIRO Mk3L-COAL ocean model. Default settings: Vary Stoich == True, Vary frac13 == False and fix == True.