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 5 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 multi-millennial timescales, running at a rate of ∼400 years per day. We show that this coarse resolution, computationally efficient model adequately reproduces water column and coretop δC and δN measurements, making it a useful tool for palaeoceanographic research. Changes to ecosystem function involve varying phytoplankton stoichiometry, 10 varying CaCO3 production based on calcite saturation state, and varying N2 fixation via iron limitation. We find that large changes in CaCO3 production have little effect on δC and δN, while changes in N2 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, N2 fixation and nutrient utilisation influence δN. Hence, there is significant scope for isotope-enabled models to provide more 15 robust interpretations of the proxy records.


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 20 against one isotope relative to the other is minuscule, the isotopic content of a sample is conventionally expressed as a δ value (δ h E), where the ratio of the heavy to light element in solution ( h E: l E) is compared to a standard ratio ( h E std : l E std ) in units of per mille (‰).
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 h E for every 1000 units of l E at a hypothetical standard ratio ( h E std : l E std ) of 1:1. At more realistic standard ratios <<< 1:1, say 0.0112372:1 for a δ 13 C value of 0 ‰, a 5 fractionation at 10 ‰ would involve ∼0.0111123 0.010 · 0.0112372 1.0112372 units of 13 C per unit of 12 C. Slightly greater preference of one isotope over another in this case involves a preference for the lighter carbon isotope ( 12 C) over the heavier ( 13 C), which enriches the remaining dissolved inorganic carbon (DIC) in 13 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 (δ 13 C) and nitrogen (δ 15 N) have been fundamental for understanding 10 these important elements cycle within the ocean (e.g. Schmittner and Somes, 2016;Menviel et al., 2017;Rafter et al., 2017;Muglia et al., 2018). We will now briefly introduce each isotope in turn.
The distribution of δ 13 C is dependent on air-sea gas exchange, ocean circulation and organic matter cycling. These contributions make the δ 13 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 12 C over 13 C by 15 biology in surface waters enforces strong horizontal and vertical gradients in δ 13 C of DIC (δ 13 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 δ 13 C DIC reservoir over longer timescales . Because air-sea and speciation fractionation are temperature-dependent, such that cooler conditions tend to elevate the δ 13 C DIC of surface waters, they also tend to smooth the gradients produced by biology 20 by working antagonistically to them. Despite this smoothing, biological fractionation drives strong gradients at the surface, which imparts unique δ 13 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., 2017;Muglia et al., 2018). δ 15 N is determined by biological processes that add or remove fixed forms of nitrogen. It therefore records the relative rates 25 of sources and sinks within the marine nitrogen cycle (Brandes and Devol, 2002). Dinitrogen (N 2 ) 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 δ 15 N of approximately -1 ‰ . Denitrification is the largest sink of fixed nitrogen and occurs in deoxygenated water columns and sediments. Denitrification fractionates strongly against 15 N at ∼25 ‰ (Cline and Kaplan, 1975). Fractionation during denitrification is most strongly expressed in the water column where ample nitrate (NO 3 ) is available, making water 30 column denitrification responsible for elevating global mean δ 15 N above the -1 ‰ of N 2 fixers (Brandes and Devol, 2002).
Meanwhile, denitrification occurring in the sediments only weakly fractionates against 15 N (Sigman et al., 2009), providing only a slight enrichment of δ 15 N above that introduced by N 2 fixation. Variations in δ 15 N can therefore tell us about global changes in the ratio of sedimentary to water column denitrification, with increases in δ 15 N associated with increases in the proportion of denitrification occurring in the water column , but it can also reflect regional changes in N 2 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 δ 15 N more complicated. Basically, when nitrogen is abundant the preference for 14 N over 15 N increases but when nitrogen is limited this preference disappears (Altabet and Francois, 2001). Complete utilisation of nitrogen therefore reduces fractionation to 0 ‰. 5 While this adds complexity, it also imbues δ 15 N as a proxy of nutrient utilisation by phytoplankton. As nitrogen supply to phytoplankton is controlled by physical delivery from below, changes in δ 15 N can be interpreted as changes in the physical supply . Phytoplankton fractionate against 15 N at ∼5 ‰ (Wada, 1980) when bioavailable nitrogen is abundant. If nitrogen is utilised to completion, which occurs in much of the low to mid latitude ocean, then no fractionation will occur and the δ 15 N of organic matter will reflect the δ 15 N of the nitrogen that was supplied . However, in the case 10 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 15 N as phytoplankton preferentially consume 14 N.
As the remaining bioavailable N is continually enriched in 15 N the organic matter that settles into sediments beneath a zone of incomplete nutrient utilisation will bear this enriched δ 15 N signal. In combination with modelling (Schmittner and Somes, 2016), the δ 15 N record is able to provide evidence for a more efficient utilisation of bioavailable nitrogen during glacial times 15 (Martinez-Garcia et al., 2014;Kemeny et al., 2018) and a less efficient one during the Holocene .
Complimentary measurements of δ 13 C and δ 15 N provide powerful, multi-focal insights into oceanographic processes. δ 13 C is largely a reflection of how water masses mix away the strong vertical and horizontal gradients enforced by biology, while δ 15 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 considera- 20 ble uncertainty because there are multiple processes that imprint on the measured values. Our goal is to equip version 1.0 of the CSIRO Mk3L-COAL Earth system model with oceanic δ 13 C and δ 15 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 25 first test of the model, we take the opportunity to document how changes in ecosystem functioning affect δ 13 C and δ 15 N.

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 Coupled Model Inter- 30 comparison Project Phase 3 (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 Lenton, 2014, 2018). These studies demonstrate that the model can 5 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 metres at the surface and 450 metres in the deep ocean (Table 1). The OGCM timestep is one hour, while the OBGCM timestep is 1 day. The ability of the OBGCM to 10 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.

Ocean biogeochemical model (OBGCM)
The OBGCM is equipped with 13 prognostic tracers ( Figure 1). These can be grouped into carbon chemistry fields, oxygen fields, nutrient fields, age tracers and nitrous oxide (N 2 O). Carbon chemistry fields include dissolved inorganic carbon (DIC), 15 alkalinity (ALK), DI 13 C and radiocarbon ( 14 C). Radiocarbon is simulated according to Toggweiler et al. (1989). Oxygen fields include dissolved oxygen (O 2 ) and abiotic dissolved oxygen (O abio 2 ), a purely physical tracer from which true oxygen utilisation (TOU) can be calculated (Duteil et al., 2013). Nutrient fields include phosphate (PO 4 ), dissolved bioavailable iron (Fe), nitrate (NO 3 ) and 15 NO 3 . Although we define the phosphorus and nitrogen tracers as their dominant species, being PO 4 and NO 3 , these tracers can also be thought of as total dissolved inorganic phosphorus and nitrogen pools. Remineralisation, 20 for instance, implicitly accounts for the process of nitrification from ammonium (NH 4 ) to NO 3 (Paulmier et al., 2009) and therefore implicitly includes NH 4 and nitrite (NO 2 ) within the NO 3 tracer. Age tracer fields include years since subduction from the surface (Age gbl ), and years since entering a suboxic zone where O 2 concentrations are less than 10 mmol m −3 (Age omz ). Finally, N 2 O in µmol m −3 is produced via nitrification and denitrification according to the temperature-dependent equations of Freing et al. (2012). All air-sea gas exchanges (CO 2 , 13 CO 2 O 2 and N 2 O) and carbon speciation reactions are 25 computed according to the Ocean Modelling Intercomparison Project phase 6 protocol (Orr et al., 2017).
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 appendix A. Default parameters for the OBGCM are further provided in appendix B. Briefly, the ecosystem model simulates the production, remineralisation and stoichiometry 30 (elemental composition) of three types of primary producers: a general phytoplankton group, diazotrophs (N 2 fixers) and calcifiers. (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 N2 fixers (diazotrophs), while particulate inorganic carbon (PIC) as calcium carbonate (CaCO3) is produced by calcifiers. Export of POC by the general (G) phytoplankton group and N2 fixers (D; diazotrophs) are herein referred to as C G org and C D org (see appendix A1), respectively. (4) Remineralisation of sinking organic matter under oxic and suboxic conditions. (5) Sedimentary oxic and suboxic remineralisation. (6) Nitrous oxide production and consumption.

δ 13 C
The OBGCM explicitly simulates the fractionation of 13 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 12 C. Fractionation occurs during air-sea gas exchange, equilibrium reactions and biological consumption in the euphotic zone.

5
The air-sea gas exchange of 13 CO 2 is calculated as the exchange of CO 2 with additional fractionation factors applied to the sea-air and air-sea components (Zhang et al., 1995;Orr et al., 2017). The flux of 13 CO 2 across the air-sea interface, F ( 13 CO 2 ), therefore takes the form of CO 2 with additional terms that convert to units of 13 C in both environments. Without any isotopic fractionation, the equation requires the gas piston velocity of carbon dioxide in m s −1 (k CO2 ), the concentration of aqueous CO 2 in both mediums at the air-sea interface in mmol m −3 (CO air 2 and CO sea 2 ), and the ratios of 13 C: 12 C in both mediums 10 (R atm and R sea ): where, A transfer of 13 C into the ocean is therefore positive and an outgassing is negative. The R atm is set to a preindustrial atmospheric δ 13 C of -6.48 ‰ (Friedli et al., 1986).
The fractionation of carbon isotopes during air-sea exchange involves three components. These are ( C k ) a kinetic fractionation that occurs during transfer of gaseous CO 2 into or out of the ocean, ( C aq←g ) a fractionation that occurs as gaseous CO 2 becomes aqueous CO 2 (is dissolved in solution), and ( C DIC←g ) an equilibrium isotopic fractionation as carbon speciates 20 into dissolved inorganic carbon (DIC) constituents (H 2 CO 2 ⇔ HCO − 3 ⇔ CO 2− 3 ). The kinetic fractionation during transfer, C k , is constant at 0.99912, thus reducing the δ 13 C of carbon entering the ocean by 0.88 ‰. Conversely, carbon outgassing increases the δ 13 C of the ocean. The fractionation during dissolution ( C aq←g ) and speciation ( C DIC←g ) are both dependent on temperature. Fractionation during speciation is also dependent on the fraction of CO 2− 3 relative to total DIC (f CO 2− 3 ). These fractionation factors are parameterised as: Dissolution of CO 2 into seawater ( C aq←g ) therefore preferences the lighter isotope and lowers δ 13 C by between 1.32 and 1.14 ‰, while the speciation of gaseous CO 2 into DIC instead preferences the heavier isotope and raises δ 13 C by between 10.7 and 6.8 ‰ for temperatures between -2 and 35 • C. 30 These fractionation factors are applied to the gaseous exchange of CO 2 (Eq. (2)) to calculate carbon isotopic fractionation.
Because fractionation to aqueous CO 2 from DIC ( C aq←DIC ) is equal to , a strong preference to hold the heavy isotope in solution exists ( C aq←DIC = -11.9 to -7.9 ‰ between -2 and 35 • C). Aqueous carbon that is transferred to the atmosphere is hence depleted in 13 C. It is therefore the equilibrium fractionation associated with carbon speciation that is largely responsible for bolstering the oceanic δ 13 C signature above the atmospheric signature, as it tends to shift 13 C towards the oxidised species (CO 2− 3 ), a tendency that strengthens under cooler conditions.

5
In the default version of CSIRO Mk3L-COAL v1.0, the fractionation of carbon during biological uptake ( 13 C bio ) 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 2 concentration (CO 2 (aq) in mmol m −3 ) and the growth rate (µ in day −1 , as a function of temperature and limiting resources) following Tagliabue and Bopp (2008): An upper bound of 25 ‰ exists within Eq. (6) when µ CO2(aq) approaches zero, but a lower bound does not. We chose to limit 13 C bio to a minimum of 15 ‰ given the reported variations of 13 C bio from culture studies (e.g. Laws et al., 1995).
Biological fractionation of 13 C is then applied to the uptake and release of organic carbon. 15 Because biological fractionation is strong for the general phytoplankton group, which dominates export production throughout most of the ocean, this imparts a negative δ 13 C signature to the deep ocean. Subsequent remineralisation releases DIC with no fractionation. Finally, the concentration of DI 13 C is converted into a δ 13 C via: 20 where 0.0112372 is the Pee Dee Belemnite standard (Craig, 1957).

δ 15 N
The OBGCM explicitly simulates the fractionation of 15 N from the pool of bioavailable nitrogen. For simplicity we treat this bioavailable pool as nitrate (NO 3 ), where total NO 3 is the sum of 15 NO 3 and 14 NO 3 . We therefore chose to ignore fractionation during reactions involving ammonium, nitrite and dissolved organic nitrogen, which can vary in their isotopic composition 25 independent of NO 3 but represent a small fraction of the bioavailable pool of nitrogen.
The isotopic signatures of N 2 fixation and atmospheric deposition, and the fractionation during water column denitrification N 2 fixation and atmospheric deposition introduce 15 NO 3 to the ocean with δ 15 N values of -1 ‰ and -2 ‰, respectively, while biological assimilation, water column denitrification and sedimentary denitrification fractionate against 15 NO 3 at 5 ‰, 20 ‰ and 3 ‰, respectively   Fig. 1).
The accepted standard 15 N: 14 N ratio used to measure variations in nature is the average atmospheric 15 N: 14 N ratio of 0.0036765. To minimise numerical errors caused by the OGCM, we set the atmospheric standard to 1. This scales up the 5 15 NO 3 such that a δ 15 N value of 0 ‰ was equivalent to an 15 N: 14 N ratio of 1:1.
Because we simulate NO 3 and 15 NO 3 as tracers, our calculations require solving for an implicit pool of 14 NO 3 during each reaction involving 15 NO 3 . The introduction of NO 3 at a fixed δ 15 N N O3 of -1 ‰ due to remineralisation of N 2 fixer biomass provides a simple example with which we can begin to describe our equations. Setting the isotopic value of newly fixed NO 3 to -1 ‰ is simple because it removes any complications associated with fractionation. We note, however, that in reality the 10 nitrogenase enzyme does fractionate during its conversion of aqueous N 2 (+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 δ 15 N of N 2 fixer biomass equal to -1 ‰, which reflects the biomass of N 2 fixers associated with the more common Mo-nitrogenase. 15 A δ 15 N N O3 of -1 ‰ is equivalent to a 15 N: 14 N ratio of 0.999 in our approach where 0 ‰ equals a 1:1 ratio of 15 N: 14 N. If the amount of NO 3 being added is known alongside its 15 N: 14 N ratio, in this case 0.999 for N 2 fixation, we are able to calculate how much 15 NO 3 is added. We begin with two equations that describe the system.
Ultimately, we need to solve for the change in 15 NO 3 associated with an introduction of NO 3 by N 2 fixation. Our knowns are the change in NO 3 , the δ 15 N of that NO 3 , and the 15 N std / 14 N std . Our two unknowns are 15 NO 3 and 14 NO 3 . We must solve 20 for 14 NO 3 implicitly by describing it according to 15 NO 3 by rearranging Eq. (11).
This allows us to replace the 14 NO 3 term in Eq. (10), such that In our example of N 2 fixation we know the δ 15 N of the newly added NO 3 as being -1 ‰. We also know 15 N std / 14 N std as equal to 1:1, or 1. Our equation is simplified.
We can now solve for 15 NO 3 by rearranging the equation.
The same calculation is applied to NO 3 addition via atmospheric deposition except at a constant fraction of 0.998 (δ 15 N = -2 ‰), and can be applied to any addition or subtraction of 15 NO 3 relative to NO 3 where the isotopic signature is known. ) involves more considerations because we must account for the preference of 14 NO 3 over 15 NO 3 . We begin with an of 5 ‰ for biological assimilation. This is equivalent to an 15 NO 3 : 14 NO 3 ratio of 0.995 when our atmospheric 5 standard is equal to 1:1 using the following equation.
Note that a positive value returns an 15 NO 3 : 14 NO 3 ratio < 1, while a negative δ 15 N in the previous example with N 2 fixation also returned an 15 NO 3 : 14 NO 3 ratio < 1. This works because the reactions are in opposite directions. N 2 fixation adds NO 3 , while assimilation removes NO 3 . This means that 0.995 units of 15 NO 3 are assimilated per unit of 14 NO 3 . As we have seen, a more useful way to quantify this is per unit of NO 3 assimilated into organic matter. Using Eq. (15), we find that ∼0.4987 However, we must also account for the effect that NO 3 availability has on fractionation. The preference of 14 NO 3 over 15 NO 3 strongly depends on the availability of NO 3 , such that when NO 3 is abundant the preference for the lighter isotope 15 will be strongest. This preference (fractionation) becomes weaker as NO 3 is depleted because cells will absorb any NO 3 that is available irrespective of its isotopic composition (Mariotti et al., 1981). Thus, as NO 3 is utilised, u, towards 100 % of its availability (u = 1), the fractionation against 15 NO 3 decreases to an of 0 ‰. This means that when u is equal to 1, no fractionation occurs and equal parts 15 N and 14 N (0.5:0.5 per unit NO 3 ) are assimilated. As we are interested in long timescales, we chose the accumulated product equations (Altabet and Francois, 2001) to approximate this process, where: 20 u = min 0.999, max 0.001, 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 15 N: 14 N of NO 3 in the reactant pool to determine the 15 N: 14 N of the product. In this case, it is the 15 N: 14 N of newly created organic matter, but could also be unused NO 3 effluxed from denitrifying cells in the case for denitrification.
We then solve for how much 15 NO 3 is assimilated into organic matter using Eq. (15) because we now know the change in NO 3

25
(∆N O 3 ) and the 15 N: 14 N of the product, which is 15 N org / 14 N org in our example of biological assimilation.
Here, the change in 15 NO 3 is equivalent to that assimilated into organic matter. Following assimilation into organic matter, the release of 15 NO 3 through the water column during remineralisation occurs with no fractionation, such that the same δ 15 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 5 , the isotopic ratio of the reactant, the amount of reactant that is used, and the total amount of reactant available.

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., 2016Buchanan et al., , 2018. Rather than reproduce these studies, we concentrate here on how the biogeochemical model performs relative to measurements 10 of δ 13 C and δ 15 N in the water column (Eide et al., 2017, δ 15 N N O3 data courtesy of The Sigman Lab at Princeton University) and in the sediments 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 appendix B. Each experiment was run towards steady-state under preindustrial atmospheric 15 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.  Menviel et al. (2017). The isotope-enabled Community Earth System Model (iCESM) fields for δ 13 C (low resolution) provided by Alexandra Jahn (Jahn et al., 2015) and those for δ 15 N (high resolution) provided by Simon Yang (Yang and Gruber, 2016 4.1 δ 13 C of dissolved inorganic carbon (δ 13 C DIC ) The recent reconstruction of preindustrial δ 13 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 δ 13 C DIC from CSIRO Mk3L-COAL broadly replicated the preindustrial distribution. The predicted 5 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 RMS error of 0.42 ‰, while a greater degree of disagreement in the values of δ 13 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 δ 13 C DIC was too low in the tropics of the major basins by ∼0.2 ‰, and too high in the 10 North Pacific and North Atlantic by 0.4 to 0.6 ‰ (Fig. 3).
These inconsistencies were likely related to physical and biological limitations of CSIRO Mk3L-COAL. δ 13 C DIC in subsurface tropical waters was too low because restricted horizontal mixing and high carbon export drove very negative δ 13 C DIC values. The very negative δ 13 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 15 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 (appendix A), where remineralisation is prevented when O 2 and NO 3 are unavailable. Organic matter therefore falls deeper into the interior through oxygen-deficient zones, leading to their vertically expansion. Almost certainly, however, it was the poorly represented dynamics within the Pacific basin that were responsible for high δ 13 C DIC in the subsurface North 20 Pacific, which contains low O 2 , low δ 13 C DIC water due to northward transport from the tropics.
Another inconsistency was a positive bias in the upper 200-500 metres, with values exceeding 2 ‰in many areas of the lower latitudes. However, values as high as 2 ‰ have been measured in the upper 500 metres 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 2 ) it is possible that the upper ocean values of Eide et al. (2017) underestimate the preindustrial δ 13 C DIC surface field. 25 It is also equally possible that a fixed biological fractionation ( 13 C bio ) of 21 ‰ may have driven unrealistic enrichment in the simulated field. High growth rates, such as occurs in the tropical regions, 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 13 C bio at 21 ‰, we implemented biological fractionation that is dependent on phytoplankton growth rate and aqueous CO 2 concentration (Eq. 6). We found the implementation of a variable 13 bio reduced high values in the upper part of the low latitude ocean, but that 30 this reduction was small (Fig. 4). The overwhelming effect was an increase in δ 13 C DIC throughout the interior, itself caused by weaker fractionation in the tropical ocean. Global mean δ 13 C DIC subsequently increased by 0.   It is helpful to place our predicted δ 13 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 δ 13 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., 2017;Tagliabue et al., 2009;Schmittner and Somes, 2016). Predicted δ 13 C DIC performs adequately in CSIRO Mk3L-COAL relative to these state-of-the-art models.

5
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 2 , is similar to CSIRO Mk3L-COAL when this form of fractionation was activated (vary-13 C bio ). PISCES and iCESM-low were the best performing models, equally demonstrating high correlations, low biases, accurate 10 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 δ 13 C DIC ≥ 2.0 ‰ (Supplementary Figures S1, S2, S3 and S4) 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 δ 13 C once more suggests that the upper ocean values between 200 and 500 metres of Eide et al. (2017) may be too low. The underestimation of δ 13 C DIC may be due to a neglect of biology introducing anthropogenic, isotopically-5 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 δ 13 C (Schmittner et al., 2017) is perhaps a more accurate representation of preindustrial δ 13 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 10 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 δ 13 C DIC .
We extended our assessment of modelled δ 13 C DIC by comparing it to a compilation of benthic δ 13 C measured within the cal- 15 cite of foraminifera from the genus Cibicides (Schmittner et al., 2017), a genus on which much of the palaeoceanographic δ 13 C records are based. For this comparison, we adjusted our predicted δ 13 C DIC to predicted δ 13 C Cib using the linear dependence on carbonate ion concentration and depth suggested by Schmittner et al. (2017): This adjustment accounts for slight fractionation during incorporation of DIC into foraminiferal calcite and is found to be partly explained by the concentration of CO 2− 3 ions and pressure. A one to one comparison between δ 13 C DIC and δ 13 C Cib hence 20 introduces some degree of error since this fractionation is not accounted for. Because we are interested in applying simulated δ 13 C DIC to a palaeoceanographic context, we must first be able to convert our simulated δ 13 C DIC to δ 13 C Cib in an effort to make better comparisons, particularly as the distribution of CO 2− 3 is subject to change. By adjusting our three-dimensional δ 13 C DIC output using Eq. (21), we attain predicted δ 13 C Cib (see inset entitled "Calibration" in Fig. 5). For good measure, we also computed measures of statistical fit for a traditional one to one comparison between δ 13 C DIC and δ 13 C Cib to assess the 25 benefit of the calibration.
Measured δ 13 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 1,763 to 690, lowered the mean of measured δ 13 C Cib from 0.76 ‰ to 0.52 ‰ and reduced the absolute range from -0.9→2.1 to -0.7→2.1. 30 Adjusted δ 13 C Cib using Eq. (21) showed good fit to measured δ 13 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 δ 13 C DIC and δ 13 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 ( 13 C bio ) had little effect except to increase values and slightly improve measures of 5 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 δ 13 C DIC with δ 13 C Cib measurements (Schmittner et al., 2017). Some notable over and underestimation occurred in the adjusted δ 13 C Cib output that more or less mirrored those inconsistencies previously discussed for δ 13 C DIC . Values as low as -1.9 ‰, well below measured δ 13 C Cib minima of -0.7 ‰, existed  We produced univariate measures of fit by comparing measurements of δ 15 N N O3 with equivalent values from CSIRO Mk3L-COAL at the nearest point ( Fig. 6; Table 4). Measured δ 15 N N O3 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.   (Marconi et al., 2017). A basin-wide rate of Atlantic N 2 fixation equal to ∼33 Tg N yr −1 lowered Atlantic values below 5 ‰ and was fundamental for reproducing the observations. Outside the Atlantic  (e and f). Colour shading represents the density of data, such that the darker a mass of data points is the more data is represented there.
where data is more sparse, the model successfully reproduced the meridional gradients across the Antarctic, Subantarctic and Subtropical zones, the subsurface δ 15 N N O3 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). 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 δ 13 C DIC that were transported too deeply into the North Atlantic  . 6). Meanwhile, the deep (> 1,500 metres) Eastern Tropical Pacific tended to overestimate the data, owing to a large, deep, unimodal suboxic zone. These physically-driven inconsistencies in the oxygen field are common to other coarse resolution models Schmittner et al., 2008), and like the δ 13 C distribution, were the main cause of the misfit between simulated and observed δ 15 N N O3 . 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 5 iCESM-high (Table 1). This comparison demonstrated that the modelled distribution of δ 15 N N O3 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 Figure 6). Atlantic δ 15 N N O3 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 ( Figure 6). UVic-MOBI underestimated the data, but produced high 10 correlations in the Southern Ocean and globally. Regionally, PISCES was best correlated to the measurements of δ 15 N N O3 of the three models, although it had a consistent positive bias. iCESM-high was acceptably correlated to 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 δ 15 N N O3 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-

15
COAL should therefore take these limitations into account. Overall, however, we find that CSIRO Mk3L-COAL broadly column denitrification is shut off (section A2.3). In CSIRO Mk3L-COAL this is set at 30 mmol m −3 , which is an arbitrary limit that was implemented to prevent water column denitrification from reducing NO 3 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 10 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 δ 15 N of organic matter (δ 15 N org ) CSIRO Mk3L-COAL tracks the δ 15 N signature of organic matter (δ 15 N org ) that is deposited in the sediments. We compared the simulated δ 15 N org to the coretop compilation of Tesdal et al. (2013) with 2,176 records of δ 15 N org . These records were binned and averaged onto the CSIRO Mk3L-COAL ocean grid, such that the 2,176 records became 592. When comparing sediment coretop measurements of δ 15 N to that of the model, it is necessary to consider how δ 15 N org is altered by early burial.

5
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 δ 15 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 10 is made, alongside an attempt to account for the diagenetic offset using two depth-dependent corrections (Table 5 and Fig. 10): The first correction (δ 15 N cor:1 org ) is taken from Robinson et al. (2012), while the second (δ 15 N cor:2 org ) 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).
Following binning and averaging onto the model grid, the raw comparison immediately showed a consistent underestimation 15 of the coretop 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 coretop data in Fig. 10. Like the nitrogen isotope model of Somes et al. (2010), we find that the offset between simulated and observed coretop bulk δ 15 N org is roughly equivalent to the observed average diagenetic offset of ∼2.3 ± 1.8 ‰. This indicates that diagenetic alteration of δ 15 N org is active during early burial in the coretop data. 20 Including a diagenetic offset therefore improved agreement between our predicted δ 15 N org and the coretop data considerably (Table 5 and Fig. 10). Both corrections accounted for the enrichment of δ 15 N in deeper regions and the minor diagenetic alteration in areas of high sedimentation that typically occur in shallower sediments. The average δ 15 N org increased to 4.5 ‰ for δ 15 N cor:1 org and 5.2 ‰ for δ 15 N cor:2 org . 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 δ 15 N org were well replicated 25 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 coretop data is encouraging for future study. We suggest that future palaeoceanographic model-data comparisons of δ 15 N org use the depth-correction of Schmittner and Somes (2016) as it provided the best correlations and reproduced Southern Ocean δ 15 N org at 0.5 ‰ greater than the global mean (see Table 5).  22) and (23)) that account for diagenetic alteration.

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

Variable versus Redfieldian stoichiometry
Enabling variable stoichiometry (see appendix A3) of the general phytoplankton group (P G org ) over a Redfieldian ratio (C:N:P:O rem 2 :NO rem 3 = 106:16:1:-138:-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 4 , such as the Atlantic Ocean (Fig. 11a), which elevated O 2 and NO 3 demand 10 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 −1 . Approximately 0.1 Pg C yr −1 , or 25 % of the increase, was attributed purely to organic carbon export from N 2 fixation, which increased from 107 to 122 Tg N yr −1 as higher N:P ratios in the tropics broadened their competitive niche. The total contribution of N 2 fixation to the increase 15 in carbon export was likely greater than 25 %, as NO 3 also became more available to NO 3 -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 20 δ 13 C DIC increased from 0.52 to 0.54 ‰, and δ 15 N N O3 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 δ 13 C DIC by up to 0.1 ‰ and δ 15 N N O3 Figure 11. Simulated difference in the C:P ratio of exported organic matter due to variable stoichiometry as compared to Redfield stoichiometry (top) and the resulting change in carbon export out of the euphotic zone (bottom). \U 1 by up to 1 ‰, while the deep ocean, particularly the Pacific, was enriched in both isotope to a similar degree (Fig. 12). Depletion of both isotope in waters subducted between 40 • S 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 3 in the heavier isotope to the same degree as surface waters travelled north. Enrichment of δ 13 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 δ 15 N in the deep ocean was the result of increased tropical production that 5 increased water column denitrification ( 15 N wc = 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 2 < 10 mmol m −3 ) by 0.5 %. The increase in water column denitrification caused by the expansion of suboxia increased δ 15 N N O3 , while the local increase in carbon export that drove the increase in water column 10 denitrification reduced δ 13 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 δ 15 N N O3 increased and δ 13 C DIC decreased.

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 3 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 3 yr −1 . Now we investigate how spatial variations in the CaCO 3 :C org ratio (R CaCO3 in Eq. (A17)) affected δ 13 C DIC and δ 13 C Cib (see appendix A1.3).

5
We applied three different values of η to Eq. (A18) to alter the quantity of CaCO 3 produced per unit of organic carbon (C G org ) 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).
Mean R CaCO3 was 4.5, 6.6 and 9.5 % and annual CaCO 3 production was 0.32, 0.47 and 0.68 Pg CaCO 3 yr −1 in the three experiments. Although different in total CaCO 3 production, the three experiments shared the same spatial patterns. Low 10 latitude waters were high in R CaCO3 , 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 CaCO3 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 3 when compared to an enforced 8 % rain Table 6. Summary of the biogeochemical effects of the different treatments of the ecosystem in CSIRO Mk3L-COAL. Corg is the total organic carbon exported from the euphotic zone composed of both general and diazotrophic phytoplankton groups (C G org + C D org ; see appendix A1), while CCaCO 3 is the total export of CaCO3 out of the euphotic zone. The sum of Corg and CCaCO 3 equal 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 δ 13 CDIC is higher than reported in Table 2  ratio, while CaCO 3 production between 40 • S and 40 • N relative to a fixed R CaCO3 of 8 % was dependent on η. The highest η coefficient of 1.09 achieved greater export of CaCO 3 in the mid to lower latitude regions of high export production (Fig. 13).
The consequence of increasing CaCO 3 production in the mid-lower latitudes was a loss of upper ocean alkalinity, subsequent outgassing of CO 2 and losses in the DIC inventory. Losses in global DIC were 95 and 130 Pg C as R CaCO3 increased from 4.6 → 6.6 → 9.5 % (Table 6), equivalent to 1 5 th of the glacial increase in oceanic carbon (Ciais et al., 2011).

5
Despite the significant changes associated with the implementation of Ω ca -dependent CaCO 3 production, effects were negligible on both δ 13 C DIC and δ 13 C Cib . Global mean δ 13 C DIC was 0.51 ‰, when R CaCO3 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 δ 13 C Cib was 0.59 ‰, when R CaCO3 was fixed at 8 %, and this changed to 0.60, 0.58 and 0.55 ‰. Minimal change in δ 13 C Cib indicated minimal change in the CO 2− 3 concentration (see Eq. (21)), which varied by ≤ 2 mmol m −3 between experiments. Visual inspection of the change 10 in δ 13 C DIC and δ 13 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 3 production, magnified the enrichment. Enrichment of δ 13 C DIC and δ 13 C Cib was caused by outgassing of CO 2 as surface alkalinity decreased in response to greater CaCO 3 production (Fig. 14). The change, however, was at most 0.1 ‰, which lies well within one standard deviation of variability known in the proxy data Figure 13. Global distribution of CaCO3 export as a percentage of organic carbon (Corg) export (top), and the change in the CaCO3 production field as a result of making CaCO3 production dependent on calcite saturation state (η = 1.09) compared to when it was a fixed 8 % of Corg (bottom). Areas where export production does not occur due to severely nutrient limited conditions are masked out.  (Schmittner et al., 2017). We therefore find little scope for recognising even large variations in global CaCO 3 production (0.32 to 0.68 Pg CaCO 3 yr −1 ) 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 3 burial or dissolution from the sediments according the calcite saturation state of overlying water (Boudreau, 2013). To neglect ocean-sediment CaCO 3 cycling is to neglect of an important aspect of the global carbon cycle active on millennial timescales (Sigman et al., 2010). Changes in 5 CaCO 3 burial and dissolution could have a non-negligible effect on δ 13 C through altering whole ocean alkalinity and thereby air-sea gas exchange of CO 2 , which would in turn affect surface δ 13 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.

Strength of coupling between N 2 fixation and denitrification
The degree to which N 2 fixers are spatially coupled to the tropical denitrification zones is controlled by altering the degree to 10 which N 2 fixers are limited by iron (K D F e ) in Eq. (A12) (see appendix A1.2). Decreasing K D F e ensures that N 2 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 3 :PO 4 ratios. Areas of low NO 3 :PO 4 exist in the tropics proximal to water column denitrification Figure 14. Changes in the distribution of carbon isotopes (δ 13 CDIC and δ 13 C Cib ; top) and carbon chemistry (dissolved inorganic carbon and alkalinity; bottom) as a result of increasing CaCO3 production in surface waters between 40 • S and 40 • N.
zones. Releasing N 2 fixers from Fe limitation therefore increases the spatial coupling between N 2 fixation and water column denitrification and increases the global rate of N 2 fixation.
We steadily decreased iron limitation (K D F e ) to increase the strength of spatial coupling between N 2 fixers and the tropical denitrification zones (Fig. 15). As N 2 fixers coupled more strongly to regions of low NO 3 :PO 4 , the rate of N 2 fixation increased from 122 to 144 to 154 Tg N yr −1 (Table 6). An expansion of suboxia from 2.1 to 2.5 to 2.7 % of global ocean volume in 5 the tropics accompanied the increase in N 2 fixation, as did a decrease in global mean δ 13 C DIC of 0.06 and 0.1 ‰, since greater rates of N 2 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 −1 .
The overarching consequence for δ 15 N N O3 due to an expansion of the suboxic zones was an increase in the sedimentary to 10 water column denitrification ratio from 1.5 to 1.9 to 2.2, which decreased mean δ 15 N N O3 from 5.6 to 5.2 to 5.0 ‰ ( Table 6).
The increase in N 2 fixation (δ 15 N org = -1 ‰) and sedimentary denitrification ( 15 N sed = 3 ‰) in the tropics was felt globally for δ 15 N N O3 (Fig. 16). Lower δ 15 N N O3 by 0.5 and 0.9 ‰ permeated water columns in the Southern Ocean and tropics, respectively. Meanwhile, δ 15 N N O3 was up to 10 ‰ lower in surface waters of the tropical and subtropical Pacific, which is where the greatest increase in N 2 fixation and sedimentary denitrification occurred. The dramatic reduction in surface δ 15 N N O3 was subsequently conveyed to the sediments as δ 15 N org ± 1-2 ‰.
These simple experiments demonstrate that the insights garnered from sedimentary records of δ 15 N are open to multiple lines of interpretation. An expansion of the suboxic zones, normally associated with an increase in δ 15 N N O3 , could instead cause a decrease in δ 15 N N O3 if more organic matter reached the sediments to stimulate sedimentary 5 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 δ 15 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 δ 15 N org might also be caused by an increase in local N 2 fixation and sedimentary denitrification. The decrease in δ 15 N N O3 associated with more sedimentary denitrification 10 and local N 2 fixation demonstrates the complexity of interpreting sedimentary δ 15 N org records in the lower latitudes.

Conclusions
The stable isotopes of carbon (δ 13 C) and nitrogen (δ 15 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 15 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 δ 13 C during formation of foraminiferal calcite does not jeopardise simple one to one comparisons with simulated δ 13 C of DIC, while diagenetic alteration of bulk organic δ 15 N during early burial must be accounted for in model-data comparisons. Third, changes in how marine ecosystems function 20 can have significant and complex effects on δ 13 C and δ 15 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.

25
Data availability. All model output is provided for download on Australia's National Computing Infrastructure (NCI) at https://geonetwork.nci.org.au/geon and is citable with doi:10.25914/5c6643f64446c. Nitrogen isotope data are available by request to Dario M. Marconi and Daniel M. Sigman at Princeton University. LOVECLIM data is freely available for download at https://researchdata.ands.org.au/loveclim-glacial-maximum-d13c-d14c/792249. UVic-MOBI data was provided by Christopher Somes, PISCES data by Laurent Bopp, iCESM-high data from Simon Yang and iCESM-low data by Alexandra Jahn.

NO3
Code availability. The source code for CSIRO Mk3L-COAL is shared via a repository located at http://svn.tpac.org.au/repos/CSIRO_Mk3L/branches/CSIR COAL/. Access to the repository may be obtained by following the instructions at https://www.tpac.org.au/csiro-mk3l-access-request/.
Access to the source code is subject to a bespoke license that does not permit commercial usage, but is otherwise unrestricted. An "outof-the-box" run directory is also available for download with all files required to run the model in the configuration used in this study, although users will need to modify the runscript according to their computing infrastructure. The production of organic matter by the general phytoplankton group (P G org ) is measured in units of mmol phosphorus (P) m −3 day −1 , and is dependent on temperature (T), nutrients (PO 4 , NO 3 , and Fe) and irradiance (I): where, In the above, S E:P converts growth rates in units of day −1 to mmol PO 4 m −3 day −1 . 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 (doublings day −1 ), 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 −2 ), α, the initial slope of the productivity versus 10 radiance curve (day −1 /(W m −2 )), and P AR, the fraction of shortwave radiation that is photosynthetically active.
The nutrient limitation terms (P G lim , N G lim , and F e G lim ) may be calculated in two ways. If the option for static nutrient limitation is true, then Michaelis-Menten kinetics (Dugdale, 1967) is used: Half-saturation coefficients (K G nutrient ) show a large range across phytoplankton species (e.g. Timmermans et al., 2004), and so for simplicity, we set K G P O4 = 0.1 mmol PO 4 m −3 (Smith, 1982), K G N O3 = 0.75 mmol NO 3 m −3 (Eppley et al., 1969;15 Carpenter and Guillard, 1971) and K G F e = 0.1 µmol Fe m −3 (Timmermans et al., 2001).
If the option for variable nutrient limitation is true (default), then Optimal Uptake kinetics (Smith et al., 2009) is used: where, 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 2 fixers)
Organic matter produced by diazotrophs (P D org ) is also measured in units of mmol phosphorus (P) m −3 day −1 , 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. P D org is calculated via: where, The half saturation values for PO 4 and Fe limitation are set at 0.1 mmol m −1 and 0.5 µmol m −1 , respectively, in the default parameterisation. The motivation for making N 2 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 P D org 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 15 of sea ice (ico) is included to ensure that cold water N 2 fixation (Sipler et al., 2017) does not occur under ice, since a light dependency is omitted.
The calcifying group produces calcium carbonate (CaCO 3 ) in units of mmol carbon (C) m −3 day −1 . The production of CaCO 3 is always a proportion of the organic carbon export of the general phytoplankton group (C G org ), according to: The ratio of CaCO 3 to C G org (R CaCO3 ) can be calculated in two ways. If the option for fixed R CaCO3 is true (default), then R CaCO3 is set to 0.08 as informed by the experiments of Yamanaka 5 and Tajika (1996). The production of CaCO 3 is thus 8 % of C G org everywhere. If the option for variable R CaCO3 is true, then R CaCO3 varies as a function of the saturation state of calcite (Ω ca ) according to Ridgwell et al. (2007), where: The exponent (η) is easily modified consistent with the parameterisations of Zhang and Cao (2016) and controls the rate of CaCO 3 production at a given value of Ω ca .

A2.1 General phytoplankton group ( G )
Organic matter produced by the general phytoplankton group (in units of phosphorus: P G org ) at the surface is instantaneously remineralised each timestep 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 (P G,z org ) as a function of organic matter 15 at the surface (P G,0 org ) and depth itself (z). Its form is as follows: Where z rem in the denominator represents the depth at which remineralisation begins and is set to be 100 metres everywhere.
The OBGCM therefore does not consider sinking speeds, nor 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. 20 Remineralisation of P G org 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 ef f ) of organic matter from 25 beneath the euphotic zone to 1,000 metres 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 (C G org ), calculate T ef f using the parameterisation of Weber et al. (2016), and subsequently calculate the b exponent: Remineralisation of diazotrophs (P D org ) is calculated in the same way as the general phytoplankton group (P G org ), with the exception that the depth at which remineralisation occurs is raised from 100 to 25 metres in Eq. (A19). This alteration emulates the release of NO 3 from N 2 fixers well within the euphotic zone, which in some cases can exceed the physical supply from 5 below (Capone et al., 2005). Release of their N and C-rich organic matter (see Stoichiometry section A3.2) therefore occurs higher in the water column than the general phytoplankton group.

A2.3 Suboxic environments
The remineralisation of P G org and P D org will typically require O 2 to be removed, except for in regions where oxygen concentrations are less than a particular threshold (Den O2 lim ), which is set to 7.5 mmol O 2 m −3 and represents the onset of suboxia. 10 In these regions, the remineralisation of organic matter begins to consume NO 3 via the process of denitrification. We calculate the fraction of organic matter that is remineralised by denitrification (F den ) via: Such that F den rises and plateaus at 100 % in a sigmoidal function as O 2 is depleted from 7.5 to 0 mmol m −3 .
Following this, the strength of denitrification is reduced if the ambient concentration of NO 3 is deemed to be limiting.
Denitrification within the modern oxygen minimum zones only depletes NO 3 towards concentrations between 15 and 40 mmol 15 m −3 (Codispoti and Richards, 1976;Voss et al., 2001). Without an additional constraint that weakens denitrification as NO 3 is drawn down, here defined as r den , NO 3 concentrations quickly go to zero in simulated suboxic zones . We weaken denitrification by prescribing a lower bound at which NO 3 can no longer be consumed via denitrification, Den N O3 lim , which is set at 30 mmol NO 3 m −3 .
if F den > r den , then F den = r den (A25) F den is therefore reduced if NO 3 is deemed to be limiting, and subsequently applied against both P G org and P D org to get the 20 proportion of organic matter to be remineralised by O 2 and NO 3 .
If the availability of O 2 and NO 3 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.
The dissolution of CaCO 3 is calculated using an e-folding depth-dependent decay, where the amount of CaCO 3 at a given depth z is defined by: Where z dis represents the depth at which e −1 of CaCO 3 (∼0.37) produced at the surface remains undissolved.
Calcifiers are not susceptible to oxygen-limited re-mineralisation nor the concentration of carbonate ion because the disso-5 lution of CaCO 3 depends solely on the this depth-dependent decay. All CaCO 3 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, 10 O 2 , and nutrients (PO 4 , NO 3 and Fe), while the calcifiers only affect carbon chemistry tracers (DIC, DI 13 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 3 , alkalinity will increase at 1 mmol Eq m −3 (Wolf- Gladrow et al., 2007).

A3.1 General phytoplankton group ( G )
The stoichiometry of the general phytoplankton group is calculated in two ways. 15 If the option for static stoichiometry is true, then the C:N:Fe:P ratio is set according to the Redfield ratio of 106:16:0.00032:1 (Redfield et al., 1937).
If the option for variable stoichiometry is true (default), then the C:N:P ratio of P G org is made dependent on the ambient nutrient concentration according to Galbraith and Martiny (2015): The stoichiometry of diazotrophs is fixed at a C:N:P:Fe ratio of 331:50:1: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 3 , rather they consume N 2 , which is assumed to be of unlimited supply, and release NO 3 during remineralisation.

5
Calcifying organisms produce CaCO 3 , which includes DIC, DI 13 C and ALK, and these tracers are consumed and released at a ratio of 1:0.998:2, respectively, relative to organic carbon. Thus, the ratio of C:DI 13 C:Alk relative to each unit of phosphorus consumed by the general phytoplankton group is equal to the rain ratio of CaCO 3 to organic phosphorus multiplied by 106:105.8:212. This group has no effect on nutrient tracers or oxygen values. The calculation of O rem 2 :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 stoi- :P D org = 294.8 These numbers change dynamically alongside C:N:P ratios when the stoichiometry of organic matter is allowed the vary.

A4 Sedimentary processes
The remineralisation of organic matter within the sediments is provided as an option in the OBGCM. Sedimentary denitrifica- the sediments and the ambient concentrations of oxygen and nitrate. In the following, we assume that the concentrations of NO 3 and O 2 that are available in the sediments are 2 3 of the concentration in overlying water column based on observations of transport across the diffusive boundary layer (Gundersen and Jorgensen, 1990).
∆N O 3 (sed) = α + β · 0.98 O2−N O3 · C G org + C D org (A34) where, α = 0.04 and β = 0.1 (A35) In the above, both the α and β values were halved from the values of Bohlen et al. (2012) to raise global mean NO 3 concentrations and lower the sedimentary to water column denitrification ratio to between 1 and 2. If NO 3 is not available, the 15 remaining organic matter is remineralised using oxygen if the environment is sufficiently oxygenated. An additional limitation is set for sediments underlying hypoxic waters (O 2 < 40 mmol m −3 ), where oxic remineralisation is weakened towards zero according to a hyperbolic tangent function (0.5 + 0.5 · tanh(0.2 · O 2 − 5)). 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. 20 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 sub-grid scale bathymetry to the course resolution 25 OGCM following the methodology of Somes et al. (2013) using the ETOPO5 1 12 th of a degree 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 sub-grid scale bathymetry can be used to remineralise all forms of exported matter (C G org , C D org and CaCO 3 ) via sedimentary processes.
Also following the methodology of Somes et al. (2013), we included an option to amplify sedimentary denitrification in the 5 upper 250 metres 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 sub-grid scale bathymetry: For those grids with a low fraction covered by the sub-grid scale bathymetry (F sgb ), the amplification of sedimentary denitrifi-10 cation 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 A1, A2, and A3. The values presented in these tables are required as input when running the ocean model. Finally, the lead author is indebted to an Australian Fulbright postgraduate scholarship, which supported him at the Princeton Geosciences department during the writing of this manuscript.