Articles | Volume 17, issue 8
Biogeosciences, 17, 2219–2244, 2020
Biogeosciences, 17, 2219–2244, 2020

Research article 22 Apr 2020

Research article | 22 Apr 2020

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

Variable C∕P composition of organic production and its effect on ocean carbon storage in glacial-like model simulations
Malin Ödalen1,2, Jonas Nycander1, Andy Ridgwell3,4, Kevin I. C. Oliver5, Carlye D. Peterson3, and Johan Nilsson1 Malin Ödalen et al.
  • 1Department of Meteorology, Bolin Centre for Climate Research, Stockholm University, 106 91 Stockholm, Sweden
  • 2Department of Geosciences, University of Arizona, Tucson, AZ 85721, USA
  • 3Department of Earth Sciences, University of California–Riverside, Riverside, CA 92521, USA
  • 4School of Geographical Sciences, Bristol University, Bristol BS8 1SS, UK
  • 5National Oceanography Centre, Southampton, University of Southampton, Southampton SO14 3ZH, UK

Correspondence: Malin Ödalen (


During the four most recent glacial maxima, atmospheric CO2 has been lowered by about 90–100 ppm with respect to interglacial concentrations. It is likely that most of the atmospheric CO2 deficit was stored in the ocean. Changes in the biological pump, which are related to the efficiency of the biological carbon uptake in the surface ocean and/or of the export of organic carbon to the deep ocean, have been proposed as a key mechanism for the increased glacial oceanic CO2 storage. The biological pump is strongly constrained by the amount of available surface nutrients. In models, it is generally assumed that the ratio between elemental nutrients, such as phosphorus, and carbon (C∕P ratio) in organic material is fixed according to the classical Redfield ratio. The constant Redfield ratio appears to approximately hold when averaged over basin scales, but observations document highly variable C∕P ratios on regional scales and between species. If the C∕P ratio increases when phosphate availability is scarce, as observations suggest, this has the potential to further increase glacial oceanic CO2 storage in response to changes in surface nutrient distributions. In the present study, we perform a sensitivity study to test how a phosphate-concentration-dependent C∕P ratio influences the oceanic CO2 storage in an Earth system model of intermediate complexity (cGENIE). We carry out simulations of glacial-like changes in albedo, radiative forcing, wind-forced circulation, remineralization depth of organic matter, and mineral dust deposition. Specifically, we compare model versions with the classical constant Redfield ratio and an observationally motivated variable C∕P ratio, in which the carbon uptake increases with decreasing phosphate concentration. While a flexible C∕P ratio does not impact the model's ability to simulate benthic δ13C patterns seen in observational data, our results indicate that, in production of organic matter, flexible C∕P can further increase the oceanic storage of CO2 in glacial model simulations. Past and future changes in the C∕P ratio thus have implications for correctly projecting changes in oceanic carbon storage in glacial-to-interglacial transitions as well as in the present context of increasing atmospheric CO2 concentrations.

1 Introduction

During the last four glacial maxima, atmospheric CO2 (henceforth pCO2atm) 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 CO2 that was removed from the atmosphere was stored in the glacial ocean. In addition, studies of paleo-proxy records indicate that carbon storage in the glacial terrestrial biosphere was smaller compared to in interglacial climate (Shackleton1977; Duplessy et al.1988; Curry et al.1988; Crowley1995; Adams and Faure1998; Ciais et al.2012; Peterson et al.2014). During deglaciation, radiocarbon evidence indicates that CO2 was rapidly released 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 contributors to increased glacial oceanic storage. As glacial climate was substantially colder than interglacial climate, the global averages of surface and ocean temperature at the Last Glacial Maximum (LGM) are estimated to have been 3–8 and 2.0–3.2 C colder respectively, than the pre-industrial averages (Stocker2014; Headly and Severinghaus2007; Bereiter et al.2018). Due to the temperature effect on solubility, a colder ocean can hold more carbon. However, glacial changes in salinity partly offset the temperature effect on solubility (reviewed by Kohfeld and Ridgwell2009). 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 hypothesized that the addition of dust contributed to increased iron availability in the surface ocean and that this contributed to a strengthening of the biological sequestration of carbon in the glacial ocean compared to interglacials (Martin1990). The addition of iron would allow for more complete usage of other nutrients in 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 CO2 storage (Broecker1982a; Sarmiento and Toggweiler1984; Archer et al.2000; Sigman and Boyle2000). Other studies have pointed to changes in carbonate preservation in coral reefs and deep-sea marine sediments (Berger1982; Broecker1982a; Archer and Maier-Reimer1994). 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 Keigwin1987; Duplessy et al.1988; Stephens and Keeling2000; Marchitto and Broecker2006; Adkins2013; 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 δ13C . In addition, it has been shown that strengthening of the winds over the Southern Ocean was a likely contributor to deglacial outgassing of CO2 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), and Sigman et al. (2010). Despite the efforts of identifying the responsible processes, models have been struggling to achieve the full lowering of pCO2atm expected for a glacial.

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 may affect the glacial drawdown of pCO2atm. 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 each other. This is modelled using the average ratio of C∕P of the ocean organic matter originally observed by Redfield (1963) (C∕P=106∕1) or adjustments to these suggested in follow-up studies (Takahashi et al.1985; Anderson and Sarmiento1994).

We investigate whether allowing for a flexible C∕P stoichiometric ratio increases model ocean CO2 storage in a glacial-like climate. This possibility was suggested in studies by Broecker (1982b), Archer et al. (2000), and Galbraith and Martiny (2015), but the implications of Redfield versus flexible C∕P for glacial ocean carbon storage have not previously been tested in an Earth system model. This type of non-Redfieldian dynamics was applied in the model used in Eggleston and Galbraith (2018) and Galbraith and de Lavergne (2018), but their results were not analysed in terms of difference from a Redfield model version. In addition, Buchanan et al. (2018) explored the importance of dynamic response of ocean biology, such as flexible stoichiometry, for modelled ocean biogeochemistry in pre-industrial simulations. They found that the dynamic response was fundamental for stabilizing the response of ocean dissolved inorganic carbon (DIC) to changes in the physical circulation state.

We conduct a sensitivity study, where we apply glacial-like changes to an interglacial control state. Changes in radiative forcing, albedo, wind-forced circulation, remineralization depth, and dust are applied separately and in combination. Here, changes in radiative forcing and albedo serve to cool the climate and to mimic glacial temperature and ice conditions. The surface wind stress is reduced in the polar regions, in order to pursue reduced AABW ventilation as suggested by paleo-proxy evidence (Menviel et al.2017; Skinner et al.2017). Ocean cooling reduces the degradation rate of sinking particulate organic carbon, which increases the average depth of remineralization of organic carbon (Matsumoto2007). Drier glacial climate resulted in increased dust deposition, and thereby iron flux, to the ocean (Martin1990). All these changes act to increase ocean carbon storage and thereby reduce pCO2atm. The applied changes are not expected to induce a full glacial maximum model state, but they allow us to explore several important effects on the biological and solubility carbon pumps and produce a state with glacial-like climate conditions.

We apply the perturbations in two different versions of the Earth system model cGENIE: an original version using fixed Redfield stoichiometry of C∕P (Ridgwell et al.2007) plus a co-limitation by iron (Tagliabue et al.2016) for biological production and a modified version using the non-Redfieldian, nutrient-concentration-dependent C∕P stoichiometry suggested by Galbraith and Martiny (2015) and iron co-limitation.

We show that flexible C∕P stoichiometry allows a larger glacial ocean CO2 storage, as predicted by the box-model study of Galbraith and Martiny (2015), and that flexible stoichiometry has the largest impact for perturbations in remineralization depth and dust forcing. Additionally, we show that flexible stoichiometry allows for increased ocean carbon storage without decreasing the storage of preformed nutrients in the deep ocean.

2 Methods

2.1 Model description

cGENIE is an Earth system model of intermediate complexity, with a 3D frictional-geostrophic ocean (36×36 equal-area horizontal grid, 16 depth levels), 2D energy–moisture balance atmosphere with prescribed wind fields, and interactive atmospheric chemistry and ocean biogeochemistry. The model code and user handbook can be found in the cGENIE GitHub repository (cGENIE GitHub repository2019). 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 Sect. 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.52019, see Code availability for details).

2.2 Stoichiometry

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/O2=1/106/-138, where O2 is dissolved oxygen (nitrogen is assumed only implicitly for the purpose of accounting for organic matter creation and remineralization-related alkalinity transformations in the ocean; Ridgwell et al.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.2005; Galbraith and Martiny2015; Yvon-Durocher et al.2015; Tanioka and Matsumoto2017; Garcia et al.2018a; Moreno et al.2018). We test the importance of this variability for glacial ocean CO2 storage by running the same experiments with the fixed Redfield stoichiometry version of cGENIE and with a model version where we have implemented the linear regression model presented by Galbraith and Martiny (2015) (Eq. 1). These two model versions are henceforth denoted RED and GAM.

The flexible stoichiometry in GAM depends on the ambient concentration of dissolved phosphate ([PO4]) in the water:

(1) P / C = 1 / PO 4 144.9 µ molL - 1 + 0.0060 - 1 .

This relation shows that, when [PO4] is low, organisms bind more C per atom of P than they do under high-[PO4] conditions (Fig. 1). Equation (1) is applied in cGENIE in the calculations of biological C uptake at the surface ocean based on the surface concentration of PO4. In Sect. 3.1.2, Eq. (1) is also used to translate model surface PO4 fields to the corresponding surface C∕P ratios for the organic matter produced in each grid cell. Note that, while we change the ratio of C∕P, the ratio of C∕O2 remains the same in all experiments. As a result, the P∕O2 ratio changes between experiments.

Figure 1Flexible stoichiometry C∕P (y axis) dependent on the P concentration (µmol L−1) (x axis), as described by Eq. (1). Here, we extend the relationship beyond the observational interval 0–1.7 µmol L−1 (bounded by dashed line) which forms the basis of the relation derived by Galbraith and Martiny (2015).


2.3 Experiments

We start all experiments from an interglacial–modern (IG-M) control state, which has been run for 10 000 years to steady state, using either Redfield (CtrlRED) or variable (CtrlGAM) stoichiometry. The control states have a prescribed pCO2atm of 278 ppm and the same climate (Table 1). However due to the differences in C∕P, they have different ocean carbon inventories (Table S1 in the Supplement). In CtrlGAM, 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 CtrlRED. This also suggests 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 Ctrl121).

Table 1List of experiments. Ensemble member acronyms, short descriptions of what the simulation tests, and specifications of parameter settings for each ensemble member. The pCO2atm is either prescribed to pre-industrial level (PI) = 278 ppm or freely varying with changes in climate and ocean circulation. The radiative forcing is either coupled to the pCO2atm of the atmospheric chemistry module of the model or fixed at a value corresponding to pCO2atm=185 ppm. The zonal albedo profile is either representative of IG-M conditions or of the LGM. The wind stress is either IG-M or has an adjusted peak in wind stress at ±50 N. The IG-M remineralization length scale (RLS) is 590 m. If RLS is changed, it is multiplied by a factor (fr). Dust forcing is either IG-M or representative of LGM. Each experiment is conducted using model versions with C∕P fixed at 106∕1 (Redfield1963,  denoted RED), C∕P variable with surface ocean PO4 concentration (Galbraith and Martiny2015, denoted GM15), and C∕P fixed at 121∕1 (denoted 121; see Sect. 2.3).

a Calculated from the HADCM3 LGM (21 ka) simulation of Davies-Barnard et al. (2017). bSee Lauderdale et al. (2013) for example of reduced peak wind profile for the Southern Hemisphere. c Re-gridded LGM dust flux from Mahowald et al. (2006). d fr denotes multiplication factor for remineralization length scale. We test multiplication factors between 0.75 and 1.75, corresponding to a change in RLS between 25 % and +75 %.

Download Print Version | Download XLSX

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 Sect. 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) of

    • radiative forcing corresponding to LGM CO2= 185 ppm and

    • zonal albedo profile representative of LGM (calculated from the LGM climate simulation of Davies-Barnard et al.2017);

  • physical perturbations (weaker overturning) of

  • biological perturbations of

By applying the above perturbations, we aim to approach, but not fully resolve, some of the characteristics of the LGM ocean, which appears to have had a global average ocean temperature (Toce) 2.57 ±0.24C colder than the 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 Toggweiler1984; Martin1990; Sigman and Boyle2000). We also aim to increase carbon retention in the deep ocean (Muglia et al.2018).

The physical perturbations serve to achieve a colder climate (Toce cools by 2.1 C compared to Ctrl, thus 80 % of the observed 2.6 C; see Sect. 3.2.1) and weaker overturning (see Sect. 3.2.2 and Fig. 2) with a longer residence time of the Antarctic 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 also 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).

Figure 2The Eulerian component of the global (a, b), Atlantic (c, d), and Pacific (e, f) ocean meridional overturning stream function (1 Sv =106 m3 s−1) of CtrlGAM (a, c, e) and GLcombGAM (b, d, f). Note that the eddy-induced transport of tracers is taken into account through a skew-diffusive flux (Griffies1998) that is present in the velocity fields used to compute the Eulerian stream function.


The biological perturbations serve to achieve a more efficient biological pump, which is connected with increased retention of nutrients and carbon in the deep ocean, and lower surface nutrient concentrations in productive regions (see Sect. 3.2.3 and 3.2.4). With flexible stoichiometry, lower surface nutrient concentrations result in higher C∕P ratios, further increasing the export production and thereby 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 Sect.

2.3.1 Physical perturbations

We change the physical conditions for climate by changing radiative forcing and albedo to LGM-like conditions and denote these changes as LGMphy. We set the radiative forcing in the model to correspond to an atmosphere with 185 ppm CO2 instead of 278. However, we allow the pCO2atm 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 PMIP LGM simulations, allow 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 cause the global ocean average temperature (Toce) to decrease by 2.1 C compared to Ctrl (see Sect. 3.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 WNA ×0.5). We use the Southern Ocean wind profile of Lauderdale et al. (2013), where the peak westerly wind strength at 50°S has been halved compared to the control state (see Sect. 2.3.3). The winds north and south of the peak are reduced accordingly to give a continuous profile (see Fig. 2a 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 Sect. 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 et al.2013, 2016; Kohfeld et al.2013). To avoid an expansion of the North Atlantic Deep Water (NADW) overturning cell that would be inconsistent with the glacial ocean (Duplessy et al.1988; Lynch-Stieglitz et al.1999; Curry and Oppo2005; Marchitto and Broecker2006; 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 exchange with the atmosphere.

Table 2Atmospheric CO2 (pCO2atm, ppm); global ocean averages of temperature (Toce, C), P*, and O2 (µmol kg−1); Atlantic overturning stream function (ψ, sverdrups (Sv)); and maximum and minimum north of 30 N, for observations, and for selected ensemble members in each model version (RED and GAM). Climatological modern-day Toce and O2 were computed using the World Ocean Atlas 2018 (Locarnini et al.2018; Garcia et al.2018c). P* for the modern-day ocean is estimated by Ito and Follows (2005). Average modern-day AMOC strength is estimated by McCarthy et al. (2015) from the RAPID-MOCHA array at 26 N (corresponding to Atlantic ψmax in the model). Note that observed pCO2atm is given for pre-industrial (PI) climate, as we do not model anthropogenic release of CO2.

Download Print Version | Download XLSX

2.3.2 Biological perturbations

In the ocean, phytoplankton growth rates and remineralization of particulate organic carbon are processes that both work more slowly at colder temperatures (Eppley1972; 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 production (i.e. the amount of C captured by primary production that leaves the surface ocean without being remineralized) (Matsumoto2007). However, Matsumoto (2007) shows that the effect of slower remineralization dominates the effect on export production. It has therefore been hypothesized that the cooling of the glacial ocean led to a deepening of the remineralization length scale (henceforth denoted RLS) in the ocean and thereby more efficient retention of organic carbon in the deep ocean (Matsumoto2007; Chikamoto et al.2012), which in turn caused a lowering of pCO2atm. 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 remineralization. Deeper remineralization also results in increased nutrient retention 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 fr (RLS × fr; see Sect. 2.3.3 ).

Increased dust forcing leads to increased iron (Fe) 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 (Martin1990; Moore et al.2013). The variable stoichiometry 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 (Martin1990; Moore et al.2013). We apply the regridded LGM dust fields of Mahowald et al. (2006) and denote this change as LGMdust.

2.3.3 Sensitivity experiments and combined simulations

In the sensitivity study, for each of the three C∕P parameterizations, we first change one forcing at a time (see Table 1). We run simulations where we individually apply the LGM boundary conditions for radiative forcing (LGMrf), albedo (LGMalb), and dust (LGMdust) and one simulation with halved wind stress near the poles (50 and 50 N) (WNS ×0.5). For the remineralization length scale (RLS), we test a range of values of the multiplication factor RLS ×fr, where fr={0.75,1.25,1.75}. This allows us to test the sensitivity to deep-ocean retention of organic carbon.

We then run simulations where we combine several changes in forcing. We get a colder climate simulation (LGMphy) by combining LGMrf and LGMalb. In simulation Acomb, we combine LGMrf, LGMalb, LGMdust, and RLS ×1.25 (see Table 1). Kwon et al. (2009) show that small changes in remineralization depth can cause substantial changes in pCO2atm. With the RLS deepening of 25 %, we keep the corresponding changes in pCO2atm from exceeding the ∼20–30 ppm obtained in other studies (Matsumoto2007; Menviel et al.2012). We finally run a glacial-like simulation GLcomb (see Sect. 3.3), which is similar to Acomb but also includes the change in wind stress WSN ×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.

2.4 Observations

For comparison and validation of model results, we use records of ocean-state variables from climatological mean fields of modern observations and proxy observations from the LGM.

Modern data of ocean temperature, oxygen, and nutrients are retrieved from the World Ocean Atlas 2018 (Locarnini et al.2018; Garcia et al.2018c, b) and we use the proxy estimates of LGM ocean temperature from Bereiter et al. (2018). Average modern-day strength of the Atlantic meridional overturning circulation (AMOC) is estimated by McCarthy et al. (2015) from the RAPID-MOCHA array at 26 N.

We use model–data comparison of benthic δ13C to assess the statistical similarity (correlation) between both the model control state and glacial-like state (see Sect. 2.3) to benthic δ13C data representing the late Holocene (0–6 ka, HOL) and Last Glacial Maximum (19–23 ka, LGM) (Peterson et al.2014). Locations of core sites can be seen in Fig. S1 in the Supplement (see also Fig. 1 in Peterson et al.2014). Note that LGM benthic δ13C is more 13C-depleted than that in the Holocene due to the addition of 13C-depleted terrestrial carbon to the glacial ocean (Shackleton1977; 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 δ13C of 0.32±0.20 ‰ 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.

2.5 Nutrient utilization efficiency

The extent to which biology succeeds in using the available nutrients can be determined by calculating the nutrient utilization efficiency P* (Ito and Follows2005; Ödalen et al.2018),

(2) P * = P rem P tot ,

which is the fraction of remineralized (Preg) to total (Ptot) nutrients (in this case, PO4) in the ocean. Overlines denote global averages. Remineralized nutrients have been transported from the surface to the interior ocean by the biological pump, and Prem is given by

(3) P rem = P tot - P pre .

Here, Ppre is preformed PO4 – the concentration of PO4 that was present in the water parcel as it sank, and 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

(4) DIC rem = DIC tot - DIC pre .

DICrem is used to compute the ocean storage of remineralized acidic carbon (ACrem; see Appendix A). ACrem is biological carbon that entered the ocean in the form of CO2 in soft tissue (as opposed to carbonates in hard tissue), measured independently of oxygen consumption and/or remineralized phosphate.

In cGENIE, Ppre and DICpre are modelled as passive tracers (Ödalen et al.2018). Hence, we can use the model output for Ppre in Eq. (3) to compute P* (Eq. 2).

In a model with a fixed Redfield ratio, P* determines the effect of the biological pump on pCO2atm. For example, it has been found that a higher P* in the initial state gives a lower potential for drawdown of pCO2atm 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 carbon retained in the deep ocean is not necessarily proportional to P*.

Figure 3Surface PO4 concentration (µM) (a, c) and corresponding C∕P as calculated using parametrization of Galbraith and Martiny (2015) (b, d). Panels show CtrlGAM (a, b) and climatological mean fields (World Ocean Atlas 2018, Garcia et al.2018b) (c, d).

3 Results

3.1 Control states

3.1.1 Ocean temperature and circulation

As the three control states CtrlRED, CtrlGAM, and Ctrl121 are driven by the same physical forcings and have the same pCO2atm, they have the same ocean circulation pattern (Fig. 2a, c, e and Tables 2 and S2) and climate (exemplified by global ocean average temperature (Toce) in Tables 2 and S2). The surface ocean nutrient fields are fairly similar, with small differences due to the different C∕P parameterizations (compare Figs. 3a and S2). The strength of the Atlantic meridional overturning circulation (AMOC), diagnosed in the model as the maximum of the Atlantic meridional overturning stream function deeper than 1000 m, is 14 Sv (1 Sv =1×106 m3 s−1) in all control states (Tables 2, S2). Results from the RAPID-MOCHA array at 26 N suggest an average AMOC strength of 17.2±0.9 Sv (McCarthy et al.2015); thus our control-state AMOC is a little bit weaker than in the present-day climate. The observational estimate for Toce according to the World Ocean Atlas 2018 (Locarnini et al.2018) is 3.49 C, thus comparable to the 3.56 C of our Ctrl simulations. The surface nutrient concentrations of our control-state CtrlGAM (Fig. 3a) compare reasonably well with observed surface ocean concentrations of PO4 (Fig. 3c), with some underestimation in the Pacific equatorial region, the North Pacific Ocean, and the Labrador Sea. The agreement with observations is better for CtrlGAM than for CtrlRED (Fig. S2).

3.1.2 Surface nutrient distribution and C∕P ratios

In Fig. 3 we see that surface PO4 fields (left-hand column) and the corresponding fields of surface C∕P ratios (as given by Eq. (1), right-hand column) of CtrlGAM (Fig. 3a, b) and of climatological mean PO4 fields (Fig. 3c, d) are similar in their pattern as well as in the magnitudes of values. Note that high concentrations of PO4 correspond to low C∕P ratios, and vice versa. The highest climatological mean PO4 concentrations in the northern and equatorial Pacific are not fully reproduced by the model, but the pattern is reproduced well. In the surface C∕P field of CtrlGAM (Fig. 3b), we see the signature of very high nutrient concentrations in the Southern Ocean (>1µmolL-1, Fig. 3b) as a band of low ratios, with the most extreme values near the Antarctic continent, as seen in the climatological mean fields (Fig. 3d).

The nutrient utilization efficiency P* (Eq. 2) in the three control states differs by a few percent: 0.43, 0.46, and 0.42 in CtrlRED, CtrlGAM, and Ctrl121 respectively (Table 2). The fraction of DICrem in DICtot (see Sect. 2.5, Eq. 4, Table S1) is 0.065, 0.077, and 0.072 in CtrlRED, CtrlGAM, and Ctrl121 respectively.

3.1.3 Ocean dissolved O2

The most apparent difference between CtrlRED and CtrlGAM is in deep-ocean oxygen concentrations, where the global ocean average dissolved O2 concentration (O2) in CtrlGAM (144 µmol kg−1) is lower than in CtrlRED and Ctrl121 (166 and 152 µmol kg−1 respectively). Compared to climatological mean fields (World Ocean Atlas 2018, Fig. 4a–b), both CtrlRED (Fig. 4c–d) and CtrlGAM (Fig. 4e–f) agree reasonably well with the real ocean. CtrlGAM appears to better capture the equatorial oxygen minimum in the Atlantic basin than CtrlRED but is too low in the North Pacific. In CtrlGAM (Fig. 4f), the North Pacific is markedly lower in oxygen than in CtrlRED (Fig. 4d) and even anoxic in the oxygen minimum zone (OMZ). This should be kept in mind when analysing the oxygen sections of the glacial-like states GLcombRED (Fig. 4g–h) and GLcombGAM (Fig. 4i–j). Global averages for dissolved O2 are given in Table 2.

Figure 4Sections of O2 concentration (µmol kg−1) along 25 W in the Atlantic basin (left-hand column) and along 135 W in the Pacific basin (right-hand column). Panels show climatological mean fields (World Ocean Atlas 2018, Garcia et al.2018c) (a, b) and model states CtrlRED (c, d), CtrlGAM (e, f), GLcombRED (g, h), and GLcombGAM (i, j).


3.1.4 Ocean δ13C

By comparing CtrlRED δ13C and Holocene (0–6 ka, HOL) benthic δ13C values, we estimate a global model–data correlation of 0.78 (Table S3). The modern-day Atlantic Ocean has a distinctive spatial δ13C pattern (Figs. 5a, S1) with 13C-enriched values in the intermediate-depth (<2 km) North Atlantic and Nordic Seas and 13C-depleted values in the deep (>2.5 km) South Atlantic. While the model produces a weaker gradient than the observed HOL Atlantic Ocean (corr. 0.50, Table S3) the model correlates well with eastern Atlantic δ13C records (Fig. S1). For the Indo-Pacific, the weaker benthic δ13C gradient is well represented by the model (Fig. 5e). This pattern emerges mainly due to 13C-depleted, biologically sourced carbon that is accumulated in the weak circulation region of the interior North Pacific (Matsumoto et al.2002). However, Indo-Pacific δ13C values of CtrlRED are overall lower than the HOL observations. The overall model–data correlation for the Indo-Pacific is 0.39 (Table S3). Comparing the control states of the RED and GAM model versions, δ13C patterns (Figs. 5a, e and S3a, e) and model–data correlations with HOL observations (Table S3) are similar between the model versions, with somewhat lower correlations for GAM.

Figure 5Model ocean δ13C (contours) compared to the two proxy record time slices (HOL and LGM) of benthic δ13C (circles) of Peterson et al. (2014). The upper half of the figure shows the Atlantic Ocean (a–d), while the lower half shows the Pacific Ocean (e–h). 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 (a, c, e, g) and the right-hand column shows GLcombRED (b, d, f, h). The rows show, from top to bottom, (a, b) HOL Atlantic, (c, d) LGM Atlantic, (e, f) HOL Pacific, and (g, h) LGM Pacific. Note that, before we compare GLcombRED to LGM observations (d, h), a constant of 0.32 ‰ is subtracted from the simulated δ13C, to account for terrestrial release of δ13C-depleted terrestrial carbon which is not modelled. The corresponding comparison for model version GAM is shown in Fig. S3.


3.2 Sensitivity experiments

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 pCO2atm. The resulting steady-state global average values for temperature (Toce), dissolved oxygen (O2), and nutrient utilization efficiency P*, as well as the maximum and minimum of the Atlantic meridional overturning stream function, are listed in Table 2.

3.2.1 Radiative forcing and albedo

In the simulations where radiative forcing and albedo are changed to represent LGM conditions (LGMrf + LGMalb = LGMphy), the reductions in pCO2atm are similar in the RED and the GAM model versions. In LGMphy, the resulting pCO2atm is 245.4 and 244.9 ppm, thus showing 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. 7a), and the associated changes in C∕P (Fig. 8a), 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 pCO2atm can mainly be attributed to the increase in solubility carbon (Csat,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 (Cdis) (Ödalen et al.2018). Ocean cooling amounts to 2.1 C in LGMphy compared to Ctrl (Toce in Table 2). In cGENIE, the increase in Csat,T associated with ocean cooling corresponds to 7±1.5 ppm C−1 (Ödalen et al.2018, Supplement Fig. S1). Thus, in LGMphy Csat,T and Cdis should contribute roughly 40 % and 60 % respectively of the change in pCO2atm. Note that cGENIE underestimates the true effect of ocean cooling on solubility, due to a temperature restriction on the solubility constants (Ödalen et al.2018). This temperature restriction limits solubility from changing below 2.0 C. In this case, the solubility effect on pCO2atm of reducing Toce by 2.1 C (from 3.6 to 1.5 C, ∼15 ppm) is thus comparable to only 1.6 C of cooling (from 3.6 to 2.0 C, ∼11 ppm), and the solubility effect is underestimated by ∼4 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 Ridgwell2009). This effect is consistent for any choice of C∕P parametrization and is therefore not explored further.

Figure 6Resulting CO2 anomaly, with respect to the control-state 278 ppm, of the sensitivity experiments LGMphy (plus symbol), RLS ×1.25 (times symbol), LGMdust (upward arrowheads), and WSN ×0.5 (circles) and of the combined experiments Acomb (stars) and GLcomb (downward arrowheads). Results of the different model versions RED, 121, and GAM are shown in red, yellow, and blue respectively. The vertical dashed line separates simulations without (left) and with (right) wind perturbation.


Figure 7Surface PO4 anomaly (µM), with respect to surface concentration of PO4 in CtrlGAM (Fig. 3a), for (a) LGMphyGAM, (b) WSN ×0.5GAM (c) RLS ×1.25GAM, (d) LGMdustGAM, (e) AcombGAM, and (f) GLcombGAM. Changes in surface nutrient fields are similar for all three model versions (RED, GAM, 121); thus only GAM is shown.

3.2.2 Reduced wind forcing

When the peak of Southern Ocean (henceforth SO) winds is reduced, the strength of the overturning circulation of AABW decreases (see difference in Southern Hemisphere overturning stream function between Fig. 2a 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 near the surface for longer and loses more nutrients before being subducted. This decreases the SO concentration of preformed phosphate in WNS ×0.5RED compared to Ctrl, as seen in Fig. 7b, and increases the nutrient utilization efficiency P* (Table 2, Fig. 9). This leads to a drawdown of pCO2atm of 12.9 ppm compared to CtrlRED (Fig. 6). As the nutrient concentration in the SO decreases (Fig. 7b), the flexible C∕P ratio (Fig. 8b) leads to an increased carbon capture efficiency in GAM compared to RED (see GLcomb-Ctrl of biologically sourced carbon (ACrem) in Fig. 9), which is partly compensated for by a reverse effect in the Pacific equatorial region. Consequently, in WSN×0.5GAM, we get a reduction of pCO2atm of 16.3 ppm compared to CtrlGAM (Fig. 6). Hence, for halved peak wind stress at ±50 N, the flexible stoichiometry increases the drawdown by ∼26 %.

Figure 8Surface C∕P anomaly, with respect to C∕P of CtrlGAM (Fig. 3b), for (a) LGMphyGAM, (b) WSN ×0.5GAM (c) RLS ×1.25GAM, (d) LGMdustGAM, (e) AcombGAM, and (f) GLcombGAM.

Figure 9Remineralized acidic carbon (ACrem=DICrem-(1/2)×ALKrem; see Appendix A) as a function of P*. Simulations using model versions RED and GAM are shown in red and blue respectively. Different symbols indicate the sensitivity experiments, listed in the panel on the right-hand side. Red and blue lines show linear least-squares fits to the separate ensembles for RED and GAM. Panels illustrate (a) the deviation from the least-squares fit of the GAM ensemble as opposed to the RED ensemble, (b) the linear fits extrapolated to the origin, and (c) a zoom-in around the origin, where the RED (red) line goes through the origin, while the GAM (blue) does not.


3.2.3 Remineralization length scale

When the remineralization length scale (RLS) increases, the biological material reaches deeper before it is remineralized, and it takes longer for it to be returned to the surface. Therefore, more of the biologically sourced carbon (ACrem) and nutrients are present in the deep ocean at any given time, leading to an increase in P* (Table 2, Fig. 9) and a decrease in pCO2atm (Fig. 6). The deeper we make the RLS, the bigger the drawdown of pCO2atm – in RLS ×1.25RED and RLS ×1.75RED pCO2atm decreases by 14 and 33 ppm respectively compared to CtrlRED. In GAM, the drawdown in each experiment is increased by an additional ∼30 %; thus RLS ×1.25GAM and RLS ×1.75GAM see a reduction of pCO2atm of 18 and 44 ppm respectively compared to CtrlGAM (Table S2). Our changes in RLS cause very small, but global, changes in PO4 concentrations (global average anomaly =0.016 µM, RLS ×1.25 in Fig. 7c), which, through the small resulting changes in C∕P (Fig. 8c), still contribute to the additional drawdown of pCO2atm in GAM.

In a sensitivity test where we make the RLS 25 % shallower (which would be representative of a warmer climate compared with Ctrl), the pCO2atm increases by 18 and 23 ppm in RED and GAM respectively compared to their control states (see RLS ×0.75RED and RLS ×0.75GAM in Table S2). Interestingly, the response in pCO2atm is again ∼30 % larger in GAM. The variable stoichiometry thus amplifies the effect on pCO2atm by any change in RLS. The potential implications of this result for warm climate scenarios are further discussed in Sect. 4.1.

3.2.4 Dust forcing

The simulations with LGM dust forcing (LGMdust, Table 1) show the largest difference in pCO2atm between the RED and the GAM. In LGMdustRED, pCO2atm decreases by 16 ppm compared to CtrlRED, whereas LGMdustGAM sees a reduction of 21 ppm compared to CtrlGAM (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 organic material that is exported out of the surface ocean (see Sect. 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 Prem (thus, P* increases, Table 2, Fig. 9). This reduces the surface nutrient concentrations in these areas (Fig. 7d). In the GAM model version, this is followed by increased C∕P ratios in these areas (Fig. 8d), resulting in a lower pCO2atm in LGMdustGAM than in LGMdustRED. The largest anomalies in PO4 concentrations, and consequently in C∕P, are observed in the 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 (10Be, 230Th, 231Pa) from the LGM (e.g. Kumar et al.1995; Kohfeld et al.2005).

3.3 Combined experiments

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

3.3.1 Ocean temperature and circulation

In the glacial-like simulations, GLcombRED and GLcombGAM, the global average ocean temperature (Toce) is 1.7 C lower than in the respective control states (Table 2). Headly and Severinghaus (2007) estimate LGM Toce to have been 2.7±0.6C colder than the modern ocean, while a more recent estimate by Bereiter et al. (2018) constrains Toce to 2.57±0.24. GLcomb is thus just outside the 1 standard deviation limit of the warm end of the Headly and Severinghaus (2007) estimate for the LGM. In Acomb, Toce is 2.1 C cooler than Ctrl; thus this simulation falls within the uncertainty of the Headly and Severinghaus (2007) estimate. Compared to the Bereiter et al. (2018) estimate, our combined experiments GLcomb and Acomb achieve 64 % and 82 % of the glacial–interglacial difference in Toce respectively. As anticipated, our combined forcings do not induce a full glacial maximum state but rather 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 reduction. For example, the AMOC (here measured as the maximum of the Atlantic overturning stream function) reduces in strength by ∼15 %, from 14 to 12 Sv (Table 2). The global meridional overturning stream function reveals that the SO overturning cell sees a reduction in transport (Fig. 2d), which is associated with weaker upwelling and thus longer residence time for AABW, as hypothesized 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).

3.3.2 Surface nutrient distribution and C∕P ratios

In the surface nutrient anomalies (GLcomb-Ctrl, shown for GAM in Fig. 7f), 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 CO2 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 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. 7f), which coincides with an increase in sea ice in this area (Fig. S4). 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. 7a). North of the biogeochemical divide, increased aeolian dust flux increases the productivity of the biology, which reduces P in the surface compared to the control (Fig. 7d). In combination with circulation changes, resulting from the reduced SO wind stress (Fig. 7b), and deeper remineralization (Fig. 7c), P concentrations in the subantarctic region are strongly reduced (Fig. 7f). In the equatorial and North Pacific, there is also a reduction of P (Fig. 7f), mainly due to the increased dust flux (Fig. 7d). 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 bands (Fig. 8f). The organic matter that is exported out of the upper layer (henceforth referred to as export production) in GLcombGAM has a global average C∕P ratio of 134∕1.

3.3.3 Ocean dissolved O2

Despite colder conditions, which allow for more dissolution of O2, the reduction in O2 is evident in the GLcomb simulations (Fig. 4). This mirrors the increase in ACrem (Fig. 9). The O2 reduction is about 50 % larger in GAM compared to RED. As the initial-state CtrlGAM is already lower in oxygen than CtrlRED (144 µmol kg−1 compared to 166 µmol kg−1), and variable stoichiometry allows for additional ocean storage of organic carbon, the end-state O2 is drastically lower in GLcombGAM (74 µmol kg−1) compared to GLcombRED (122 µmol kg−1).

3.3.4 Ocean δ13C

In GLcombRED (contours in Fig. 5b, d), the Atlantic north–south gradient in δ13C is stronger than in CtrlRED (contours in Fig. 5a, c). This strong gradient is not observed in the Holocene Atlantic δ13C data (dots in Fig. 5a, b), but it is prominent in the LGM time slice (dots in Fig. 5c, d). The LGM observations are reproduced well in GLcombRED (corr. 0.62, Table S3), especially in the east Atlantic (corr. 0.77). When we correct GLcombRED for the absence of injected terrestrial carbon, we see clear similarities with LGM observations (Fig. 5d), though the southernmost cores still indicate more 13C-depleted conditions than the model. In the Indo-Pacific, GLcombRED is too 13C-depleted compared to LGM observations, particularly with the correction for the absence of a terrestrial signal (Fig. 5h), and the model–data correlation is poor (0.05). For the Indo-Pacific, the model–data correlation for Holocene data is similar between GLcombRED and CtrlRED (0.24 and 0.39 respectively, Table S3). This suggests that the poor correlation with LGM benthic records is simply due to our changes in forcings being insufficient to achieve the required rearrangements in Indo-Pacific circulation patterns. In GLcomb, similarly as in Ctrl, there is very little difference between the RED and GAM model versions in terms of δ13C patterns (compare Fig. 5d, h to Fig. S3d, h) and model–data correlations (Table S3). However, there are overall lower values of δ13C in GAM, reflecting the larger storage of biologically sourced carbon in the deep ocean in this model version.

3.3.5 Atmospheric CO2

In the combined experiments Acomb and GLcomb, pCO2atm decreases strongly compared to the control state (reduction of between −56 and −80 ppm, Fig. 6). This is partly a result of colder conditions (Table 2), which lead to increased solubility for CO2 in seawater and an expanded sea-ice cover (Fig. S4) which restricts air–sea gas exchange. The latter slows down the equilibration of the CO2-rich upwelling water with the atmosphere before it is subducted into the deep ocean and thus causes increased Cdis. In Acomb and GLcomb, changes in biological production (see Sect. 3.3.2) and storage of biologically sourced carbon (Fig. 9) also contribute strongly to the reduced pCO2atm.

The combined experiments show a striking difference in pCO2atm between the model versions RED and GAM; in GLcombRED and GLcombGAM, we achieve drawdown of pCO2atm 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 and 173 PgC respectively (Table S1). The drawdown is thus 25 % larger in GAM than in RED. In Acomb, where the perturbation in wind stress is omitted, drawdown of pCO2atm is smaller than in GLcomb and the difference between model versions is less pronounced, but there is still a 14 % difference between AcombRED and AcombGAM (−56 and −63 ppm respectively (Fig. 6).

4 Discussion

4.1 Accounting for variable C∕P in ocean carbon cycle models

The representation of ocean biology in general circulation models (GCMs) tends to be oversimplified (Le Quéré et al.2005) 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 parametrization, or possibly a power law such as that described by Tanioka and Matsumoto (2017), an additional facet of the complexity of ocean biology can be implicitly accounted for at a relatively small computational cost. Previous model ensemble studies have shown that this type of dynamical response of the biology to changes in the modelled ocean state can improve the model's ability to realistically simulate ocean biogeochemistry (Buchanan et al.2018). In pre-industrial and future simulations respectively, Buchanan et al. (2018) and Tanioka and Matsumoto (2017) find that the flexible stoichiometry acts to stabilize the response of ocean DIC to changes in the physical (circulation) state. In our glacial-like simulations, we find that the response of ocean DIC, and thus pCO2atm, to the combined perturbations is greater in the simulations with flexible stoichiometry. Nonetheless, our study confirms the potential importance of dynamical biological response for the outcome of model studies.

The approach taken by both Galbraith and Martiny (2015) and Tanioka and Matsumoto (2017) is to adapt a function to species-independent C : P observations. This means that they account for the adaptation of plankton to the surrounding conditions, in terms of both 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 is what we use here. Ganopolski and Brovkin (2017) appear to succeed with full glacial–interglacial CO2 cycles 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 remineralization of nutrients (e.g. Kolowith et al.2001; Letscher and Moore2015), 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 [PO4]. As such, in a high-surface-[PO4] 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-Durocher et al.2015; Tanioka and Matsumoto2017; Moreno et al.2018). Up to [PO4]=1.7 µM, the GAM parametrization fits the binned observational data well (see Figs. 1 and S2 in Galbraith and Martiny2015). To account for the lack of observational data at higher [PO4] in Galbraith and Martiny (2015), we have tested the effect of saturation of the C∕P ratio at higher [PO4] in Ctrl and GLcomb simulations where we kept C/P=55/1 at [PO4]>1.7 µM. The increase in ocean carbon storage and decrease in pCO2atm between Ctrl and GLcomb are nearly identical with GAM; thus saturation of the C∕P ratio at very high [PO4] causes no noticeable impact on our results.

Our sensitivity experiments RLS ×0.75 and RLS ×1.25 reveal that the response in pCO2atm to the perturbation is enhanced in GAM compared to RED for both increased and decreased RLS. While increased RLS would be an effect of ocean cooling, and thus of interest for glacial studies, reduced RLS would be a consequence of ocean warming (Matsumoto2007). Matsumoto (2007) describes how decreased RLS would have a positive feedback on pCO2atm in a future warming climate. Our results imply that flexible C∕P could have a further reinforcing effect on this feedback. It would therefore be of interest to apply a parametrization of flexible C∕P in models used for simulations of future climate feedbacks.

4.2 Decoupling of biologically sourced carbon and nutrient utilization efficiency

Many studies have suggested increased ocean storage of organic carbon as a potentially important contributor to the low glacial pCO2atm (e.g. Sarmiento and Toggweiler1984; Martin1990; Archer et al.2000; Sigman and Boyle2000). However, biological production depends on water temperature (Eppley1972) and decreases in cold conditions. This temperature effect is parameterized in cGENIE as a local temperature-dependent uptake rate modifier proportional to e(T/15.9C), and we see overall reduced productivity in our cold climate simulations (see e.g. LGMphy, Table S1). As the climate cools, the temperature effect leads to a decrease in biological productivity and a subsequent decrease in P* (see LGMphy 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 (increased residence time of the AABW; Menviel et al.2017) and deeper remineralization (Matsumoto2007; Chikamoto et al.2012), which we apply through our perturbations in winds (Sect. 3.2.2) and RLS (Sect. 3.2.3). From our results, it is evident that variable stoichiometry can be another contributing factor. In our simulations, global average export production decreases, in terms of both POP (particulate organic phosphorus) and POC (particulate organic carbon), by 15 % between CtrlRED and GLcombRED 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 POP decreases by 17 %, the corresponding decrease in POC is only 10 %. Thus, in GLcombGAM we achieve an increase in ocean inventory of remineralized carbon, which exceeds the one of GLcombRED (Fig. 9).

Due to the competing effects between decreased export production and increased retention that are introduced by the different forcing components, the net change in Prem, 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 GLcomb compared to Ctrl, which is linearly related to a small increase in storage of ACrem 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 ACrem as well. Instead, we see a reasonably large increase in global average ACrem. In summary, we suggest this is a result of changes in surface P fields (see Fig. 7).

The reduced pCO2atm in GLcombGAM, compared to GLcombRED, can be explained by the non-linearity introduced by the local variability in C∕P. When changes in ocean circulation, remineralization 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 (Fig. S5). 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, resulting in a non-linear relationship between storage of remineralized phosphate and biologically sourced carbon (Fig. 9). In CtrlGAM, the average elemental C∕P composition is 121∕1. In GLcombGAM, this average is 134∕1 (Table S1, Fig. S5). 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 remineralization. In CtrlGAM, the global average concentration of Prem is 1.16 µmol kg−1 (compared with 1.17 µmol kg−1 in GLcombGAM). By increasing the average C∕P composition of 1.16 µmol kg−1 organic molecules from 121 to 134 (i.e. by 13 units), this causes an increase in Csoft by ∼15µmol kg−1, which corresponds to the observed increase in ACrem. From this, we conclude that, in a system where stoichiometry is variable on a local scale, ocean storage of biologically sourced CO2 can change while the amount of remineralized nutrients remains constant.

Note also that CtrlGAM has a larger inventory of DIC as well as Csoft compared to CtrlRED. Ödalen et al. (2018) found that drawdown of pCO2atm in response to a perturbation is larger when the control state has a smaller inventory of DIC and Csoft. Yet, the effect of applying the same perturbation results in a larger drawdown of pCO2atm in GAM than in RED. This is opposite of the conclusions of Ödalen et al. (2018). The reason is that the flexible stoichiometry in effect increases the drawdown potential, which more than compensates for the increased carbon inventory in the control state. In GLcombGAM, the C∕P in global average export production increases from the control 121∕1 to 134∕1, which reflects the increased storage of organic carbon allowed by the flexible stoichiometry.

4.3 Implications of flexible C∕P for deep-ocean oxygen

Studies have shown that deep-ocean O2 concentrations were lower during the LGM than in the Holocene (e.g. Bradtmiller et al.2010; Jaccard and Galbraith2012; Galbraith and Jaccard2015). 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, CtrlGAM (Fig. 4e) reproduces the climatological mean extent of the oxygen minimum zone (OMZ) (Fig. 4a) better than CtrlRED (Fig. 4c) does. In the Pacific Ocean, the O2 gradient in the climatological mean fields (Fig. 4b) is more gradual compared to that of the control states (Fig. 4d, 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 GLcombRED (Fig. 4h) compared to CtrlRED (Fig. 4c), in agreement with the findings of Hoogakker et al. (2018). In GLcombGAM, oxygen depletion is too extensive (see below), but the tendency of vertical expansion compared to the control state is present here as well.

In our glacial-like states, we see a significant reduction of O2 compared to the control state in both RED and GAM (Fig. 4d and e). This response is expected when we apply increased dust deposition and deeper remineralization (Galbraith and Jaccard2015), and the direction of the overall response of deep-ocean O2 to the glacial-like forcings is in line with observations. The reduction in O2 is stronger in GAM, due to the larger storage of respired carbon in this model version.

Finally, it should be noted that CtrlGAM (Fig. 4c) 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 CtrlRED (Fig. 4b) and in climatological mean fields. In CtrlGAM, a large part of the interior North Pacific is anoxic (Fig. 4c), while the climatological mean fields (Fig. 4a) 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. Among other factors, model ocean oxygen conditions are also dependent on deep water formation characteristics of the model (Galbraith and de Lavergne2018). The deep water formation characteristics of a model affect the amount of time available for remineralization 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 Lavergne2018). In cGENIE, this effect is small enough not to cancel the increased O2 consumption caused by the higher average C∕P in CtrlGAM compared to CtrlRED.

In summary, the GAM version of the model quantitatively reproduces the observed O2 patterns of both the LGM and the Holocene, but with too low concentrations. If this parametrization of variable stoichiometry is to be used in cGENIE in future studies, we suggest some retuning, for example by reducing the global average concentration of nutrients, to improve the representation of observed modern-day ocean O2 concentrations in the GAM control state.

4.4 Effect of modified but fixed C∕P

Part of the observed difference in pCO2atm between GLcombRED and GLcombGAM results from a difference in global average C∕P in the control states (CtrlRED and CtrlGAM). In CtrlGAM, the average C∕P in the export production is close to 121∕1, instead of 106∕1 as in CtrlRED. We illustrate the consequences of this difference by running parallel simulations with fixed C∕P of 121∕1 (model version 121).

In a simulation with reduced wind stress (WNS ×0.5121), pCO2atm is reduced by −14.5 ppm compared to CtrlRED (Table S2). The corresponding pCO2atm drawdown in WNS ×0.5RED and WNS ×0.5GAM is −12.9 and −16.3 ppm respectively. This indicates that, for the reduced wind stress case, about half of the enhanced drawdown achieved by the flexible stoichiometry can be attributed to a difference in the mean C∕P composition of the organic material between the control states CtrlRED and CtrlGAM. Similarly, for the enhanced dust deposition experiments, LGMdust121 explains about one-third of the difference in pCO2atm between LGMdustRED and LGMdustGAM (Table S2). In the combined experiment GLcomb, we can attribute about half (−8 ppm of −17 ppm) of the observed difference in drawdown between RED and GAM to the difference in control-state average C∕P. In Acomb, the control-state difference in C∕P accounts for about two-thirds of the difference between model versions (−5 ppm of −8 ppm). Evidently, the individual perturbation simulations and GLcomb have smaller fractions of the change in pCO2atm that are due to changed average C∕P. As shown above, the simulations with 121 indicate that, depending on the change in forcing, between one-third and two-thirds 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  S2). From this, we conclude that the effects of the perturbations do not add linearly.

As outlined in Sect. 4.3, an increase in C∕P reduces deep-ocean O2 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 O2 concentration, O2, in Ctrl121 is 152 µmol kg−1, which is lower than CtrlRED (166 µmol kg−1), but higher than CtrlGAM (144 µmol kg−1). Similarly, GLcomb121 has a lower O2 than GLcombRED (96 compared to 122 µmol kg−1), but higher than GLcombGAM (74 µmol kg−1).

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 Sarmiento1994 or Takahashi et al.1985 stoichiometries, compared to Redfield1963). 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 a different model response to the same perturbation.

4.5 What can we learn from the model–data comparison of δ13C ?

Proxy records of benthic δ13C indicate a change in ocean δ13C across the deglaciation. The whole ocean deglacial change has been estimated to ∼0.35 ‰ (Peterson et al.2014; Peterson and Lisiecki2018), and the surface-to-deep gradient weakened (shown in numerous studies; see e.g. Curry et al.1988; Duplessy et al.1988; Curry and Oppo2005; Herguera et al.2010; Peterson and Lisiecki2018, and references therein). Here we compare model ocean δ13C of our simulations to the benthic δ13C records, to see how well the simulations capture the observed patterns and whether there is a difference in model–data correlation between RED and GAM.

The model–data comparison in δ13C (Sect. 3.1.4) suggests that the Ctrl simulations are overall well correlated with Holocene benthic δ13C data (HOL in Table S3). For the Atlantic, the correlation of the Ctrl simulations is higher with the LGM benthic δ13C (LGM in Table S3) than with HOL. However, the standard deviations (SDs) suggest that the Atlantic 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 δ13C gradient in the Atlantic, allowing for a closer match with LGM data in terms of SDs. A stronger gradient in δ13C and more depleted values 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 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 δ13C are, for example, brine rejection (Bouttes et al.2010, 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 we do not achieve a glacial Pacific Ocean circulation consistent with observed δ13C patterns.

Each of the two observational datasets (HOL and LGM) displays similar correlations across the two model simulations. The correlation of the HOL δ13C records with CtrlRED, CtrlGAM, GLcombRED, and GLcombGAM is in all cases between 0.76 and 0.78. Meanwhile, the correlation of LGM δ13C records with the same four simulations is in all cases between 0.55 and 0.58. That our GLcomb simulations still correlate so well with the HOL dataset 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 δ13C 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 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 CO2 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.

How δ13C 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 δ13C 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.

5 Conclusions

In this paper, we examine the potential role of variable stoichiometry in biological production for glacial ocean CO2 storage. We show that flexible C∕P composition of organic matter allows a stronger response of pCO2atm to glacial-like changes in climate, remineralization length scale, and aeolian dust flux. We conclude that variable stoichiometry may be important for glacial ocean CO2 storage and for achieving the full extent of drawdown of atmospheric CO2 in model simulations. In the experiment GLcombRED, with glacial-like climate and Redfield stoichiometry (Redfield1963), ocean carbon storage increases by 139 PgC and atmospheric CO2 decreases by 64 ppm. In GLcombGAM, with glacial-like climate changes and variable stoichiometry, the corresponding numbers are 173 PgC and 80 ppm. Hence, the drawdown of atmospheric CO2 increases by 25 % when C∕P is variable.

About half of the increased drawdown of CO2 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 ratio of remineralized to total nutrients in the deep ocean. With fixed stoichiometry, an increase in remineralized carbon is inevitably tied to a corresponding increase in remineralized nutrients.

We apply variable C∕P parameterized as a simple function of the surface water concentration of PO4, 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 [PO4] and therefore suggest a more complex power-law function, which takes into account saturation of the C∕P ratio at high concentrations of PO4. However, we found that saturation of the C∕P ratio at concentrations higher than the observational upper bound of 1.7 µM causes no noticeable impact on our results.

The representation of flexible stoichiometry used in this study (Galbraith and Martiny2015) 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 Wing2011). We show here that, for glacial–interglacial cycles, this complexity contributes to changes in atmospheric CO2 through flexible C∕P ratios. Flexible C∕P also has the potential to be an additional positive feedback of ocean warming on pCO2atm in future climate.

Appendix A: Regenerated acidic carbon

Carbon enters the ocean mainly in the form of CO2 and dissolved carbonate. Despite this, the major fraction of carbon in the ocean is 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 CO2 (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 CO2 that was supplied to the ocean via biological soft tissue, we make use of the separation of DIC and ALK into their preformed and remineralized fractions (see Sect. 2.5). Thus, we then compute the remineralized acidic carbon (ACrem) as

(A3) AC rem = DIC rem - 1 2 ALK rem .
Appendix B:δ13C in cGENIE

cGENIE represents 13C 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 δ13C in per mille relative to the standard Vienna Pee Dee Belemnite (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 Δ14C) in cGENIE, in Turner and Ridgwell (2016).

In the aqueous phase, the isotopic partitioning of carbon between CO2(aq), HCO3-, and CO32- is resolved and follows Zeebe and Wolf-Gladrow (2001) (their Sect. 3.2). The empirical fractionation factors used are from Zhang et al. (1995). The air–sea 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 (δ13CPOC), the model of Rau et al. (1996) is adapted, assuming that the isotopic signature of exported POC reflects that of phytoplankton biomass. Following Ridgwell (2001), the full equation of Rau et al. (1996) is simplified to

(B1) δ 13 C POC = δ 13 C CO 2 ( aq ) - ϵ f + ϵ f - ϵ d K Q CO 2 ( aq ) ,

where [CO2(aq)] is the ambient concentration of aqueous CO2, and δ13C(aq) is its isotopic composition. KQ is a temperature-only-dependent approximation of the full cell-dependent size and growth rate parameterization in the Rau et al. (1996) model (see Ridgwell2001). We take an intermediate value for the enzymatic isotope fractionation factor associated with intracellular C fixation (ϵf) of −25 ‰ following Rau et al. (1996, 1997), and we 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 δ13CPOC 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).

For 13C 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 δ13C of aqueous HCO3- and calcite.

Code and data availability

The code for the cGENIE.muffin model is hosted on GitHub. The specific version used in this paper, tagged as release v0.9.5, can be obtained at, last access: 12 April 2020 and is assigned a DOI (cGENIE release v0.9.52019). Configuration files for the specific experiments presented in the paper can be found in the directory cgenie.muffin/genie-userconfigs/MS/odalenetal.BG.2019. Details of the different experiments, plus the command line needed to run each one, are given in readme.txt. A corresponding user manual detailing software installation and configuration, plus cGENIE.muffin model tutorials, is available from, last access: 12 April 2020 and is assigned a DOI (cGENIE release 1.9.1b2018). Datasets are available upon request (e-mail to


The supplement related to this article is available online at:

Author contributions

MÖ, JonN, KICO, and AR designed the model experiments. AR developed the original cGENIE model code. MÖ and KICO adapted the model code and forcings for the experimental design. MÖ performed the model simulations and produced the tables and figures. CDP provided expertise on ocean δ13C analysis. MÖ prepared the manuscript with contributions from all co-authors.

Competing interests

The authors declare that they have no conflict of interest.


The authors thank Pearse James Buchanan and the one anonymous referee for their helpful comments which improved the paper. Malin Ödalen is grateful to Eric Galbraith for helpful discussions in the early stages of this work. Model simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputer Centre (NSC), Sweden. Malin Ödalen would like to acknowledge support from the Bolin Centre for Climate Research, Research Areas 1 and 6.

Financial support

Andy Ridgwell was supported by a Heising-Simons Foundation award and by EU (grant ERC 2013-CoG-617313).

The article processing charges for this open-access
publication were covered by Stockholm University.

Review statement

This paper was edited by Katja Fennel and reviewed by Pearse Buchanan and one anonymous referee.


Adams, J. M. and Faure, H.: A new estimate of changing carbon storage on land since the last glacial maximum, based on global land ecosystem reconstruction, Glob. Planet. Change, 16, 3–24, 1998. a

Adkins, J. F.: The role of deep ocean circulation in setting glacial climates, Paleoceanography, 28, 539–561, 2013. a

Anderson, L. A. and Sarmiento, J. L.: Redfield ratios of remineralization determined by nutrient data analysis, Global Biogeochem. Cy., 8, 65–80, 1994. a, b

Archer, D. and Maier-Reimer, E.: Effect of deep-sea sedimentary calcite preservation on atmospheric CO2 concentration, Nature, 367, 260–263, 1994. a

Archer, D., Winguth, A., Lea, D., and Mahowald, N.: What caused the glacial/interglacial atmospheric pCO2 cycles?, Rev. Geophys., 38, 159–189, 2000. a, b, c

Bereiter, B., Shackleton, S., Baggenstos, D., Kawamura, K., and Severinghaus, J.: Mean global ocean temperatures during the last glacial transition, Nature, 553, 39–44, 2018. a, b, c, d, e

Berger, W.: Deglacial CO2 buildup: constraints on the coral-reef model, Palaeogeogr. Palaeocl., 40, 235–253, 1982. a

Bouttes, N., Paillard, D., and Roche, D. M.: Impact of brine-induced stratification on the glacial carbon cycle, Clim. Past, 6, 575–589,, 2010. a

Bouttes, N., Paillard, D., Roche, D. M., Brovkin, V., and Bopp, L.: Last Glacial Maximum CO2 and δ13C successfully reconciled, Geophys. Res. Lett., 38, L02705,, 2011. a

Bouttes, N., Roche, D. M., and Paillard, D.: Systematic study of the impact of fresh water fluxes on the glacial carbon cycle, Clim. Past, 8, 589–607,, 2012. a

Boyle, E. A. and Keigwin, L.: North Atlantic thermohaline circulation during the past 20,000 years linked to high-latitude surface temperature, Nature, 330, 35–40, 1987. a

Bradtmiller, L., Anderson, R., Sachs, J., and Fleisher, M.: A deeper respired carbon pool in the glacial equatorial Pacific Ocean, Earth Planet. Sc. Lett., 299, 417–425, 2010. a

Broecker, W. S.: Glacial to interglacial changes in ocean chemistry, Prog. Oceanogr., 11, 151–197,, 1982a. a, b

Broecker, W. S.: Ocean chemistry during glacial time, Geochim. Cosmochim. Ac., 46, 1689–1705, 1982b. a

Brovkin, V., Ganopolski, A., Archer, D., and Rahmstorf, S.: Lowering of glacial atmospheric CO2 in response to changes in oceanic circulation and marine biogeochemistry, Paleoceanography, 22, PA4202,, 2007. a

Buchanan, P., Matear, R., Chase, Z., Phipps, S., and Bindoff, N.: Dynamic biological functioning important for simulating and stabilizing ocean biogeochemistry, Global Biogeochem. Cy., 32, 565–593, 2018. a, b, c

cGENIE GitHub repository: available at: (last access: 12 April 2020), 2019. a

cGENIE release 1.9.1b:, available at: (last access: 12 April 2020), 2018. a

cGENIE release v0.9.5:, available at: (last access: 12 April 2020), 2019. a, b

Chikamoto, M., Abe-Ouchi, A., Oka, A., and Smith, S. L.: Temperature-induced marine export production during glacial period, Geophys. Res. Lett., 39, L21601,, 2012. a, b, c, d

Ciais, P., Tagliabue, A., Cuntz, M., Bopp, L., Scholze, M., Hoffmann, G., Lourantou, A., Harrison, S. P., Prentice, I. C., Kelley, D., Koven, C., and Piao, S. L.: Large inert carbon pool in the terrestrial biosphere during the Last Glacial Maximum, Nat. Geosci., 5, 74–79, 2012. a

Ciais, P., Sabine, C., Bala, G., Bopp, L., Brovkin, V., Canadell, J., Chhabra, A., DeFries, R., Galloway, J., Heimann, M., Jones, C., Le Quéré, C., Myeni, R. B., Piao, S., and Thornton, P.: Carbon and other biogeochemical cycles, in: Climate Change 2013: The Physical Science Basis, Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G. K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M., chap. 6, Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 465–570, 2013. a

Crowley, T. J.: Ice age terrestrial carbon changes revisited, Global Biogeochem. Cy., 9, 377–389, 1995. a

Curry, W. B. and Oppo, D. W.: Glacial water mass geometry and the distribution of δ13C of ΣCO2 in the western Atlantic Ocean, Paleoceanography, 20, PA1017,, 2005. a, b

Curry, W. B., Duplessy, J.-C., Labeyrie, L., and Shackleton, N. J.: Changes in the distribution of δ13C of deep water ΣCO2 between the last glaciation and the Holocene, Paleoceanogr. Paleocl., 3, 317–341, 1988. a, b, c

Davies-Barnard, T., Ridgwell, A., Singarayer, J., and Valdes, P.: Quantifying the influence of the terrestrial biosphere on glacial–interglacial climate dynamics, Clim. Past, 13, 1381–1401,, 2017. a, b, c

Duplessy, J., Shackleton, N., Fairbanks, R., Labeyrie, L., Oppo, D., and Kallel, N.: Deepwater source variations during the last climatic cycle and their impact on the global deepwater circulation, Paleoceanography, 3, 343–360, 1988. a, b, c, d, e

Eggleston, S. and Galbraith, E. D.: The devil's in the disequilibrium: multi-component analysis of dissolved carbon and oxygen changes under a broad range of forcings in a general circulation model, Biogeosciences, 15, 3761–3777,, 2018. a

Eppley, R. W.: Temperature and phytoplankton growth in the sea, Fish. Bull, 70, 1063–1085, 1972. a, b

Galbraith, E. and de Lavergne, C.: Response of a comprehensive climate model to a broad range of external forcings: relevance for deep ocean ventilation and the development of late Cenozoic ice ages, Clim. Dynam., 52, 653–679, 2018. a, b, c, d

Galbraith, E. D. and Jaccard, S. L.: Deglacial weakening of the oceanic soft tissue pump: global constraints from sedimentary nitrogen isotopes and oxygenation proxies, Quaternary Sci. Rev., 109, 38–48, 2015. a, b

Galbraith, E. D. and Martiny, A. C.: A simple nutrient-dependence mechanism for predicting the stoichiometry of marine ecosystems, P. Natl. Acad. Sci. USA, 112, 8199–8204, 2015. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o

Ganopolski, A. and Brovkin, V.: Simulation of climate, ice sheets and CO2 evolution during the last four glacial cycles with an Earth system model of intermediate complexity, Clim. Past, 13, 1695–1716,, 2017. a

Garcia, C. A., Baer, S. E., Garcia, N. S., Rauschenberg, S., Twining, B. S., Lomas, M. W., and Martiny, A. C.: Nutrient supply controls particulate elemental concentrations and ratios in the low latitude eastern Indian Ocean, Nat. Commun., 9, 4868,, 2018a. a

Garcia, H. E., Weathers, K., Paver, C. R., Smolyar, I., Boyer, T. P., Locarnini, R. A., Zweng, M. M., Mishonov, A. V., Baranova, O., Seidov, D., and Reagan, J. R.: World Ocean Atlas 2018, Volume 4: Dissolved Inorganic Nutrients (phosphate, nitrate, silicate), Tech. Rep. 84, NOAA Atlas NESDIS, 35 pp., 2018b. a, b

Garcia, H. E., Weathers, K., Paver, C. R., Smolyar, I., Boyer, T. P., Locarnini, R. A., Zweng, M. M., Mishonov, A. V., Baranova, O. K., Seidov, D., and Reagan, J. R.: World Ocean Atlas 2018, Volume 3: Dissolved Oxygen, Apparent Oxygen Utilization, and Oxygen Saturation), Tech. Rep. 82, NOAA Atlas NESDIS, 38 pp., 2018c. a, b, c

Gebbie, G., Peterson, C. D., Lisiecki, L. E., and Spero, H. J.: Global-mean marine δ13C and its uncertainty in a glacial state estimate, Quaternary Sci. Rev., 125, 144–159, 2015. a, b

Griffies, S. M.: The Gent–McWilliams skew flux, J. Phys. Oceanogr., 28, 831–841, 1998. a

Hain, M. P., Sigman, D. M., and Haug, G. H.: Carbon dioxide effects of Antarctic stratification, North Atlantic Intermediate Water formation, and subantarctic nutrient drawdown during the last ice age: Diagnosis and synthesis in a geochemical box model, Global Biogeochem. Cy., 24, GB4023,, 2010. a

Headly, M. A. and Severinghaus, J. P.: A method to measure Kr∕N2 ratios in air bubbles trapped in ice cores and its application in reconstructing past mean ocean temperature, J. Geophys. Res.-Atmos., 112, D19105,, 2007. a, b, c, d

Herguera, J., Herbert, T., Kashgarian, M., and Charles, C.: Intermediate and deep water mass distribution in the Pacific during the Last Glacial Maximum inferred from oxygen and carbon stable isotopes, Quaternary Sci. Rev., 29, 1228–1245, 2010. a

Hesse, T., Butzin, M., Bickert, T., and Lohmann, G.: A model-data comparison of δ13C in the glacial Atlantic Ocean, Paleoceanography, 26, PA3220,, 2011. a

Heuzé, C., Heywood, K. J., Stevens, D. P., and Ridley, J. K.: Southern Ocean bottom water characteristics in CMIP5 models, Geophys. Res. Lett., 40, 1409–1414, 2013. a

Hewitt, C. D., Broccoli, A., Crucifix, M., Gregory, J., Mitchell, J., and Stouffer, R.: The effect of a large freshwater perturbation on the glacial North Atlantic Ocean using a coupled general circulation model, J. Clim., 19, 4436–4447, 2006. a

Ito, T. and Follows, M. J.: Preformed phosphate, soft tissue pump and atmospheric CO2, J. Mar. Res., 63, 813–839, 2005. a, b

Jaccard, S. L. and Galbraith, E. D.: Large climate-driven changes of oceanic oxygen concentrations during the last deglaciation, Nat. Geosci., 5, 151–159, 2012. a

Kohfeld, K., Graham, R., De Boer, A., Sime, L., Wolff, E., Le Quéré, C., and Bopp, L.: Southern Hemisphere westerly wind changes during the Last Glacial Maximum: paleo-data synthesis, Quaternary Sci. Rev., 68, 76–95, 2013. a

Kohfeld, K. E. and Ridgwell, A.: Glacial-Interglacial Variability in Atmospheric CO2, in: Surface Ocean-Lower Atmosphere Processes, edited by: Le Quéré, C. and S., S. E., American Geophysical Union, Washington, DC, 251–286, 2009. a, b, c

Kohfeld, K. E., Le Quéré, C., Harrison, S. P., and Anderson, R. F.: Role of marine biology in glacial-interglacial CO2 cycles, Science, 308, 74–78, 2005. a

Kolowith, L. C., Ingall, E. D., and Benner, R.: Composition and cycling of marine organic phosphorus, Limnol. Oceanogr., 46, 309–320, 2001. a

Kumar, N., Anderson, R., Mortlock, R., Froelich, P., Kubik, P., Dittrich-Hannen, B., and Suter, M.: Increased biological productivity and export production in the glacial Southern Ocean, Nature, 378, 675–680, 1995. a

Kwon, E. Y., Primeau, F., and Sarmiento, J. L.: The impact of remineralization depth on the air–sea carbon balance, Nat. Geosci., 2, 630–635, 2009. a

Lauderdale, J. M., Garabato, A. C. N., Oliver, K. I. C., Follows, M. J., and Williams, R. G.: Wind-driven changes in Southern Ocean residual circulation, ocean carbon reservoirs and atmospheric CO2, Clim. Dynam., 41, 2145–2164, 2013. a, b, c, d

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

Le Quéré, C., Harrison, S. P., Colin Prentice, I., Buitenhuis, E. T., Aumont, O., Bopp, L., Claustre, H., Cotrim Da Cunha, L., Geider, R., Giraud, X., Klaas, C., Kohfeld, K. E., Legendre, L., Manizza, M., Platt, T., Rivkin, R. B., Sathyendranath, S., Uitz, J., Watson, A. J., and Wolf-Gladrow, D.: Ecosystem dynamics based on plankton functional types for global ocean biogeochemistry models, Glob. Change Biol., 11, 2016–2040, 2005. a, b

Letscher, R. T. and Moore, J. K.: Preferential remineralization of dissolved organic phosphorus and non-Redfield DOM dynamics in the global ocean: Impacts on marine productivity, nitrogen fixation, and carbon export, Global Biogeochem. Cy., 29, 325–340, 2015. a

Locarnini, R. A., Mishonov, A. V., Baranova, O. K., Boyer, T. P., Zweng, M. M., Garcia, H. E., Reagan, J. R., Seidov, D., Weathers, K., Paver, C. R., and Smolyar, I.: World Ocean Atlas 2018, Volume 1: Temperature, Tech. Rep. 81, NOAA Atlas NESDIS, 52 pp., 2018. a, b, c

Lüthi, D., Le Floch, M., Bereiter, B., Blunier, T., Barnola, J.-M., Siegenthaler, U., Raynaud, D., Jouzel, J., Fischer, H., Kawamura, K., and Stocker, T. F.: High-resolution carbon dioxide concentration record 650,000–800,000 years before present, Nature, 453, 379–382,, 2008. a

Lynch-Stieglitz, J., Curry, W. B., Slowey, N., and Schmidt, G. A.: The overturning circulation of the glacial Atlantic, in: Reconstructing Ocean History, Springer, 7–31, 1999. a

Mahowald, N. M., Muhs, D. R., Levis, S., Rasch, P. J., Yoshioka, M., Zender, C. S., and Luo, C.: Change in atmospheric mineral aerosols in response to climate: Last glacial period, preindustrial, modern, and doubled carbon dioxide climates, J. Geophys. Res.-Atmos., 111, D10202,, 2006. a, b, c, d

Marchal, O., Stocker, T. F., and Joos, F.: A latitude-depth, circulation-biogeochemical ocean model for paleoclimate studies, Development and sensitivities, Tellus B, 50, 290–316, 1998. a

Marchitto, T. M. and Broecker, W. S.: Deep water mass geometry in the glacial Atlantic Ocean: A review of constraints from the paleonutrient proxy Cd∕Ca, Geochem. Geophy. Geosy., 7, Q12003,, 2006. a, b

Marchitto, T. M., Lehman, S. J., Ortiz, J. D., Flückiger, J., and van Geen, A.: Marine radiocarbon evidence for the mechanism of deglacial atmospheric CO2 rise, Science, 316, 1456–1459, 2007. a

Marinov, I., Gnanadesikan, A., Toggweiler, J., and Sarmiento, J.: The southern ocean biogeochemical divide, Nature, 441, 964–967, 2006. a, b

Marinov, I., Gnanadesikan, A., Sarmiento, J. L., Toggweiler, J. R., Follows, M., and Mignone, B. K.: Impact of oceanic circulation on biological carbon storage in the ocean and atmospheric pCO2, Global Biogeochem. Cy., 22, GB3007,, 2008. a

Marsh, R., Müller, S. A., Yool, A., and Edwards, N. R.: Incorporation of the C-GOLDSTEIN efficient climate model into the GENIE framework: “eb_go_gs” configurations of GENIE, Geosci. Model Dev., 4, 957–992,, 2011. a

Martin, J. H.: Glacial-interglacial CO2 change: The iron hypothesis, Paleoceanography, 5, 1–13, 1990. a, b, c, d, e, f

Matsumoto, K.: Biology-mediated temperature control on atmospheric pCO2 and ocean biogeochemistry, Geophys. Res. Lett., 34, L20605,, 2007. a, b, c, d, e, f, g, h, i

Matsumoto, K., Oba, T., Lynch-Stieglitz, J., and Yamamoto, H.: Interior hydrography and circulation of the glacial Pacific Ocean, Quaternary Sci. Rev., 21, 1693–1704, 2002. a

Mayr, C., Lücke, A., Wagner, S., Wissel, H., Ohlendorf, C., Haberzettl, T., Oehlerich, M., Schäbitz, F., Wille, M., Zhu, J., and Zolitschka, B.: Intensified Southern Hemisphere Westerlies regulated atmospheric CO2 during the last deglaciation, Geology, 41, 831–834, 2013. a

McCarthy, G., Smeed, D., Johns, W. E., Frajka-Williams, E., Moat, B., Rayner, D., Baringer, M., Meinen, C., Collins, J., and Bryden, H.: Measuring the Atlantic meridional overturning circulation at 26 N, Prog. Oceanogr., 130, 91–111, 2015. a, b, c

McInerney, F. A. and Wing, S. L.: The Paleocene-Eocene Thermal Maximum: A perturbation of carbon cycle, climate, and biosphere with implications for the future, Ann. Rev. Earth Planet. Sc., 39, 489–516, 2011. a

Menviel, L., Joos, F., and Ritz, S.: Simulating atmospheric CO2, 13C 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 Sci. Rev., 56, 46–68, 2012. a, b, c

Menviel, L., Yu, J., Joos, F., Mouchet, A., Meissner, K., and England, M.: Poorly ventilated deep ocean at the Last Glacial Maximum inferred from carbon isotopes: A data-model comparison study, Paleoceanogr. Paleocl., 32, 2–17, 2017. a, b, c, d, e, f, g

Mook, W.: 13C in atmospheric CO2, Neth. J. Sea Res., 20, 211–223, 1986. a

Moore, C., Mills, M., Arrigo, K., Berman-Frank, I., Bopp, L., Boyd, P., Galbraith, E., Geider, R., Guieu, C., Jaccard, S. L., Jickells, T. D., La Roche, J., Lenton, T. M., Mahowald, N. M., Maranon, E., Marinov, I., Moore, J. K., Nakatsuka, T., Oschlies, A., Saito, M. A., Thingstad, T. F., Tsuda, A., and Ulloa, O.: Processes and patterns of oceanic nutrient limitation, Nat. Geosci., 6, 701–710, 2013. a, b

Moreno, A. R., Hagstrom, G. I., Primeau, F. W., Levin, S. A., and Martiny, A. C.: Marine phytoplankton stoichiometry mediates nonlinear interactions between nutrient supply, temperature, and atmospheric CO2, Biogeosciences, 15, 2761–2779,, 2018. a, b

Muglia, J., Skinner, L. C., and Schmittner, A.: Weak overturning circulation and high Southern Ocean nutrient utilization maximized glacial ocean carbon, Earth Planet. Sc. Lett., 496, 47–56, 2018. a

Ödalen, M., Nycander, J., Oliver, K. I. C., Brodeau, L., and Ridgwell, A.: The influence of the ocean circulation state on ocean carbon storage and CO2 drawdown potential in an Earth system model, Biogeosciences, 15, 1367–1393,, 2018. a, b, c, d, e, f, g, h, i

Paulmier, A., Kriest, I., and Oschlies, A.: Stoichiometries of remineralisation and denitrification in global biogeochemical ocean models, Biogeosciences, 6, 923–935,, 2009. a

Peterson, C. D. and Lisiecki, L. E.: Deglacial carbon cycle changes observed in a compilation of 127 benthic δ13C time series (20–6 ka), Clim. Past, 14, 1229–1252,, 2018. a, b

Peterson, C. D., Lisiecki, L. E., and Stern, J. V.: Deglacial whole-ocean δ13C change estimated from 480 benthic foraminiferal records, Paleoceanography, 29, 549–563, 2014. a, b, c, d, e

Petit, J.-R., Jouzel, J., Raynaud, D., Barkov, N. I., Barnola, J.-M., Basile, I., Bender, M., Chappellaz, J., Davis, M., Delaygue, G., Delmotte, M., Kotlyakov, V. M., Legrand, M., Lipenkov, V. Y., Lorius, C., Pépin, L., Ritz, C., Saltzman, E., and Stievenard, M.: Climate and atmospheric history of the past 420,000 years from the Vostok ice core, Antarctica, Nature, 399, 429–436, 1999. a

Rau, G. H., Riebesell, U., and Wolf-Gladrow, D.: A model of photosynthetic 13C fractionation by marine phytoplankton based on diffusive molecular CO2 uptake, Mar. Ecol. Prog. Ser., 133, 275–285, 1996. a, b, c, d

Rau, G. H., Riebesell, U., and Wolf-Gladrow, D.: CO2aq-dependent photosynthetic 13C fractionation in the ocean: A model versus measurements, Global Biogeochem. Cy., 11, 267–278, 1997. a

Redfield, A. C.: The influence of organisms on the composition of sea–water, The Sea, 2, 26–77, 1963. a, b, c, d, e

Ridgwell, A., Hargreaves, J. C., Edwards, N. R., Annan, J. D., Lenton, T. M., Marsh, R., Yool, A., and Watson, A.: Marine geochemical data assimilation in an efficient Earth System Model of global biogeochemical cycling, Biogeosciences, 4, 87–104,, 2007. a, b, c, d, e, f, g

Ridgwell, A. J.: Glacial-interglacial perturbations in the global carbon cycle, Ph.D. thesis, University of East Anglia, 146 pp., 2001. a, b, c, d

Sarmiento, J. and Toggweiler, J.: A new model for the role of the oceans in determining atmospheric pCO2, Nature, 308, 621–624, 1984. a, b, c

Schmittner, A., Meissner, K., Eby, M., and Weaver, A.: Forcing of the deep ocean circulation in simulations of the Last Glacial Maximum, Paleoceanography, 17, 1015,, 2002. a

Shackleton, N.: Carbon-13 in Uvigerina: Tropical rain forest history and the equatorial Pacific carbonate dissolution cycle, in: The fate of fossil fuel CO2 in the oceans, Plenum Press, 401–428, 1977. a, b

Sigman, D. M. and Boyle, E. A.: Glacial/interglacial variations in atmospheric carbon dioxide, Nature, 407, 859–869, 2000. a, b, c

Sigman, D. M., Hain, M. P., and Haug, G. H.: The polar ocean and glacial cycles in atmospheric CO2 concentration, Nature, 466, 47–55, 2010. a

Sime, L. C., Kohfeld, K. E., Le Quéré, C., Wolff, E. W., de Boer, A. M., Graham, R. M., and Bopp, L.: Southern Hemisphere westerly wind changes during the Last Glacial Maximum: model-data comparison, Quaternary Sci. Rev., 64, 104–120, 2013. a

Sime, L. C., Hodgson, D., Bracegirdle, T. J., Allen, C., Perren, B., Roberts, S., and de Boer, A. M.: Sea ice led to poleward-shifted winds at the Last Glacial Maximum: the influence of state dependency on CMIP5 and PMIP3 models, Clim. Past, 12, 2241–2253,, 2016. a

Skinner, L., Fallon, S., Waelbroeck, C., Michel, E., and Barker, S.: Ventilation of the deep Southern Ocean and deglacial CO2 rise, Science, 328, 1147–1151, 2010. a

Skinner, L., Primeau, F., Freeman, E., de la Fuente, M., Goodwin, P., Gottschalk, J., Huang, E., McCave, I., Noble, T., and Scrivner, A.: Radiocarbon constraints on the glacial ocean circulation and its impact on atmospheric CO2, Nat. Commun., 8, 16010,, 2017. a, b, c, d

Stephens, B. B. and Keeling, R. F.: The influence of Antarctic sea ice on glacial–interglacial CO2 variations, Nature, 404, 171–174, 2000. a

Stocker, T.: Climate change 2013: the physical science basis: Working Group I contribution to the Fifth assessment report of the Intergovernmental Panel on Climate Change, Cambridge University Press, 2014. a

Tagliabue, A., Aumont, O., DeAth, R., Dunne, J. P., Dutkiewicz, S., Galbraith, E., Misumi, K., Moore, J. K., Ridgwell, A., Sherman, E., Stock, C., Vichi, M., Völker, C., and Yool, A.: How well do global ocean biogeochemistry models simulate dissolved iron distributions?, Global Biogeochem. Cy., 30, 149–174, 2016. a, b, c

Takahashi, T., Broecker, W. S., and Langer, S.: Redfield ratio based on chemical data from isopycnal surfaces, J. Geophys. Res.-Ocean., 90, 6907–6924, 1985. a, b

Tanioka, T. and Matsumoto, K.: Buffering of ocean export production by flexible elemental stoichiometry of particulate organic matter, Global Biogeochem. Cy., 31, 1528–1542, 2017. a, b, c, d, e, f

Turner, S. K. and Ridgwell, A.: Development of a novel empirical framework for interpreting geological carbon isotope excursions, with implications for the rate of carbon injection across the PETM, Earth Planet. Sc. Lett., 435, 1–13, 2016. a

Walin, G., Hieronymus, J., and Nycander, J.: Source-related variables for the description of the oceanic carbon system, Geochem. Geophy. Geosy., 15, 3675–3687, 2014. a

Wanninkhof, R.: Relationship between gas exchange and wind speed over the ocean, J. Geophys. Res., 97, 7373–7382, 1992. a

Watson, A. J., Bakker, D., Ridgwell, A., Boyd, P., and Law, C.: Effect of iron supply on Southern Ocean CO2 uptake and implications for glacial atmospheric CO2, Nature, 407, 730–733, 2000. a

Yvon-Durocher, G., Dossena, M., Trimmer, M., Woodward, G., and Allen, A. P.: Temperature and the biogeography of algal stoichiometry, Glob. Ecol. Biogeogr., 24, 562–570, 2015. a, b

Zeebe, R. E. and Wolf-Gladrow, D. A.: CO2 in seawater: equilibrium, kinetics, isotopes, 65, Gulf Professional Publishing, 2001. a

Zhang, J., Quay, P., and Wilbur, D.: Carbon isotope fractionation during gas-water exchange and dissolution of CO2, Geochim. Cosmochim. Ac., 59, 107–114, 1995. a, b

Short summary
In glacial periods, ocean uptake of carbon is likely a key player for achieving low atmospheric CO2. In climate models, ocean biological uptake of carbon (C) and phosphorus (P) are often assumed to occur in fixed proportions. In this study, we allow the ratio of C : P to vary and simulate, to first approximation, the complex biological changes that occur in the ocean over long timescales. We show here that, for glacial–interglacial cycles, this complexity contributes to low atmospheric CO2.
Final-revised paper