Variable C/P composition of organic production and its effect on ocean carbon storage in glacial–like model simulations

ocean carbon storage in glacial–like model simulations Malin Ödalen1, Jonas Nycander1, Andy Ridgwell2,3, Kevin I. C. Oliver4, Carlye D. Peterson2, and Johan Nilsson1 1Department of Meteorology, Bolin Centre for Climate Research, Stockholm University, 106 91 Stockholm, Sweden 2Department of Earth Sciences, University of California–Riverside, Riverside, CA 92521, USA 3School of Geographical Sciences, Bristol University, Bristol BS8 1SS, UK 4National Oceanography Centre, Southampton, University of Southampton, Southampton SO14 3ZH, United Kingdom Correspondence: Malin Ödalen (malin.odalen@misu.su.se)

1 Introduction 20 During the last four glacial maxima, atmospheric CO 2 (henceforth pCO atm 2 was lowered by ∼ 90-100 ppm compared to the interglacials (Petit et al., 1999;Lüthi et al., 2008). Due to the difference in size between the oceanic, terrestrial and atmospheric carbon reservoirs, where the oceanic reservoir is by far the largest with >90% of their summed carbon contents (reviewed by, Ciais et al., 2013), it is likely that most of the CO 2 that was removed from the atmosphere was stored in the glacial ocean. In addition, studies of paleoproxy records indicate that carbon storage in the glacial terrestrial biosphere was smaller compared to in interglacial climate (Shackleton, 1977;Duplessy et al., 1988;Curry et al., 1988;Crowley, 1995;Adams and Faure, 1998;Ciais et al., 2012;Peterson et al., 2014). During deglaciation, radiocarbon evidence indicate that CO 2 was rapidly released 5 from the ocean back to the atmosphere (Marchitto et al., 2007;Skinner et al., 2010).
Numerous processes, both physical and biological, have been identified as possible contributers to increased glacial oceanic storage. As glacial climate was substantially colder than interglacial climate; the global averages of surface, and ocean temperature at the LGM are estimated to have been 3-8°C and 2.0-3.2°C colder, respectively, than the pre-industrial (Stocker, 2014;Headly and Severinghaus, 2007;Bereiter et al., 2018). Due to the temperature effect on solubility, a colder ocean can 10 hold more carbon. However, glacial changes in salinity partly offset the temperature effect on solubility (reviewed by, Kohfeld and Ridgwell, 2009). A colder climate is also drier, and the dry conditions led to increased glacial dust deposition compared to interglacial climate (Mahowald et al., 2006). It has been hypothesised that the addition of dust contributed to increased iron availability in the surface ocean, and that this contributed to a strengthening the biological sequestration of carbon in the glacial ocean compared to interglacials (Martin, 1990). The addition of iron would allow for more complete usage of other nutrients in 15 regions where iron is limiting for biological production. Such strengthening of the retention of biologically-sourced carbon in the deep ocean (or the so-called biological pump), through changes in nutrient availability, light conditions, and/or ocean circulation, has long been considered an important player in the glacial increase in ocean CO 2 storage (Broecker, 1982a;Sarmiento and Toggweiler, 1984;Archer et al., 2000;Sigman and Boyle, 2000). Other studies have pointed to changes in carbonate preservation in coral reefs and deep-sea marine sediments (Berger, 1982;Broecker, 1982a;Archer and Maier-Reimer, 1994). 20 It is likely that reduced ventilation of the deep water, through changes in ocean circulation and expanded sea ice cover acting as a barrier for air-sea gas exchange, contributed to increasing the glacial ocean carbon retention (e.g., Boyle and Keigwin, 1987;Duplessy et al., 1988;Stephens and Keeling, 2000;Marchitto and Broecker, 2006;Adkins, 2013;Menviel et al., 2017;Skinner et al., 2017). Model studies by (Menviel et al., 2017) show that reduced Southern Hemisphere westerly winds produce reduced ventilation of Antarctic Bottom Water (AABW) in line with evidence from proxy records of δ 13 C . In addition, it has 25 been shown that strengthening of the winds over the Southern Ocean was a likely contributer to deglacial outgassing of CO 2 from the ocean to the atmosphere (Mayr et al., 2013). Extensive summaries of the processes responsible for high glacial ocean carbon storage, and examples of their interactions, are given by Brovkin et al. (2007); Kohfeld and Ridgwell (2009) ;Hain et al. (2010); Sigman et al. (2010). Despite the efforts of identifying the responsible processes, models have been struggling to achieve the full lowering of pCO atm 2 expected for a glacial. 30 In this paper, we focus on the biological pump and how it responds to glacial-like changes in climate. Our aim is to investigate how the level of simplification of the biological carbon uptake in an Earth system model model may affect the glacial drawdown of pCO atm 2 . Most biogeochemical models used in glacial climate studies have a simple representation of biological production, which assumes that carbon, C, and inorganic nutrients such as phosphorus, P, are taken up in fixed proportion to GitHub repository, 2019). We run a version of cGENIE with the same phosphorus plus iron (Fe) co-limitation scheme as used in the iron cycle model inter-comparison study of Tagliabue et al. (2016). The model branch enabled for use with flexible C:P ratios (see 2.2) is tagged as release v0.9.5, and the model configurations used in this paper are included in this release (cGENIE release v0.9.5, 2019, see Code availability for details).

5
In the original version of the cGENIE Earth system model (Ridgwell et al., 2007), as well as in the version of Tagliabue et al. (2016), the stoichiometric ratios are based on Redfield (1963). Thus, there is a fixed relationship between the number of moles of the elements that are taken up (positive) or released (negative) during production of organic matter in the ocean. This relationship is P : C : O 2 = 1 : 106 : −138, where O 2 is dissolved oxygen (nitrogen is assumed only implicitly for the purpose of accounting for organic matter creation and remineralisation related alkalinity transformations in the ocean (Ridgwell et al., 10 2007)). An exception is iron, where Fe:C varies as a function of iron availability as described in Watson et al. (2000).
Although the average elemental composition of organic matter in the ocean is close to the Redfield ratios, the stoichiometry of production of new organic material has shown high in-situ variability. Variability occurs between species, but also within the same species, and has been shown to depend on environmental factors such as nutrient availability, water temperature and light (e.g., Le Quéré et al., 2005b;Galbraith and Martiny, 2015;Yvon-Durocher et al., 2015;Tanioka and Matsumoto, 2017; 15 Garcia et al., 2018a;Moreno et al., 2018). We test the importance of this variability for glacial ocean CO 2 storage by running (1) This relation shows that, when [P O 4 ] is low, organisms bind more C per atom of P than they do under high [P O 4 ] conditions ( Fig. 1). Eq. 1 is applied in cGENIE in the calculations of biological C uptake at the surface ocean based on the surface concentration of P O 4 . In Section 3.1.2, Eq. 1 is also used to translate model surface P O 4 fields to the corresponding surface C/P ratios for the organic matter produced in each grid cell. Note that, while we change the ratio C/P, the ratio C/O 2 remains 25 the same in all experiments. As a result, the P/O 2 ratio changes between experiments.

Experiments
We start all experiments from an interglacial/modern control state, which has been run for 10,000 years to steady state, using either Redfield (Ctrl RED ) or variable (Ctrl GAM ) stoichiometry. The control states have a prescribed pCO atm 2 of 278 ppm and the same climate (Table 1), but due to the differences in C/P, they have different ocean carbon inventories (Table S.1.). 30 In Ctrl GAM , the export flux of organic matter (see Ridgwell et al., 2007) has a global average C/P composition of 121/1 and thus the global ocean carbon storage is larger than in Ctrl RED . This also suggest that a perturbation, which increases ocean storage of P through the biological pump, could cause storage of 15 (i.e. 121-106) more carbon atoms in simulations using GAM compared to RED, simply because the average composition of the formed biological material is different. To distinguish between the role of the flexibility of the stoichiometry and the change in the mean composition of organic material, we add a control state with fixed stoichiometry where C/P = 121/1 (henceforth denoted Ctrl 121 ).

5
In order to explore the effects of variable stoichiometry, we make a sensitivity study where we apply changes to boundary conditions, individually and in combination (see Section 2.3.3), that may be representative of changes that occurred during glacial periods. All experiments are listed in Table 1.
The applied changes in boundary conditions are physical perturbations (colder climate):

10
radiative forcing corresponding to LGM CO 2 = 185 ppm zonal albedo profile representative of LGM (calculated from the LGM climate simulation of Davies-Barnard et al.,

15
biological perturbations: changed remineralisation length scale (e.g., Matsumoto, 2007;Chikamoto et al., 2012;Menviel et al., 2012) increased dust forcing, as simulated for LGM (regridded from Mahowald et al., 2006) By applying the above perturbations, we aim to approach, but not fully resolve, some of the characteristics of the Last Glacial Maximum (LGM) ocean, which appears to have had a global average ocean temperature (T oce ) 2.57 ±0.24°C colder than the 20 Holocene (Bereiter et al., 2018), a weakly ventilated deep ocean (e.g., Menviel et al., 2017) and a more efficient biological pump (e.g., Sarmiento and Toggweiler, 1984;Martin, 1990;Sigman and Boyle, 2000). We also aim to increase carbon retention in the deep ocean (Muglia et al., 2018). Bottom Water (AABW) cell compared to Ctrl. Colder conditions achieve a stronger solubility pump, thereby strengthening the retention of carbon in the deep ocean. As the physical perturbations affect the ocean circulation and temperature, they thereby affect the nutrient distribution, and the rates of nutrient upwelling and biological growth (slower growth in colder water). They thereby affect the biological productivity (Ridgwell et al., 2007).
The biological perturbations serve to achieve a more efficient biological pump, which is connected with increased retention 30 of nutrients and carbon in the deep ocean, and lower surface nutrient concentrations in productive regions (see Sections 3.2.3 and 3.2.4). With flexible stoichiometry, lower surface nutrient concentrations results in a higher C/P ratio, further increasing the export production, and therby the carbon retention in the deep ocean. In our experiments, we show that the flexible stoichiometry amplifies the response of the biological pump to both physical and biological perturbations.
The perturbations and the experiments are described in detail in Sections 2.3.1-2.3.3.

Physical perturbations 5
We change the physical conditions for climate by changing radiative forcing and albedo to LGM-like conditions and denote these changes LGM phy. We set the radiative forcing in the model to correspond to an atmosphere with 185 ppm CO 2 instead of 278. However, we allow the pCO atm 2 to freely evolve (starting from the value of 278 ppm of the Ctrl state atmosphere) in response to the cooler climate. For albedo, we apply a zonal LGM albedo profile (calculated from the LGM climate simulation of Davies-Barnard et al., 2017). Assumptions of a simple zonal profile, instead of a 2D field re-gridded from 10 PMIP LGM simulations, allows for a better consistency with the original zonal mean albedo profile developed for the modern configuration of GENIE (Marsh et al., 2011). Together, the changes in radiative forcing and albedo causes the global ocean average temperature (T oce ) to decrease by 2.1°C compared to Ctrl (see Section3.2.1).
To achieve a longer residence time of the AABW water mass, and an associated increase in carbon and nutrient retention, we apply weaker winds (denoted W N A × 0.5). We use the Southern Ocean wind profile of Lauderdale et al. (2013), where the 15 peak westerly wind strength at 50°S has been halved compared to the control state (see Section 2.3.3). The winds north and south of the peak are reduced accordingly to give a continuous profile (see Fig. 2 a of Lauderdale et al., 2013). The result is a weaker overturning (see Table 2) and a longer residence time of the AABW as may be expected for the glacial ocean (Menviel et al., 2017;Skinner et al., 2017) (see also Section 4.5). Thus, this approach is justifiable in a model of reduced complexity.
However, there are studies suggesting that the Southern Ocean winds may in fact have been stronger during glacial times (Sime 20 et al., 2013;Kohfeld et al., 2013;Sime et al., 2016). To avoid an expansion of the NADW overturning cell that would be inconsistent with the glacial ocean Lynch-Stieglitz et al., 1999;Curry and Oppo, 2005;Marchitto and Broecker, 2006;Hesse et al., 2011), winds north of 35°N are also gradually reduced so that the wind strength north of 50°N is reduced by half compared to the control state. In cGENIE, gas transfer velocities are calculated as a function of wind speed (described in Ridgwell et al., 2007), and following Wanninkhof (1992). Consequently, weaker winds also lead to reduced gas 25 exchange with the atmosphere.

Biological perturbations
In the ocean, phytoplankton growth rates and remineralisation of particulate organic carbon are processes that both work more slowly at colder temperatures (Eppley, 1972;Laws et al., 2000). Cooling of the ocean would thus lead to decreased production of particulate organic matter (POC), and simultaneously to a slower degradation of POC, with competing effects on export 30 production (i.e. the amount of C captured by primary production that leaves the surface ocean without being remineralised) (Matsumoto, 2007). However, Matsumoto (2007) shows that the effect of slower remineralisation dominates the effect on export production. It has therefore been hypothesised that the cooling of the glacial ocean led to a deepening of the remineralisation length scale (henceforth denoted RLS) in the ocean, and thereby more efficient retention of organic carbon in the deep ocean (Matsumoto, 2007;Chikamoto et al., 2012), which in turn caused a lowering of pCO atm 2 . Menviel et al. (2012) find that such a deepening results in model changes in export production in poor agreement with paleo-proxies, while Chikamoto et al.
(2012) find improved model agreement with the glacial proxy records of export production and stable carbon isotopes for temperature-dependent growth rates and remineralisation. Deeper remineralisation also results in increased nutrient retention 5 in the deep ocean, thus causing changes in surface nutrient fields and in C/P ratios of GAM . We test the effect of changes in RLS by multiplying the model default RLS by a factor Y (RLS × Y,see Section 2.3.3 ).
Increased dust forcing leads to increased iron (F e) availability. This allows for increased productivity (and hence more efficient usage of other nutrients) in the high-nutrient, low-chlorophyll (HNLC) regions in the North Pacific, Equatorial Pacific and Southern Ocean, where iron (Fe) is the limiting micronutrient (Martin, 1990;Moore et al., 2013). The variable stoichiometry 10 in GAM is expected to be influential if the concentrations of P decrease in such regions as a result of increased Fe availability.
This process may hence be of importance in a glacial scenario where dust forcing increases as a result of the drier conditions (Martin, 1990;Moore et al., 2013). We apply the regridded LGM dust fields of Mahowald et al. (2006) and denote this change LGM dust.

15
In the sensitivity study, for each of the three C/P parametrisations, we first change one forcing at a time (see Table 1). We run This allows us to test the sensitivity to deep ocean retention of organic carbon. 20 We then run simulations where we combine several changes in forcing. We get a colder climate simulation (LGM phy) by combining LGM rf and LGM alb. In simulation Acomb, we combine LGM rf , LGM alb, LGM dust, and RLS × 1.25 (see Table 1). Kwon et al. (2009) show that small changes in remineralisation depth can cause substantial changes in pCO atm 2 .
With the RLS deepening of 25%, we keep the corresponding changes in pCO atm 2 from exceeding the ∼ 20 − 30 ppm obtained in other studies (Matsumoto, 2007;Menviel et al., 2012). We finally run a glacial-like simulation GLcomb (see Section 3.3), 25 which is similar to Acomb but also includes the change in wind stress W SN × 0.5. The achieved GLcomb model state does not represent a full glacial maximum state, but is more glacial-like compared to the control state; it has a colder climate (see Table 2), reduced deep ocean ventilation and more carbon retention in the deep ocean.

Observations
For comparison and validation of model results, we use records of ocean state variables from observations of modern data and We use model-data comparison of benthic δ 13 C to assess the statistical similarity (correlation) between both the model control state and glacial-like state (see Section 2.3) to benthic δ 13 C data representing the Late Holocene (0-6 ka, HOL) and Last Glacial Maximum (19-23 ka, LGM), respectively Peterson et al. (2014). Locations of core sites can be seen in Fig Peterson et al. (2014)). Note that LGM benthic δ 13 C is more 13 C-depleted than the Holocene due to the addition of 13 C-depleted terrestrial carbon to the glacial ocean (Shackleton, 1977;Curry et al., 1988;Duplessy et al., 1988), which is not simulated in our model experiments. Therefore, to compare our glacial-like simulations (GLcomb) to LGM observations, we subtract a Holocene-LGM global average difference of 0.32 ‰ (Gebbie et al., 2015) from the GLcomb experiments. Gebbie et al. (2015) state that the wide range of error for the estimate of glacial-to-modern change in benthic δ 13 C of of 0.32 ±0.20‰ 10 suffers from a lack of observations in all ocean basins but the Atlantic. Therefore, we place more emphasis on the results of the model-data comparison in the Atlantic than in the Indo-Pacific sector.

Nutrient utilisation efficiency
The extent to which biology succeeds to use the available nutrients can be determined by calculating the nutrient utilisation efficiency P * (Ito and Follows, 2005;Ödalen et al., 2018), which is the fraction of remineralised (P reg ) to total (P tot ) nutrients (in this case, PO 4 ) in the ocean. Overlines denote global averages. Remineralised nutrients have been transported from the surface to the interior ocean by the biological pump, and P rem is given by 20 Here, P pre is preformed P O 4 -the concentration of P O 4 that was present in the water parcel as it sank, thus the fraction that was not used by biology in the surface ocean. In cGENIE, the concentration of preformed tracers is set in the surface ocean and then passively advected through the ocean interior (Ödalen et al., 2018). The biological pump also captures carbon, and a similar relationship can be used for concentrations of DIC, where 25 DIC rem is used to compute the ocean storage of remineralised acidic carbon (AC rem , see Appendix A ). AC rem is biological carbon that entered the ocean in the form of CO 2 in soft tissue (as opposed to carbonates in hard tissue), measured independent of oxygen consumption and/or remineralised phosphate.
In cGENIE, P pre and DIC pre are modelled as passive tracers (Ödalen et al. (2018)). Hence, we can use the model output for P pre in Eq. 3 to compute P * (Eq. 2).
In a model with fixed Redfield ratio, P * determines the effect of the biological pump on pCO atm 2 . For example, it has been found that a higher P * in the initial state gives a lower potential for drawdown of pCO atm 2 in response to similar perturbations (Marinov et al., 2008;Ödalen et al., 2018). However, with variable stoichiometry this is no longer true, since the amount of 5 carbon retained in the deep ocean is not necessarily proportional to P * .

Ocean temperature and circulation
As the three control states, Ctrl RED , Ctrl GAM and Ctrl 121 are driven by the same physical forcings and have the same 10 pCO atm 2 , they have the same ocean circulation pattern ( Fig. 2 a, c, e and Table 2, Table S.2) and climate (exemplified by global ocean average temperature (T oce ) in Table 2 and Table S.2). The surface ocean nutrient fields are fairly similar, with small differences due to the different C/P parametrisations (compare Fig. 3   The agreement with observations is better for Ctrl GAM than for Ctrl RED (Fig. S.2).

Surface nutrient distribution and C/P ratios
In Fig. 3 we see that surface P O 4 fields (left hand column) and the corresponding fields of surface C/P ratios (as given by Eq. 1, right hand column) of Ctrl GAM (panels a, b) and of observations (panels c, d) are similar in their pattern as well as in the magnitudes of values. Note that high concentrations of PO 4 correspond to low C/P ratios, and vice versa. The highest 25 observed P O 4 concentrations in the Northern and equatorial Pacific are not fully reproduced by the model, but the pattern is well reproduced. In the surface C/P field of Ctrl GAM (Fig. 3 b) we see the signature of very high nutrient concentrations in the Southern Ocean (> 1µmolL −1 , Fig. 3 b) as a band of low ratios, with the most extreme values near the Antarctic continent, as seen in observations (Fig. 3 d).
The nutrient utilisation efficiency P * (Eq. 2) in the three control states differs by a few percent; 0.43, 0.46 and 0.42 in Ctrl RED , Ctrl GAM and Ctrl 121 respectively ( Table 2). The fraction of DIC rem in DIC tot (see Section 2.5, Eq. 4, Table   S.1) is 0.065, 0.077 and 0.072 in Ctrl RED , Ctrl GAM and Ctrl 121 respectively.

Ocean dissolved O 2
The most apparent difference between Ctrl RED and Ctrl GAM is in deep ocean oxygen concentrations, where the global  Table 2.

Ocean δ 13 C
By comparing Ctrl RED δ 13 C and Holocene (0-6 ka, HOL) benthic δ 13 C values, we estimate a global model-data correlation of 0.78 (Table S.3). The modern-day Atlantic Ocean has a distinctive spatial δ 13 C pattern ( Atlantic. While the model produces a weaker gradient than the observed HOL Atlantic Ocean (corr. 0.50, Table S.3) the model correlates well with Eastern Atlantic δ 13 C records ( Fig. S.1). For the Indo-Pacific, the weaker benthic δ 13 C gradient is well represented by the model (Fig. 5 e). This pattern emerges mainly due to 13 C-depleted, biologically sourced carbon that is accumulated in the weak circulation region of the interior North Pacific (Matsumoto et al., 2002). However, Indo-Pacific δ 13 C 20 values of Ctrl RED are overall lower than the HOL observations. The overall model-data correlation for the Indo-Pacific is 0.39 (Table S. a, e) and model-data correlations with HOL observations (Table S.3) are similar between the model versions, with somewhat lower correlations for GAM .

Sensitivity experiments 25
The applied changes listed in Table 1 cause changes in ocean characteristics such as overturning circulation, temperature, surface nutrient distributions and biological productivity, which result in changed pCO atm 2 . The resulting steady state global average values for temperature (T oce ), dissolved oxygen (O 2 ), and nutrient utilisation efficiency P * , as well as the maximum and minimum of the Atlantic meridional overturning streamfunction, are listed in Table 2.

Radiative forcing and albedo
In the simulations where radiative forcing and albedo are changed to represent LGM conditions (LGM rf + LGM alb = LGM phy), the reductions in pCO atm 2 are similar in the RED and the GAM model versions. In LGM phy, the resulting pCO atm 2 is 245.4 and 244.9 ppm respectively, thus a reduction of 33 ppm compared to the Ctrl 278 ppm (Fig. 6). Here, variable C/P does not impact the results, because changes in the surface nutrient distribution (Fig. 7 a), and the associated 5 changes in C/P (Fig. 8 a), are limited to very high latitudes where productivity is already low in the control state, due to low temperatures and a lack of light and iron. The drawdown of pCO atm 2 can mainly be attributed to the increase in solubility carbon (C sat,T ) due to ocean cooling, and to an increase in sea ice, which prevents air-sea gas exchange and therefore causes an increase in disequilibrium carbon (C dis ) (Ödalen et al., 2018). Ocean cooling amounts to 2.1°C in LGM phy compared to Ctrl (T oce in Table 2). In cGENIE, the increase in C sat,T associated with ocean cooling corresponds to ∼ 7±1. ppm. As we do not change salinity, we are simultaneously likely to overestimate the increase in solubility between Ctrl and a glacial-like state, by ∼ 6 ppm (Kohfeld and Ridgwell, 2009). This effect is consistent for any choice of C/P parametrisation, and is therefore not explored further.

Reduced wind forcing
When the peak of Southern Ocean (henceforth SO) winds is reduced, the strength of the overturning circulation of AABW 20 decreases (see difference in Southern Hemisphere overturning streamfunction between Fig. 2 a and b). Thus, given that the volume of AABW does not change, its residence time increases. This also means that the upwelling nutrient-rich water in the SO stays a longer time near the surface and loses more nutrients before being subducted. This decreases the SO concentration of preformed phosphate in W N S ×0.5 RED compared to Ctrl, as seen in Fig. 7 b, and increases the nutrient utilisation efficiency P * (Table 2, Fig. 9). This leads to a drawdown of pCO atm 2 of 12.9 ppm compared to Ctrl RED (Fig. 6). As the nutrient 25 concentration in the SO decreases ( Fig. 7 b), the flexible C/P ratio (Fig. 8 b) leads to an increased carbon capture efficiency in GAM compared to RED (see GLcomb−Ctrl of biologically sourced carbon (AC rem ) in Fig. 9), which is partly compensated by a reverse effect in the Pacific equatorial region. Consequently, in W SN × 0.5 GAM , we get a reduction of pCO atm 2 of 16.3 ppm compared to Ctrl GAM (Fig. 6). Hence, for halved peak wind stress at ±50°N, the flexible stoichiometry increases the drawdown by ∼ 26 %.

Remineralisation length scale
When the remineralisation length scale (RLS) increases, the biological material reaches deeper before it is remineralised, and it takes longer for it to be returned to the surface. Therefore, more of the biologically sourced carbon (AC rem ) and nutrients are present in the deep ocean at any given time, leading to an increase in P * (  Fig. 7 c), which, through the small resulting changes in C/P (Fig. 8

Dust forcing
The simulations with LGM dust forcing (LGM dust, Table 1) show the largest difference in pCO atm 2 between the RED and the GAM . In LGM dust RED , pCO atm 2 decreases by 16 ppm compared to Ctrl RED , whereas LGM dust GAM sees a reduction of 21 ppm compared to Ctrl GAM (Fig. 6). The drawdown is thus ∼ 30 % larger with variable stoichiometry. About a third (∼ 10 % of ∼ 30 %) of the increase in drawdown can be explained by a change in average C/P composition of the 20 organic material that is exported out of the surface ocean (see Section 4.4).
As anticipated, the iron added by the dust forcing allows more efficient usage of P in the HNLC-regions, which increases the ocean storage of biologically sourced carbon P rem (thus, P * increases, Table 2, Fig. 9). This reduces the surface nutrient concentrations in these areas (Fig. 7 d). In the GAM model version, this is followed by increased C/P ratios in these areas ( Fig.   8 d), resulting in a lower pCO atm 2 in LGM dust GAM than in LGM dust RED . The largest anomalies in P O 4 concentrations, 25 and consequently in C/P, are observed in subantarctic zone of the Southern Ocean, particularly in the Atlantic and Indian sectors. This subantarctic increase in biological efficiency is consistent with radionuclide proxy data ( 10 Be, 230 Th, 231 Pa) from the LGM (e.g., Kumar et al., 1995;Kohfeld et al., 2005).

Combined experiments
We show the results of two different combined simulations; Acomb and GLcomb. GLcomb is the "glacial-like" simulation, 30 which combines all the sensitivity experiments (Table 1). Acomb omits the reduction in wind stress.

Ocean temperature and circulation
In the glacial-like simulations, GLcomb RED and GLcomb GAM , the global average ocean temperature (T oce ) is 1.7°C lower than in the respective control states ( and 82 % of the glacial-interglacial difference in T oce , respectively. As anticipated, our combined forcings do not induce a full glacial maximum state, but a state with glacial-like climate conditions. Ocean overturning circulation weakens in GLcomb compared to Ctrl (Fig. 2, Table 2), mainly as a result of the wind stress 10 reduction. For example, the AMOC (here measured as the maximum of the Atlantic overturning streamfunction) reduces in strength by ∼ 15 %, from 14 Sv to 12 Sv ( Table 2). The global meridional overturning stream function reveals that the SO overturning cell sees a reduction in transport ( Fig. 2 d), which is associated with weaker upwelling and thus longer residence time for AABW, as hypothesised for the glacial ocean (e.g., Menviel et al., 2017;Skinner et al., 2017). In Acomb, where the wind stress is kept at modern values, the ocean overturning circulation remains similar to the control state (Table 2).

Surface nutrient distribution and C/P ratios
In the surface nutrient anomalies (GLcomb − Ctrl, shown for GAM in Fig. 7 f), we see the strongest response in the Southern Ocean, with different effects south and north of the so-called biogeochemical divide described by Marinov et al. (2006). Marinov et al. (2006) show that the air-sea balance of CO 2 is dominated by processes in the waters close to Antarctica, whereas global export production is instead controlled by the biological pump and circulation in the Subantarctic region. The 20 border between these two regimes is referred to as the biogeochemical divide. South of the biogeochemical divide, close to the Antarctic continent, we see an increase in GLcomb nutrient concentrations compared to Ctrl (Fig. 7 f), which coincides with an increase in sea-ice in this area (not shown). Colder conditions due to changed albedo and radiative forcing, with more sea ice than in the control state, cause a reduction in biological production, leaving more unused P in the surface layer (Fig. 7 a).
North of the biogeochemical divide, increased aeolian dust flux increases the productivity of the biology, which reduces P in 25 the surface compared to the control (Fig. 7 d). In combination with circulation changes, resulting from the reduced SO wind stress (Fig. 7 b), and deeper remineralisation (Fig. 7 c), P concentrations in the Subantarctic region are strongly reduced (Fig.   7 f). In the North and Equatorial Pacific, there is also a reduction of P (Fig. 7 f), mainly due to to the increased dust flux (Fig.   7 d). There is an increase in P seen in the Arctic, again coincident with an increase in sea ice in the same area. As a result of these changes, we see strong positive anomalies in C/P in the HNLC-regions, and negative anomalies in the highest latitude 30 bands (Fig. 8 f). The organic matter that is exported out of the upper layer (henceforth referred to as export production) in GLcomb GAM has a global average C/P ratio of 134/1.
Despite colder conditions, which allow for more dissolution of O 2 , the reduction in O 2 is evident in the GLcomb-simulations (Fig. 4). This mirrors the increase of AC rem (Fig. 9). The O 2 reduction is about 50 % larger in GAM compared to RED. As the initial state Ctrl GAM is already lower in oxygen than Ctrl RED (144µmolkg −1 compared to 166µmolkg −1 ), and variable stoichiometry allows for additional ocean storage of organic carbon, the end state O 2 is drastically lower in GLcomb GAM

Ocean δ 13 C
In GLcomb RED (contours in Fig. 5 b,d), the Atlantic North-South gradient in δ 13 C is stronger than in Ctrl RED (contours in Fig. 5 a,c). This strong gradient is not observed in the Holocene Atlantic δ 13 C-data (dots in Fig. 5 a,b), but is prominent in the LGM time slice (dots in Fig. 5 c,d). The LGM observations are well reproduced in GLcomb red (corr. 0.62, Table S.3), 10 especially in the East Atlantic (corr. 0.77). When we correct GLcomb RED for the absence of injected terrestrial carbon, we see clear similarities with LGM observations (Fig. 5 d), though the southernmost cores still indicate more 13 C-depleted conditions than the model. In the Indo-Pacific, GLcomb RED is too 13 C-depleted compared to LGM observations, particularly with the correction for the absence of a terrestrial signal (Fig. 5 h), and the model-data correlation is poor (0.05). For the Indo-Pacific, the model-data correlation for Holocene data is similar between GLcomb RED and Ctrl RED (0.24 and 0.39 respectively, 15 Table S.3). This suggests that the poor correlation with LGM data is simply due to our changes in forcings being insufficient to achieve the required rearrangements in Indo-Pacific circulation patterns. In GLcomb, similarily as in Ctrl, there is very little difference between the RED and GAM model versions in terms of δ 13 C patterns (compare Fig. 5   In the combined experiments Acomb and GLcomb, pCO atm 2 decreases strongly compared to the control state (reduction of between −56 to −80 ppm, Fig. 6). This is partly a result of colder conditions (Table 2), which lead to increased solubility for CO 2 in sea water and an expanded sea-ice cover ( Figure S.4) which restricts air-sea gas exchange. The latter slows down the equilibration of the CO 2 -rich upwelling water with the atmosphere before it is subducted into the deep ocean, and thus 25 causes increased C dis . In Acomb and GLcomb, changes in biological production (see Section 3.3.2) and storage of biologically sourced carbon (Fig. 9) also contribute strongly to the reduced pCO atm 2 .
The combined experiments show a striking difference in pCO atm 2 between the model versions RED and GAM ; in GLcomb RED and GLcomb GAM , we achieve drawdown of pCO atm 2 of −64 and −80 ppm, respectively, from the 278 ppm of the control states (Fig. 6). This corresponds to an increase in ocean carbon storage of 139 PgC and 173 PgC, respectively (Table S.1). The drawdown is thus 25% larger in GAM than in RED. In Acomb, where the perturbation in wind stress is omitted, drawdown is smaller than in GLcomb and the difference between model versions is less pronounced, but there is still a 14 % difference between Acomb RED and Acomb GAM (−56 and −63ppm, respectively (Fig. 6).

Accounting for variable C/P in ocean carbon cycle models
The representation of ocean biology in General Circulation Models (GCMs) tends to be over-simplified (Le Quéré et al., 2005a) 5 and the development of the models is often held back by constraints imposed by maintaining the computational efficiency of the model. The Galbraith and Martiny (2015) empirical model is simple and based on nutrient variables that are already present in biogeochemical models (nitrate and phosphate). By implementing the GAM parametrisation, or possibly a power law as that described by (Tanioka and Matsumoto, 2017), an additional facet of the complexity of ocean biology can be implicitly The approach taken by both Galbraith and Martiny (2015) and Tanioka and Matsumoto (2017) is to adapt a function to all the available species-independent observations. This means that they account for the adaptation of plankton to the surrounding conditions, both in terms of species composition and individuals being more frugal in low nutrient conditions. This is one of the main advantages of such an approach, as it can be applied in a model without different plankton functional types, which 20 is what we use here. Ganopolski and Brovkin (2017) appear to succeed with full glacial-interglacial CO 2 cyles in CLIMBER-2, which does not have flexible stoichiometry for primary producers, but which has a temperature limitation on growth and explicit phyto-and zooplankton with different C uptake rates. This combination may perhaps achieve a similar response in carbon export in their simulations, when moving into a colder climate, as the flexible stoichiometry does in our simulations.
The next step to approach a more realistic modelling of the biological pump would be to include a representation of preferential 25 remineralisation of nutrients (e.g., Kolowith et al., 2001;Letscher and Moore, 2015), but this goes beyond the scope of the present study.
One drawback of the Galbraith and Martiny (2015) approach is that it assumes that the C/P ratio continues to increase continuously with increasing [P O 4 ]. As such, in a high surface [P O 4 ] region like the Southern Ocean, GAM does not represent the effects on C/P of temperature and light, and the associated non-linear effects, that could be of importance here (e.g., Yvon- in GAM compared to RED for both increased and decreased RLS. While increased RLS would be an effect of ocean 5 cooling, and thus of interest for glacial studies, reduced RLS would be a consequence of ocean warming (Matsumoto, 2007).
Matsumoto (2007)  (e.g., Sarmiento and Toggweiler, 1984;Martin, 1990;Archer et al., 2000;Sigman and Boyle, 2000). However, biological production depends on water temperature (Eppley, 1972) and decreases in cold conditions. This temperature effect is parametrised in cGENIE as a local temperature-dependent uptake rate modifier proportional to e (T /15.9 • C) , and we see overall reduced productivity in our cold climate simulations (see e.g. LGM phy , Table S.1). As the climate cools, the temperature 15 effect leads to a decrease in biological productivity and a subsequent decrease in P * (see LGM phy in Fig. 9). If productivity decreases, other mechanisms are needed to offset this decrease, if the total ocean storage of organic carbon is to increase.
Mechanisms that can contribute to increased deep ocean carbon retention are for example reduced SO overturning circulation is evident that variable stoichiometry can be another contributing factor. In our simulations, global average export production decreases, both in terms of P OP (particulate organic phosporus) and P OC (particulate organic carbon), by 15 % between Ctrl RED and GLcomb RED because of the colder climate. However, the variable stoichiometry in GAM partly offsets this decrease in biological carbon capture -while the export production in terms of P OP decreases by 17 %, the corresponding decrease in P OC is only 10 %. Thus, in GLcomb GAM we achieve an increase in ocean inventory of remineralised carbon, 25 which exceeds the one of GLcomb RED (Fig. 9).
In summary, we suggest this is a result of changes in surface P fields (see Fig. 7) Due to the competing effects between decreased export production and increased retention that are introduced by the different forcing components, the net change in P rem , and thus in P * , is very small when moving from the control state (Ctrl) to the glacial-like state (GLcomb) (Fig. 9). In the fixed stoichiometry case (RED), there is a small net increase in P * of 0.020 in 30 GLcomb compared to Ctrl, which is linearly related to a small increase in storage of AC rem in the ocean. In the case with variable stoichiometry (GAM ), there is instead a very small decrease in P * of 0.003. With a linear response, we would thus expect a decrease in storage of AC rem as well. Instead, we see a reasonably large increase in global average AC rem .
The reduced pCO atm 2 in GLcomb GAM , compared to GLcomb RED , can be explained by the non-linearity introduced by the local variability in C/P. When changes in ocean circulation, remineralisation depth and dust deposition cause the local nutrient availability in the surface waters to change, this affects the elemental composition of the exported organic material. In GLcomb, there is a reduction of the surface layer P concentration compared to Ctrl in some key (HNLC) regions. In GAM , this decrease in surface P results in an increase in C/P. This further strengthens the biological pump in these key regions, 5 resulting in a non-linear relationship between storage of remineralised phosphate and biologically sourced carbon (Fig. 9).
In Ctrl GAM , the average elemental C/P composition is 121/1. In GLcomb GAM , this average is 134/1. This means that even though the same amount of P is exported to the deep ocean, the organic molecules carry more carbon, which is released in the deep ocean during remineralisation. In Ctrl GAM , the global average concentration of P rem is 1.16 µmolkg −1 (c.f. 1.17 µmolkg −1 in GLcomb GAM ). By increasing the average C/P composition of 1.16 µmolkg −1 organic molecules from 121 to 10 134 (i.e. by 13 units), this causes an increase in C sof t by ∼ 15 µmolkg −1 , which corresponds to the observed increase in AC rem . From this, we conclude that, in a system where stoichiometry is variable on a local scale, ocean storage of biologically sourced CO 2 can change while the amount of remineralised nutrients remains constant.
Note also that Ctrl GAM has a larger inventory of DIC as well as C sof t compared to Ctrl RED . Ödalen et al. (2018)

Implications of flexible C/P for deep ocean oxygen
Studies have shown that deep ocean O 2 concentrations were lower during the LGM than in the Holocene (e.g., Bradtmiller et al., 2010;Jaccard and Galbraith, 2012;Galbraith and Jaccard, 2015). Here, we discuss implications of using flexible C/P for the model's ability to reproduce ocean oxygen patterns and concentrations in the different time periods.
In the Atlantic, Ctrl GAM (Fig. 4 e) reproduces the observed extent of the Oxygen Minimum Zone (OMZ) (Fig. 4 a) better 25 than Ctrl RED (Fig. 4 c) does. In the Pacific Ocean, the O 2 gradient in the observations (Fig. 4 b) is more gradual compared to that of the control states (Fig. d, f), but the core of the OMZ is well reproduced by the model. The forcings applied to GLcomb are not sufficient to reproduce a full glacial state. Still, we do get a vertical expansion of the OMZ in GLcomb RED (Fig. 4 h) compared to Ctrl RED (Fig. 4 c), in agreement with the findings of Hoogakker et al (2018). In GLcomb GAM , oxygen depletion is too extensive (see below), but the tendency of vertical expansion compared to the control state is present here as well. 30 In our glacial-like states, we see a significant reduction of O 2 compared to the control state in both RED and GAM (Fig. 4 d and e). This response is expected when we apply increased dust deposition and deeper remineralisation (Galbraith and Jaccard, 2015), and the direction of the overall response of deep ocean O 2 to the glacial-like forcings is in line with observations. The reduction in O 2 is stronger in GAM , due to the larger storage of respired carbon in this model version.
Finally, it should be noted that Ctrl GAM (Fig. 4 c) displays deeper oxygen minima in the oxygen minimum zones of both the Atlantic equatorial region and the North Pacific compared to what is seen in Ctrl RED (Fig. 4 b), and in observations.
In Ctrl GAM , a large part of the interior North Pacific is anoxic (Fig. 4 c), while observations (Fig. 4 a) indicate very low oxygen levels, but not anoxia. As illustrated by simulations in Galbraith and de Lavergne (2018), variable stoichiometry in itself is not always sufficient to achieve widespread deep ocean de-oxygenation in a model under glacial-like climate change.

5
Among other factors, model ocean oxygen conditions are also dependent on deep water formation characteristics of the model (Galbraith and de Lavergne, 2018). The deep water formation characteristics of a model affects the amount of time available for remineralisation and, consequently, the oxygen consumption. In addition, due to a lack of resolution deep water formation in climate models generally happens as open water convection, rather than as dense plumes along slopes (Heuzé et al., 2013).
This may cause too much oxygen to be entrained into the deep ocean Galbraith and de Lavergne (2018). In cGENIE, this 10 effect is small enough not to cancel the increased O 2 consumption caused by the higher average C/P in Ctrl GAM compared to In summary, the GAM version of the model reproduces quantitatively the observed O 2 patterns of both the LGM and the Holocene, but with too low concentrations. If this parametrisation of variable stoichiometry is to be used in cGENIE in future studies, we suggest some re-tuning, for example by reducing the global average concentration of nutrients, to improve the 15 representation of observed modern-day ocean O 2 concentrations in the GAM control state.

Effect of modified but fixed C/P
Part of the observed difference in pCO atm 2 between GLcomb RED and GLcomb GAM results from a difference in global average C/P in the control states (Ctrl RED and Ctrl GAM ). In Ctrl GAM , the average C/P in the export production is close to 121/1, instead of 106/1 as in Ctrl RED . We illustrate the consequences of this difference by running parallell simulations with 20 fixed C/P of 121/1 (model version 121).
In a simulation with reduced wind stress (W N S × 0.5 121 ), pCO atm 2 is reduced by −14.5 ppm compared to Ctrl RED (Table S. that are due to changed average C/P. As shown above, the simulations with 121 indicate that, depending on the change in forcing, between 1/3 and 2/3 of the difference in drawdown between RED and GAM is due to the difference in average C/P between the control states (Fig. 6, Table S.2). From this, we conclude that the effects of the perturbations do not add linearly.
As outlined in Section 4.3, an increase in C/P reduces deep ocean O 2 concentrations, through an increase in regenerated carbon. In this respect, the 121 experiments fall between the corresponding RED and GAM experiments. For example, the global ocean average O 2 concentration, O 2 , in Ctrl 121 is 152 µmolkg −1 , which is lower than Ctrl RED (166 µmolkg −1 ), but higher than Ctrl GAM (144 µmolkg −1 ). Similarly, GLcomb 121 has a lower O 2 than GLcomb RED (96 compared to 122 µmolkg −1 ), but higher than GLcomb GAM (74 µmolkg −1 ).

5
The observed effects of modified average C/P could have implications for model intercomparison projects, if they compare results from models that use different versions of fixed stoichiometry (for example, Anderson and Sarmiento (1994) or Takahashi et al. (1985) stoichiometries, compared to Redfield (1963)). The problem with different stoichiometry assumptions in models is extensively discussed by Paulmier et al. (2009). Our study shows a direct consequence of such differences, with different model response to the same perturbation. Proxy records of benthic δ 13 C indicate a change in ocean δ 13 C across the deglaciation. The whole ocean deglacial change has been estimated to ∼ 0.35 ‰ (Peterson et al., 2014;Peterson and Lisiecki, 2018), and the surface-to-deep gradient weakened (shown in numerous studies, see e.g. Curry et al., 1988;Duplessy et al., 1988;Curry and Oppo, 2005;Herguera et al., 2010;Peterson and Lisiecki, 2018, , and references therein). Here we compare model ocean δ 13 C of our simulations to the 15 benthic δ 13 C records, to see how well the simulations capture the observed patterns, and if there is a difference in model-data correlation between RED and GAM .
The model-data comparison in δ 13 C (Section 3.1.4) suggests that the Ctrl simulations are overall well correlated with Holocene benthic δ 13 C data (HOL in Table S.3). For the Atlantic, the correlation of the Ctrl simulations is higher with the LGM benthic δ 13 C (LGM in Table S.3) than with HOL. However, the standard deviations (STDs) suggest that the Atlantic 20 North-South gradient is not as strong as in an LGM ocean state and thus more similar to the modern ocean.
When we apply our combined forcings in the GLcomb simulations, we achieve a stronger δ 13 C gradient in the Atlantic, allowing for a closer match with LGM data in terms of STDs. A stronger gradient in δ 13 C and more depleted suggest weaker ventilation of the deep ocean. The poor correlation in the Indo-Pacific, which may be partly due to sparse mid-ocean observations for almost 70% of the ocean volume, makes the global statistics for LGM observations difficult to interpret. The forcings 25 applied to GLcomb are factors that are likely to be important for the glacial ocean circulation and biogeochemistry. However, these forcings are not sufficient to reproduce a full glacial state (i.e. the use of the term glacial-like simulations, rather than LGM simulation). Other forcings that have shown to be important for modelling of glacial δ 13 C are, for example, brine rejection (Bouttes et al., 2010(Bouttes et al., , 2011, and freshwater forcing (Schmittner et al., 2002;Hewitt et al., 2006;Bouttes et al., 2012). The fact that some important forcings are missing is likely the main cause for the model-data discrepancy, and the reason for why 30 we do not achieve a glacial Pacific Ocean circulation consistent with observed δ 13 C patterns.
Each of the two observational datasets (HOL and LGM) display similar correlations across the two model simulations. The correlation of the HOL δ 13 C records with Ctrl RED , Ctrl GAM , GLcomb RED , and GLcomb GAM , is in all cases between 0.76-0.78. Meanwhile, the correlation of LGM δ 13 C records with the same four simulations is in all cases between 0.55-0.58.
As our GLcomb-simulations still correlate so well with the HOL dataset, this suggests the applied forcings have not caused these simulations to be clearly different from Ctrl in terms of water mass distribution. For the same reason, the correlation with LGM δ 13 C records does not significantly improve from Ctrl to GLcomb. The water mass distribution in cGENIE is strongly constrained by the resolution of the model, especially in the vertical. Changes in temperature and salinity that should cause changes in water mass volume may not be sufficient to allow a water mass to extend to the next vertical level of the 5 model. As a consequence, while the gradient between water masses may become more or less pronounced, the interface of water masses may still remain at the same depth. The applied changes affect the chemical and biological conditions for ocean carbon storage, such as CO 2 solubility (temperature dependent) and nutrient availability, more than the physical conditions, such as water mass volume and turnover time. To achieve a full glacial state with cGENIE, with a more glacial-like water mass distribution, additional physical forcings (see above) are likely to be required.

10
How δ 13 C is represented in cGENIE is detailed in Appendix B. The very small differences between RED and GAM suggests that using variable stoichiometry does not impact our ability to represent the δ 13 C patterns seen in observational data.
However, we note that some retuning of the modern control state is recommended before cGENIE with variable stoichiometry is used in other studies.

15
In this paper, we examine the potential role of variable stoichiometry in biological production for glacial ocean CO 2 storage.
We show that flexible C/P composition of organic matter allows a stronger response of pCO atm 2 to glacial-like changes in climate, remineralisation length scale and aeolian dust flux. We conclude that variable stoichiometry may be important for glacial ocean CO 2 storage and for achieving the full extent of drawdown of atmospheric CO 2 in model simulations. In the experiment GLcomb RED , with glacial-like climate and Redfield stoichiometry (Redfield, 1963), ocean carbon storage increases by 139 20 PgC and atmospheric CO 2 decreases by 64 ppm. In GLcomb GAM , with glacial-like climate changes and variable stoichiometry, the corresponding numbers are 173 PgC and 80 ppm. Hence, the drawdown of atmospheric CO 2 increases by 25 % when C/P is variable.
About half of the increased drawdown of CO 2 results from different global average C/P in the export production. In addition, flexible stoichiometry allows increased carbon capture through the biological pump, while maintaining or even decreasing the 25 fraction of remineralised to total nutrients in the deep ocean. With fixed stoichiometry, an increase in remineralised carbon is inevitably tied to a corresponding increase in remineralised nutrients.
We apply variable C/P parametrised as a simple function of the surface water concentration of P O 4 , as suggested by Galbraith and Martiny (2015). Tanioka and Matsumoto (2017) suggest that it is unrealistic for C/P to continue to increase indefinitely with increased [P O 4 ], and therefore suggest a more complex power law function, which takes into account saturation of 30 the C/P ratio at high concentrations of P O 4 . However, we found that saturation of the C/P ratio at concentrations higher than the observational upper bound of 1.7 µM causes no noticable impact on our results.
The representation of flexible stoichiometry used in this study (Galbraith and Martiny, 2015) can be used without large increases in computational cost. It makes it possible to take into account, to first approximation, the complex biological changes that occur in the ocean during long-timescale climate change scenarios (see e.g. McInerney and Wing, 2011). We show here that, for glacial-interglacial cycles, this complexity contributes to changes in atmospheric CO 2 through flexible C/P ratios.
Flexible C/P also has the potential to be an additional positive feedback of ocean warming on pCO atm 2 in future climate.
Configuration files for the specific experiments presented in the paper can be found in the directory Appendix A: Regenerated acidic carbon 15 Carbon enters the ocean mainly in the form of CO 2 and dissolved carbonate. Despite this, the major fraction of carbon in the ocean resides as bicarbonate ions. The source related state variables acidic and basic carbon (AC and BC, respectively) allow us to separate the ocean DIC inventory into the sources of CO 2 (AC) and dissolved carbonate (BC). The concept of the sourced related state variables was first described by Walin et al. (2014).
AC and BC are defined from DIC and alkalinity (ALK) as Total AC and BC include all ocean sources of carbon, including (but not limited to) river runoff, air-sea gas exchange, and the biological pump.
To isolate the CO 2 that was supplied to the ocean via biological soft tissue, we make use of the separation of DIC and ALK 25 into their preformed and remineralised fractions (see Section 2.5). Thus, we then compute the remineralised acidic carbon (AC rem ) as Appendix B: δ 13 C in cGENIE cGENIE represents 13 C as an explicit tracer (separate from and in addition to bulk carbon) in the model, tracking its concentration in all the same gaseous, dissolved, and solid forms that carbon exists in, reporting δ 13 C in ‰ relative to the standard VPDB. The current scheme is based on that described in Ridgwell (2001) and updated as described in Ridgwell et al. (2007) and is evaluated for the modern ocean (alongside simulated ∆ 14 C) in cGENIE, in Turner and Ridgwell (2016).

5
In the aqueous phase, the isotopic partitioning of carbon between CO 2 (aq), HCO − 3 , and CO 2− 3 is resolved and follows Zeebe and Wolf-Gladrow (2001) (their Section 3.2). The empirical fractionation factors used are from Zhang et al. (1995). The airsea fractionation scheme follows that of Marchal et al. (1998) with the individual fractionation factors again taken from Zhang et al. (1995).
For the isotopic composition of organic carbon (δ 13 C P OC ), the model of Rau et al. (1996) is adapted, assuming that the 10 isotopic signature of exported POC reflects that of phytoplankton biomass. Following Ridgwell (2001), the full equation of Rau et al. (1996) is simplified to: where [CO 2 (aq)] is the ambient concentration of aqueous CO 2 and δ 13 C aq is its isotopic composition. K Q is a temperatureonly dependent approximation of the full cell-dependent size and growth rate parameterization in the Rau et al. (1996) model 15 (see Ridgwell, 2001). We take an intermediate value for the enzymatic isotope fractionation factor associated with intracellular C fixation ( f ) of -25‰ following Rau et al. (1996Rau et al. ( , 1997, and assume a temperature-invariant value for d of 0.7 ‰. The result of applying this scheme in cGENIE, is a zonal mean profile characterized by δ 13 C P OC of -22 to -21 ‰ in the tropics, declining with increasing latitude to reach -28 to -30 ‰ in the Southern Ocean. This latitudinal pattern is comparable to measurements made on suspended particulate organic matter as discussed in Ridgwell (2001). 20 For 13 C fractionation into biogenic carbonates at the ocean surface (e.g. foraminiferal tests, and coccolithophorid coccoliths), cGENIE follows Mook (1986) and employs a simple temperature-dependent fractionation between the δ 13 C of aqueous HCO − 3 and calcite.  260-263, 1994. Archer, D., Winguth, A., Lea, D., and Mahowald, N.: What caused the glacial/interglacial atmospheric pCO2 cycles?, Reviews of Geo- Menviel, L., Joos, F., and Ritz, S.: Simulating atmospheric CO2, 13 C and the marine carbon cycle during the Last Glacial-Interglacial cycle: possible role for a deepening of the mean remineralization depth and an increase in the oceanic nutrient inventory, Quaternary Science   streamfunction (1 Sv = 10 6 m 3 s −1 ) of CtrlGAM (panels a, c, e) and GLcombGAM (panels b, d, f). Note that the eddy-induced transport of tracers is taken into account through a skew-diffusive flux (Griffies, 1998) that is present in the velocity fields used to compute the eulerian stream function.

HOL Atlantic
LGM Atlantic

HOL Pacific
LGM Pacific LGM Atlantic

HOL Pacific
LGM Pacific LGM A

HOL Pacific
LGM Pacific

HOL Pacific
LGM Pacific

HOL Pacific
LGM Pacific LGM A

HOL Pacific
LGM Pacific LGM Atlantic

HOL Pacific
LGM Pacific  LGM Atlantic

HOL Pacific
LGM Pacific The columns represent the model simulations (CtrlRED or CtrlGAM ), while each row represents one of the proxy record time slices (HOL or LGM). The left hand column shows CtrlRED (panels a, c, e, g) and the right hand column shows GLcombRED (panels b, d, f, h). The rows show, from top to bottom, a-b) HOL Atlantic, c-d) LGM Atlantic, e-f) HOL Pacific, g-h) LGM Pacific. Note that, before we compare GLcombRED to LGM observations (panels d and h), a constant of 0.32 ‰ is subtracted from the simulated δ 13 C, to account for terrestrial release of δ 13 C-depleted terrestrial carbon which is not modelled. The corresponding comparison for model version GAM is shown in Fig.   S.3.    (f) Figure 8. Surface C/P anomaly, with respect to C/P of CtrlGAM (Fig. 3b) LGMphy LGMrf LGMalb RLSx0.75 RLSx1.25 RLSx1.75 LGMdust WSNx0.5 Acomb GLcomb  is either prescribed to pre-industrial (PI) = 278 ppm, or freely varying with changes in climate and ocean circulation. The radiative forcing is either coupled to the pCO atm 2 of the atmospheric chemistry module of the model, or fixed at a value corresponding pCO atm 2 = 185 ppm. The zonal albedo profile is either representative of modern/PI conditions or of the LGM. The wind stress is either modern/PI, or has an adjusted peak in wind stress at ±50 • N.
The modern/PI remineralisation length scale (RLS) is 590 m. If RLS is changed, it is multiplied by a factor f r.
Dust forcing is either modern/PI or representative of LGM. Each experiment is conducted using model versions with C/P fixed at 106/1 (Redfield, 1963, , denoted RED), C/P variable with surface ocean P O4 concentration (Galbraith and Martiny, 2015, denoted GM 15), LGM (all forcings)

ppm
Acomb "GL" with variable fixed at LGM modern/PI 590 × 1.25 LGM PI wind 185 ppm 1) See Lauderdale et al. (2013) for example of reduced peak wind profile for the Southern Hemisphere.
2) f r = multiplication factor for remineralisation length scale. We test multiplication factors between 0.75 and 1.75, corresponding to a change in RLS between -25 % to +75 %.