Articles | Volume 17, issue 1
Research article
15 Jan 2021
Research article |  | 15 Jan 2021

Sequential changes in ocean circulation and biological export productivity during the last glacial–interglacial cycle: a model–data study

Cameron M. O'Neill, Andrew McC. Hogg, Michael J. Ellwood, Bradley N. Opdyke, and Stephen M. Eggins

We conduct a model–data analysis of the marine carbon cycle to understand and quantify the drivers of atmospheric CO2 concentration during the last glacial–interglacial cycle. We use a carbon cycle box model, “SCP-M”, combined with multiple proxy data for the atmosphere and ocean, to test for variations in ocean circulation and Southern Ocean biological export productivity across marine isotope stages spanning 130 000 years ago to the present. The model is constrained by proxy data associated with a range of environmental conditions including sea surface temperature, salinity, ocean volume, sea-ice cover and shallow-water carbonate production. Model parameters for global ocean circulation, Atlantic meridional overturning circulation and Southern Ocean biological export productivity are optimized in each marine isotope stage against proxy data for atmospheric CO2, δ13C and Δ14C and deep-ocean δ13C, Δ14C and CO32-. Our model–data results suggest that global overturning circulation weakened during Marine Isotope Stage 5d, coincident with a  25 ppm fall in atmospheric CO2 from the last interglacial period. There was a transient slowdown in Atlantic meridional overturning circulation during Marine Isotope Stage 5b, followed by a more pronounced slowdown and enhanced Southern Ocean biological export productivity during Marine Isotope Stage 4 (−30 ppm). In this model, the Last Glacial Maximum was characterized by relatively weak global ocean and Atlantic meridional overturning circulation and increased Southern Ocean biological export productivity (−20 ppm during MIS 3 and MIS 2). Ocean circulation and Southern Ocean biological export productivity returned to modern values by the Holocene period. The terrestrial biosphere decreased by 385 Pg C in the lead-up to the Last Glacial Maximum, followed by a period of intense regrowth during the last glacial termination and the Holocene ( 600 Pg C). Slowing ocean circulation, a colder ocean and to a lesser extent shallow carbonate dissolution contributed −70 ppm to atmospheric CO2 in the  100 000-year lead-up to the Last Glacial Maximum, with a further −15 ppm contributed during the glacial maximum. Our model results also suggest that an increase in Southern Ocean biological export productivity was one of the ingredients required to achieve the Last Glacial Maximum atmospheric CO2 level. We find that the incorporation of glacial–interglacial proxy data into a simple quantitative ocean transport model provides useful insights into the timing of past changes in ocean processes, enhancing our understanding of the carbon cycle during the last glacial–interglacial period.

1 Introduction

Large and regular fluctuations in the concentration of atmospheric CO2 and ocean proxy signals for carbon isotopes and carbonate ion concentration during the last 800 kyr are preserved in ice and marine core records. The most obvious of these fluctuations is the repeated oscillation of atmospheric CO2 concentration over the range of  180–280 ppm every  100 kyr. The magnitude and regularity of these oscillations in atmospheric CO2, combined with proxy observations for carbon isotopes, point to the quasi-regular transfer of carbon between the main Earth reservoirs: the ocean, atmosphere, terrestrial biosphere and marine sediments (Broecker1982; Sigman and Boyle2000; Toggweiler2008; Hogg2008; Kohfeld and Ridgwell2009; Menviel et al.2012; Kohfeld and Chase2017; Ganopolski and Brovkin2017). The ocean, given its large size as a carbon store and ongoing exchange of CO2 with the atmosphere, likely plays the key role in changing atmospheric CO2 (Broecker1982; Knox and McElroy1984; Siegenthaler and Wenk1984; Sarmiento and Toggweiler1984; Sigman and Boyle2000; Kohfeld and Ridgwell2009). Ocean-centric hypotheses for variation in atmospheric CO2 concentration have been examined in great detail for the Last Glacial Maximum (LGM) and Holocene periods, supported by the abundance of paleo data from marine sediment coring and sampling activity (e.g. Sikes et al.2000; Curry and Oppo2005; Kohfeld and Ridgwell2009; Oliver et al.2010; Menviel et al.2012; Peterson et al.2014; Yu et al.2014b; Menviel et al.2016; Skinner et al.2017; Muglia et al.2018; Yu et al.2019). However, the hypotheses for variation in atmospheric CO2 across the LGM–Holocene remain debated (e.g. Kohfeld et al.2005; Martinez-Garcia et al.2014; Menviel et al.2016; Skinner et al.2017; Muglia et al.2018; Khatiwala et al.2019). Established hypotheses include those emphasizing ocean biology (e.g. Martin1990; Martinez-Garcia et al.2014), ocean circulation (e.g. Burke and Robinson2012; Menviel et al.2016; Skinner et al.2017), sea surface temperature (SST) (Khatiwala et al.2019), or the aggregate effect of several mechanisms (e.g. Kohfeld and Ridgwell2009; Hain et al.2010; Köhler et al.2010; Menviel et al.2012; Ferrari et al.2014; Ganopolski and Brovkin2017; Muglia et al.2018) to explain the LGM–Holocene carbon cycle transition. Hypotheses for an ocean biological role include the effects of iron fertilization on biological export productivity (e.g. Martin1990; Watson et al.2000; Martinez-Garcia et al.2014), the depth of remineralization of particulate organic carbon (POC) (e.g. Matsumoto2007; Kwon et al.2009; Menviel et al.2012), changes in the organic carbon : carbonate (“the rain ratio”) or carbon : silicate constitution of marine organisms (e.g. Archer and Maier-Reimer1994; Harrison2000), and increased biological utilization of exposed shelf-derived nutrients such as phosphorus (e.g. Menviel et al.2012).

Several studies have attempted to solve the problem of glacial–interglacial CO2 by modelling either the last glacial–interglacial cycle in its entirety or multiple glacial–interglacial cycles (e.g. Ganopolski et al.2010; Menviel et al.2012; Brovkin et al.2012; Ganopolski and Brovkin2017). These studies highlight the roles of orbitally forced Northern Hemisphere ice sheets in the onset of the glacial periods and important feedbacks from ocean circulation, carbonate chemistry and marine biological productivity throughout the glacial cycle (Ganopolski et al.2010; Brovkin et al.2012; Ganopolski and Brovkin2017). Menviel et al. (2012) modelled a range of physical, biological and biogeochemical mechanisms to deliver the full amplitude of atmospheric CO2 variation in the last glacial–interglacial cycle, using transient simulations with the Bern3D model. According to Brovkin et al. (2012), a  50 ppm drop in atmospheric CO2 concentration early in the last glacial–interglacial cycle was caused by lower SST, increased Northern Hemisphere ice sheet cover and the expansion of southern-sourced abyssal waters in place of North Atlantic Deep Water (NADW) formation. Ganopolski and Brovkin (2017) modelled the last four glacial–interglacial cycles with orbital forcing as the singular driver of carbon cycle feedbacks. They described the “carbon stew”, a feedback of combined physical and biogeochemical changes in the carbon cycle driving the last four glacial–interglacial cycles of atmospheric CO2.

Kohfeld and Chase (2017) also extended the LGM–Holocene CO2 debate further into the past by evaluating proxy data over the period 115 000–18 000 years before present (115–18 ka), a time that encompasses the gradual fall in atmospheric CO2 of  85–90 ppm from the last interglacial period until the last glacial termination. Kohfeld and Chase (2017) identified time periods during which CO2 decreased and aligned these with concomitant changes in proxies for SST, sea-ice extent, deep Atlantic Ocean circulation, and mixing and ocean biological productivity. Kohfeld and Chase (2017) observed that the  100 kyr transition to the LGM involved three discrete CO2 reduction events. Firstly, a drop in atmospheric CO2 of  35 ppm at  115–100 ka (Marine Isotope Stage, or MIS, 5d) was accompanied by lower SST and the expansion of Antarctic sea-ice cover. A second phase of CO2 drawdown between 72 and 65 ka (MIS 4), of  40 ppm, likely resulted from a slowdown in deep-ocean circulation (Kohfeld and Chase2017). Finally, during the period 40–18 ka (MIS 3-2) atmospheric CO2 dropped a further 5–10 ppm, which according to Kohfeld and Chase (2017) was the result of enhanced Southern Ocean biological productivity, continually intensifying deep-ocean stratification, shoaling of NADW and northward extension of Antarctic Bottom Water (AABW).

In this paper we quantitatively test the Kohfeld and Chase (2017) hypothesis by undertaking model–data experiments in each MIS across the last glacial–interglacial cycle. We extend their analysis to include Pacific and Indian Ocean modelling and proxy data. We use the SST reconstructions compiled by Kohfeld and Chase (2017) and other proxy records presented in Kohfeld and Chase (2017), covering the last glacial–interglacial cycle. We apply a carbon cycle box model (O'Neill et al.2019) constrained by available atmospheric and oceanic proxy data, to solve for optimal model–data parameter solutions for ocean circulation and biological export productivity. We also present a qualitative analysis of the compiled proxy data to place the model–data experiment results in context. We thereby further constrain the timing and magnitude of posited CO2 mechanisms operating during each MIS in the last glacial–interglacial cycle (e.g. Kohfeld and Ridgwell2009; Oliver et al.2010; Menviel et al.2012; Brovkin et al.2012; Yu et al.2013; Eggleston et al.2016; Yu et al.2016; Kohfeld and Chase2017). This longer-dated analysis complements recent multi-proxy model–data studies of the LGM and Holocene (e.g. Menviel et al.2016; Kurahashi-Nakamura et al.2017; Muglia et al.2018; O'Neill et al.2019) by testing for changes in the ocean carbon cycle in the lead-up to the LGM, in addition to the LGM to Holocene. Our modelling approach differs from other model studies of the last glacial–interglacial cycle (e.g. Ganopolski et al.2010; Menviel et al.2012; Brovkin et al.2012; Ganopolski and Brovkin2017) because we constrain several physical processes from observations (SST, sea level, sea-ice cover, salinity, coral reef fluxes of carbon) and then solve for the values of model parameters for ocean circulation and biology based on an optimization against atmospheric and ocean proxy data.

2 Materials and methods

2.1 Model description

We use the Simple Carbon Project Model (SCP-M) carbon cycle box model in our model–data experiment (O'Neill et al.2019). In summary, SCP-M contains simple parameterizations of the major fluxes in the Earth's surface carbon cycle (Fig. 1). SCP-M incorporates the ocean, atmosphere, terrestrial biosphere and marine/continental sediment carbon reservoirs, weathering and river fluxes, and a number of variables including atmospheric CO2, dissolved inorganic carbon (DIC), phosphorus, alkalinity, carbon isotopes (13C and 14C) and CO32-. SCP-M calculates ocean pCO2 using the equations of Follows et al. (2006) and applies the first and second “dissociation constants” of carbonic acid estimated by Lueker et al. (2000) to calculate HCO3- and CO32- concentrations, respectively, in units of µmol kg−1, in each ocean box. The model employs partial differential equations for determining the concentration of elements, with each box represented as a row and column in a matrix. In this paper, we extend SCP-M by incorporating a separate basin for the combined Pacific and Indian oceans (Fig. 1) following the conceptual model of Talley (2013), to incorporate modelling and proxy data for those regions of the ocean. This version of SCP-M consists of 12 ocean boxes plus the atmosphere and terrestrial biosphere. SCP-M splits out depth regions of the ocean between surface boxes (100–250 m average depth) and intermediate (1000 m average depth), deep (2500 m average depth) and abyssal depth boxes (3700 (Atlantic)–4000 m (Pacific–Indian) average depth). The Southern Ocean is split into two boxes, including a polar box which covers the latitude range 60–80 S (box 12 in Fig. 1), and subpolar Southern Ocean boxes in the Atlantic (box 7) and Pacific–Indian (box 11) basins, which cover the latitude range 40–60 S. See O'Neill et al. (2019) for a discussion of the choice of box depth and latitude dimensions.

Figure 1SCP-M configured as a 12-box ocean model plus atmosphere with marine sediments, continents and the terrestrial biosphere. Exchange of elemental concentrations occur due to fluxes between boxes. Ψ1 (red arrows) is global overturning circulation (GOC), Ψ2 (orange arrows) is Atlantic meridional overturning circulation (AMOC). GOC upwelling in both basins is set by default to a 50 % split between upwelling into the subpolar and polar Southern Ocean. Ψ3 (pink arrows) is Antarctic intermediate water (AAIW) and Subantarctic mode water (SAMW) formation in the Indian and Pacific oceans (e.g. Talley2013). Blue arrows represent mixing fluxes between boxes. γ1 and γ3 parameterize deep–abyssal and Southern Ocean deep topographically induced mixing (e.g. De Boer and Hogg2014), while γ2 is low-latitude thermohaline mixing (e.g. Liu et al.2016). Z (green downward arrows) is the biological pump, FCA (white downward arrows) is the carbonate pump, DCA (white squiggles) is carbonate dissolution and P (black, bidirectional arrows) is the air–sea gas exchange. Key to boxes: Atlantic (box 1: low-latitude or tropical surface ocean, 0–100 m; box 2: northern surface ocean, 0-250 m; box 3: intermediate ocean, 100–1000 m; box 4: deep ocean, 1000–2500 m; box 6: abyssal ocean, 2500–3700 m; box 7: subpolar southern surface ocean, 0–250 m). Pacific–Indian (box 8: low-latitude or tropical surface ocean, 0–100 m; box 9: deep ocean, 100–2500 m; box 10: abyssal ocean, 2500–4000 m; box 11: subpolar southern surface ocean, 0–250 m). Southern Ocean (box 5: intermediate-deep; box 12: surface ocean). For a more detailed model description, see O'Neill et al. (2019) and updated model code and data at


The major ocean carbon flux parameters of interest in this model–data study are global ocean circulation (GOC), Ψ1, Atlantic meridional overturning circulation (AMOC), Ψ2, and ocean biological export productivity, Z. The ocean circulation parameters Ψ1 and Ψ2 are simply prescribed in units of Sverdrups (Sv, 106 m3 s−1). Ocean biological export productivity Z is calculated using the method of Martin et al. (1987). The biological productivity flux at 100 m depth is attenuated with depth for each box according to the decay rule of Martin et al. (1987). Each subsurface box receives a biological flux of an element at its ceiling depth and loses a flux at its floor depth (lost to the boxes below it). The difference between influx and outflux is the amount of element that is remineralized into each box. The input parameter is the value of export production at 100 m depth, in units of mol C m−2 yr−1 as per Martin et al. (1987). Equation (1) shows the general form of the Martin et al. (1987) equation:

(1) F = F 100 ( d 100 ) b ,

where F is a flux of carbon in mol C m−2 yr−1, F100 is an estimate of carbon flux at 100 m depth, d is depth in metres and b is a depth scalar. In SCP-M, the Z parameter implements the Martin et al. (1987) equation. Z is an estimate of biological productivity at 100 m depth (in mol C m−2 yr−1), and coupled with the Martin et al. (1987) depth scalar, it controls the amount of organic carbon that sinks from each model surface box to the boxes below.

Air–sea gas exchange is based on the relative pCO2 between the surface ocean boxes and the atmosphere and is implemented in SCP-M by a parameter that sets its rate in metres per day: P (Fig. 1). SCP-M parameterizes shallow-water carbonate production, which is linked to the Z parameter by an assumption for the relative proportion of carbonate vs. organic matter in the biological export flux, known as the rain ratio (e.g. Archer and Maier-Reimer1994; Ridgwell2003). Carbonate dissolution is calculated based on the ocean box or marine surface sediment calcium carbonate concentration relative to a depth-dependant saturation concentration (Morse and Berner1972; Millero1983). The isotopes of carbon are calculated applying various fractionation factors associated with the biological, physical and chemical fluxes of carbon (see Table S1 and O'Neill et al.2019).

We have added a simple representation of shallow-water carbonate fluxes of carbon and alkalinity in SCP-M's low-latitude surface boxes, to cater for this feature in theories for glacial–interglacial cycle CO2 (e.g. Berger1982; Opdyke and Walker1992; Ridgwell et al.2003; Vecsei and Berger2004; Menviel and Joos2012), using

(2) d C i d t reef = C reef / V i ,

where Creef is the prescribed flux of carbon out of/into the low-latitude surface ocean boxes during net reef accumulation or dissolution, in mol C yr−1, and Vi is the volume of the low-latitude surface box i. The alkalinity flux associated with reef production or dissolution is simply Eq. (2) multiplied by 2 (e.g. Sarmiento and Gruber2006).

SCP-M contains a simple parameterization of the terrestrial carbon cycle. For continental rock weathering, we apply the simple scheme of Walker and Kasting (1992) as implemented in Toggweiler (2008), Hogg (2008) and Zeebe (2012). Weathering of silicate and carbonate rocks supplies DIC and alkalinity to the low-latitude surface ocean boxes in each basin (boxes 1 and 8 in Fig. 1) as a function of a weathering constant and atmospheric CO2, in units of mol m−3 yr−1. The parameter values used are shown in Table S1. For the SCP-M weathering equations please see O'Neill et al. (2019). δ13C fluxes for carbonate and silicate weathering are shown in Table S1. A volcanic flux of carbon (and δ13C) is also assumed, which sets the rate of volcanic CO2 outgassing roughly to the rate of silicate rock weathering (Walker and Kasting1992; Toggweiler2008; Hogg2008; Zeebe2012). Parameters for volcanic CO2 and δ13C fluxes are shown in Table S1.

The terrestrial biosphere is represented in SCP-M as a stock of carbon (a box) that fluxes with the atmosphere, governed by parameters for net primary productivity (NPP) and respiration. In SCP-M, NPP is calculated as a function of carbon fertilization, which increases NPP as atmospheric CO2 rises via a simple logarithmic relationship, using the model of Harman et al. (2011). This is a simplified approach, which omits the effects of temperature and precipitation on NPP (François et al.1999; van der Sleen et al.2015). The terrestrial biosphere module in SCP-M assumes a fixed δ13C fractionation factor of −23 ‰ (Table S1).

The major fluxes of carbon are parameterized simply in SCP-M to allow them to be solved by model–data optimization with respect to atmospheric and ocean proxy data. In this study the values for GOC, AMOC and biological export productivity at 100 m depth are outputs of the model–data experiments, as they are deduced from a data optimization routine. Their input values for the experiments are ranges, as described in Sect. 2.2.1. SCP-M's fast run time and flexibility renders it useful for long-term paleo-reconstructions involving large numbers of quantitative experiments and data integration (O'Neill et al.2019). SCP-M is a simple box model, which incorporates large regions of the ocean as averaged boxes and parameterized fluxes. It is an appropriate tool for this study, in which we evaluate many tens of thousands of simulations to explore possible parameter combinations, in conjunction with proxy data.

2.2 Model–data experiment design

We undertake series of model–data experiments to solve for the values of ocean circulation and biological parameters for each MIS during the last glacial–interglacial cycle (130–0 ka). We target these parameters due to their central role in many LGM–Holocene CO2 hypotheses (e.g. Knox and McElroy1984; Siegenthaler and Wenk1984; Toggweiler and Sarmiento1985; Martin1990; Kohfeld and Ridgwell2009; Hain et al.2010; Sigman et al.2010; Yu et al.2014a; Menviel et al.2016; Kohfeld and Chase2017; Muglia et al.2018; Menviel et al.2020). We force SST, salinity, sea volume and ice cover, and reef carbonate production, in each MIS (Sect. 2.2.1, Fig. 2), using values sourced from the literature (e.g. Opdyke and Walker1992; Key2001; Adkins et al.2002; Ridgwell et al.2003; Kohfeld and Ridgwell2009; Rohling et al.2009; Wolff et al.2010; Muscheler et al.2014; Kohfeld and Chase2017). Then, we optimize the model parameters for GOC, AMOC and Southern Ocean biological export productivity in each MIS time slice. We choose GOC and AMOC due to the prevalence of varying ocean circulation in many theories for glacial–interglacial cycles of CO2 (e.g. Sarmiento and Toggweiler1984; Siegenthaler and Wenk1984; Toggweiler1999; Kohfeld and Ridgwell2009; Burke and Robinson2012; Freeman et al.2016; Menviel et al.2016; Kohfeld and Chase2017; Skinner et al.2017; Muglia et al.2018; Menviel et al.2020) and its key role in distribution of carbon and other elements in the ocean (Talley2013). We choose to vary Southern Ocean biological export productivity due to its long-standing place and debate among theories of atmospheric CO2 during the LGM and Holocene (e.g. Martin1990; Knox and McElroy1984; Sarmiento and Toggweiler1984; Sigman and Boyle2000; Anderson et al.2002; Kohfeld and Ridgwell2009; Martinez-Garcia et al.2014; Menviel et al.2016; Kohfeld and Chase2017; Muglia et al.2018).

The GOC (Ψ1), AMOC (Ψ2) and Southern Ocean biology (Z) parameters are varied over  9000 possible combinations for each MIS, a total of  80 000 simulations across MIS 5e-1. At the end of each experiment batch, the model results are solved for the best fit to the ocean and atmosphere proxy data using a least-squares optimization, and the parameter values for Ψ1, Ψ2 and Z are returned. Our experiment time slices are the MIS of Lisiecki and Raymo (2005), with two minor modifications (see Fig. 2). MIS 2 (14–29 ka) as per Lisiecki and Raymo (2005) straddles the LGM (18–24 ka) and the last glacial termination (15–18 ka), while MIS 1 (0–14 ka) incorporates the Holocene period (0–11.7 ka) and the end of the termination. We are interested in the LGM and Holocene as discrete periods, so our experiment time slice for MIS 2 is truncated at 18 ka and our MIS 1 simply covers the Holocene, removing overlaps with the glacial termination. Therefore, our modelling excludes the last glacial termination ( 11–18 ka). The glacial termination period was highly transient with atmospheric CO2 varying by  85 ppm in < 10 kyr and large changes in carbon isotopes. Thus it is anticipated that in a model–data reconstruction, model parameters would vary substantially for this period. Joos et al. (2004), Ganopolski et al. (2010), Menviel et al. (2012), Menviel and Joos (2012), Brovkin et al. (2012), and Ganopolski and Brovkin (2017) provide coverage of the termination period with transient simulations, using intermediate-complexity models (more complex than our model). For MIS 5, we take the timing for peak glacial and interglacial substages of Lisiecki and Raymo (2005): ± 5 kyr for MIS 5c–5e and ± 2.5 kyr for MIS 5a–5b.

Figure 2Model forcings for MIS across the last glacial–interglacial cycle. (a) Sea surface temperature reconstruction of Kohfeld and Chase (2017), mean values mapped into SCP-M surface boxes (fine lines) and averaged across MIS (bold horizontal lines). (b) Proxy for Antarctic sea-ice extent using ssNa fluxes from the EPICA Dome C ice core (Wolff et al.2010), used to temporally contour MIS model forcings for (c) salinity (Adkins et al.2002) and (d) polar Southern Ocean air–sea gas exchange. Global ocean salinity is forced to a glacial maximum of +1 psu (shown in a) and the polar Southern Ocean is forced to +2 psu (not shown), as modified from Adkins et al. (2002). Ocean volume (e) forced using global relative sea level reconstruction of Rohling et al. (2009). (f) Atmospheric 14C production rate time series for 0–50 ka of Muscheler et al. (2014). Long-term values assumed for > 50 ka (Key2001). (g) Shallow-water carbonate flux of carbon from Ridgwell et al. (2003) profiled across the glacial–interglacial cycle using a curve from Opdyke and Walker (1992). Fine lines are the time series data and bold lines are the model forcings in each MIS. Data behind the figure are shown in Tables S2 and S3.


2.2.1 Model forcings and parameter variations

We take a reconstructed SST time series for the last 130 kyr (Kohfeld and Chase2017), map these to SCP-M's surface boxes and average the time series across each MIS (Fig. 2a). We extrapolate an Antarctic sea-ice cover proxy as shown in Fig. 2b (Wolff et al.2010) to the profiles for sea surface salinity (Fig. 2c) and the polar Southern Ocean box air–sea gas exchange parameter (Fig. 2d). For example, our notional reduction in the strength of the polar Southern Ocean box air–sea gas exchange due to Antarctic sea-ice cover (-30 %) is linearly (negatively) profiled with the Antarctic sea-ice proxy time series of Wolff et al. (2010). Note the polar Southern Ocean box, which is forced with reduced air–sea exchange, is separate from the subpolar Southern Ocean box in which the biological export productivity parameter is varied in the model–data experiment. Our treatment of sea-ice cover is simply as a regulator of air–sea gas exchange in the polar Southern Ocean surface boxes in each basin, not as a driver of other physical processes or biogeochemical feedbacks (e.g. Morrison and Hogg2013; Ferrari et al.2014; Jansen2017; Kohfeld and Chase2017; Marzocchi and Jansen2017). Furthermore, our linear application of the sea-ice proxy data of Wolff et al. (2010) to our air–sea gas exchange parameter (Fig. 2d) may overestimate its effect on the model results early in the glacial period (MIS 5d) and underestimate its effects during MIS 4–2 (Wolff et al.2010).

Adkins et al. (2002) reconstructed LGM deep-sea salinity for the Southern, Atlantic and Pacific oceans. They found increased salinity for the LGM at all locations across a range of +0.95-2.4 practical salinity units (psu) above modern values, with an average value of +1.5 psu. The most saline LGM waters were in the Southern Ocean (+2.4 psu), with Atlantic and Pacific waters ranging from +0.95-1.46 psu and a global ocean average of +1.2 psu. Adkins et al. (2002) also observed that within a (globally) more saline ocean, lower glacial temperatures would have caused less evaporation during the LGM, a negative feedback on salinity. We choose a forcing for LGM sea surface salinity of +1 psu for the global ocean and +2 psu for the polar Southern Ocean, relative to the interglacial period. These values conservatively reflect the hypothesis that surface evaporation may have been less in the LGM, hence resulting in a lesser magnitude of change in salinity in the surface ocean relative to the deep-ocean values estimated by Adkins et al. (2002), and also that the most voluminous parts of the ocean were less saline than the Southern Ocean (Adkins et al.2002). In our model–data experiments, the estimated glacial change in sea surface salinity (Fig. 2c) is also contoured through time with the variation in Antarctic sea-ice cover of Wolff et al. (2010). Adkins et al. (2002) observed that glacial salinity is a poor predictor of global mean sea level, due to storage of saline waters in ice shelves and groundwater reserves. Therefore, the proxy for Antarctic sea-ice cover may have a more direct linkage to sea surface salinity than using global sea level, for our purposes of estimating glacial–interglacial evolution in salinity.

Rohling et al. (2009) reconstructed global relative sea level (RSL) over the past five glacial–interglacial cycles. According to Rohling et al. (2009), the glacial RSL minimum was −115 m at  27 ka, immediately prior to the LGM. We perform a simple calculation to reduce ocean depth and volume in SCP-M, in line with the Rohling et al. (2009) time series. In a box model this is only an approximation, given the lack of topographical detail. Varying ocean box volume and surface area affects the ocean surface area available for in-gassing and degassing and the overall ocean capacity to store CO2, which impacts atmospheric CO2, δ13C and Δ14C (Köhler et al.2010; O'Neill et al.2019). Opdyke and Walker (1992) reconstructed coral reef carbonate fluxes of CaCO3 for the last glacial–interglacial cycle for the purposes of modelling the “coral reef hypothesis”. According to Opdyke and Walker (1992), reef carbon fluxes (out of the ocean) declined through the glacial cycle, with net dissolution in MIS 3 and MIS 2 leading to positive fluxes of carbon and alkalinity into the ocean in those periods. Fluxes of carbon and alkalinity out of the ocean into coral reefs rebounded from the LGM (MIS 2) into the Holocene (MIS 1), driven by increased sea level and temperature (Kleypas1997). Given that Opdyke and Walker (1992) evaluated the possibility for coral reefs to drive the entire glacial–interglacial CO2 variation, we take the more conservative modelling assumption of Ridgwell et al. (2003) of 0.5×1017 mol C for the postglacial accumulation of coral reefs. We profile this value across the glacial–interglacial cycle accumulation or dissolution curve of Opdyke and Walker (1992) as shown in Fig. 2. We apply the estimated atmospheric production rate for 14C for the last 50 kyr of Muscheler et al. (2014), with a long-term average production rate of  1.7 atoms cm−2 s−1 assumed for 130–50 ka (Key2001). Model forcing values are shown in Tables S2 and S3.

The terrestrial biosphere module in SCP-M does not explicitly represent the carbon stored in buried peat, permafrost and also cold-climate vegetation that may have expanded its footprint in the glaciation, such as tundra biomes (e.g. Tarnocai et al.2009; Ciais et al.2012; Schneider et al.2013; Eggleston et al.2016; Ganopolski and Brovkin2017; Treat et al.2019). The freezing and burial of organic matter across the glacial period sequesters carbon on land and may modify atmospheric CO2 and δ13C (Tarnocai et al.2009; Ciais et al.2012; Schneider et al.2013; Eggleston et al.2016; Ganopolski and Brovkin2017; Mauritz et al.2018; Treat et al.2019). Ganopolski and Brovkin (2017) incorporated permafrost, peat, and buried land carbon into their transient simulations of the last four glacial–interglacial cycles with the CLIMBER-2 model. Ganopolski and Brovkin (2017) observed that these features dampened the amplitude of glacial–interglacial variations in terrestrial biosphere carbon stock and its effects on atmospheric CO2. As a crude measure to account for this counter-CO2 cycle storage of carbon in the terrestrial biosphere and frozen soils or buried carbon, we force the terrestrial biosphere productivity parameter in SCP-M in the range +5-10 Pg C yr−1 throughout the last glacial–interglacial cycle, increasing into the LGM (MIS 2) and maintained in the Holocene (MIS 1). We maintain this forcing in the Holocene, as the posited effects of buried peat and permafrost storage of carbon on atmospheric CO2 and δ13C during the lead-up to the LGM were likely not reversed after the glacial termination (Tarnocai et al.2009; Eggleston et al.2016; Mauritz et al.2018; Lindgren et al.2018; Treat et al.2019). SCP-M calculates NPP using this productivity input parameter and a logarithmic function of carbon fertilization (Harman et al.2011).

More than 9000 model simulations are undertaken across the parameter ranges in Table 1 for each MIS. Parameters are varied simultaneously to allow coverage of all possible combinations of the parameter values within their respective experiment ranges. Within these ranges, values are incremented by 1 Sv for GOC (Ψ1) and AMOC (Ψ2), and  0.5 mol C m−2 yr−1 for Atlantic Southern Ocean biological export productivity (Z). Each simulation is run for 10 kyr to enable the model to achieve steady state. We show the experiment ranges for the biological export productivity parameter Z for the Atlantic and Pacific–Indian sectors of the Southern Ocean (Table 1). In SCP-M, the Pacific–Indian Southern Ocean biological export productivity parameter (in mol C m−2 yr−1) is set by default at a value of  70 % of the corresponding Atlantic sector Southern Ocean box, to align with natural observations of variations in the Southern Ocean biological export productivity (e.g. Dunne et al.2005; Sarmiento and Gruber2006; Henson et al.2011; Siegel et al.2014; DeVries and Weber2017). This variation is reflected in the values in Table 1. In the experiments, the values for Z in the Pacific–Indian Southern Ocean surface box scale linearly with the values for the Atlantic Southern Ocean surface box (Table 1). Herein we focus our presentation and discussion of the experiment results for the Z parameter on the Atlantic Southern Ocean due to its prominence in glacial–interglacial cycle hypotheses for increased biological productivity (e.g. Martinez-Garcia et al.2014; Lambert et al.2015; Shaffer and Lambert2018; Muglia et al.2018).

Table 1Free-floating parameter ranges in the model–data experiments for global overturning circulation (GOC, Ψ1), Atlantic meridional overturning circulation (AMOC, Ψ2) and Southern Ocean biological export productivity (Z). Parameters are varied simultaneously across these ranges and then optimized against proxy data in each MIS. Also shown are pre-industrial control values for GOC (Talley2013), AMOC (Talley2013) and Southern Ocean biological export productivity (Dunne et al.2005; Sarmiento and Gruber2006; Henson et al.2011; Siegel et al.2014; DeVries and Weber2017). The Pacific–Indian Southern Ocean biology parameter is set at a base value of  70 % Atlantic Southern Ocean box but scales linearly with the Atlantic Ocean parameter in the experiments. The smaller values for Pacific–Indian Southern Ocean take account of natural observations of a relatively stronger biological export productivity in the Atlantic sector of the subpolar Southern Ocean (e.g. Dunne et al.2005; Sarmiento and Gruber2006; Henson et al.2011; Siegel et al.2014; DeVries and Weber2017).

Download Print Version | Download XLSX

2.2.2 Optimization procedure

We perform a least-squares optimization of the model experiment output against MIS data for atmospheric CO2, atmospheric, deep and abyssal ocean Δ14C and δ13C, and deep and abyssal ocean carbonate ion proxy, to source the best-fit parameter values for GOC, AMOC and Southern Ocean biological productivity in each time slice – a brute force form of the gradient descent method for optimization (e.g. Strutz2016). The equation for least fit applied is

(3) Opt n = Min i , k = 1 N ( R i , k - D i , k σ i , k ) 2 ,

where Optn is the optimal value of parameters n (e.g. GOC, AMOC and Southern Ocean biological productivity), Ri,k is model output for concentration of each element i in box k, Di,k is average data concentration each element i in box k and σi,k is standard deviation of the data for each element i in box k. The standard deviation performs two roles. It normalizes for different unit scales (e.g. ppm, ‰ and µmol kg−1), which allows multiple proxies to be incorporated in the optimization, and reduces the weighting of a proxy data point with a high standard deviation and therefore an uncertain value. The weighting by proxy data standard deviation also fulfils the important role of accounting for data variance in the optimized parameter results, such that the effects of data variance are embedded in the optimized parameter values. Where proxy data are unavailable for a box, that data and box combination is automatically omitted from the optimization routine. The experiment routine returns the model run with the best fit to the data and the model's parameters and results.

Monnin et al. (2004)MacFarling Meure et al. (2006)Bereiter et al. (2012)Rubino et al. (2013)Schneider et al. (2013)Ahn and Brook (2014)Marcott et al. (2014)Bereiter et al. (2015)Elsig et al. (2009)Schmitt et al. (2012)Schneider et al. (2013)Eggleston et al. (2016)Reimer et al. (2009)Oliver et al. (2010)Govin et al. (2009)Piotrowski et al. (2009)Skinner and Shackleton (2004)Marchitto et al. (2007)Barker et al. (2010)Bryan et al. (2010)Skinner et al. (2010)Burke and Robinson (2012)Siani et al. (2013)Davies-Walczak et al. (2014)Skinner et al. (2015)Chen et al. (2015)Hines et al. (2015)Sikes et al. (2016)Ronge et al. (2016)Skinner et al. (2017)Zhao et al. (2017)Yu et al. (2010)Yu et al. (2013)Yu et al. (2014b)Yu et al. (2014a)Broecker et al. (2015)Yu et al. (2016)Qin et al. (2017)Qin et al. (2018)Chalk et al. (2019)

Table 2Ocean and atmosphere proxy data sources for the last glacial–interglacial cycle

Download Print Version | Download XLSX

2.3 Data

Our model–data optimization rests on compilations of atmospheric and ocean paleo proxy data. We compile and apply published proxy data for atmospheric CO2, δ13C and Δ14C and ocean δ13C, Δ14C and CO32- concentration. We calculate the simple mean and standard deviation of data points for each model box and MIS. The proxy data for each ocean box is binned into a model box based on depth, latitude and longitude, which assigns the data to either the Atlantic or Pacific–Indian basin. The box-mapped data are binned into MIS age groups and the sample population is then averaged and the standard deviation is calculated. The standard deviation is then used as a weighting in the model–data optimization procedure. Sources of proxy data are shown in Table 2 and data locations in Fig. 3. MIS and model box-averaged atmospheric and ocean proxy data and their respective standard deviations are shown in Tables S4–S7.

2.3.1 Ocean carbon isotopes

We gather published marine Δ14C data extending back to  40 ka (Table 2). Our dataset incorporates individual records contributed over the last 30 years and supplemented by the recent compilations of Skinner et al. (2017) and Zhao et al. (2017). The data total  75 individual location estimates for benthic and planktonic foraminifera and deep-sea corals. We have restricted our efforts to time series which contain independent calendar ages and therefore corrections for radioactive decay in the time since the sample was deposited (yielding Δ14C). Figure 3 shows the geographic distribution of the Δ14C data, which is generally concentrated on ocean basin margins. Some regions, such as the central Pacific, southern Indian and polar Southern Ocean, are devoid of data.

Figure 3Δ14C, δ13C and CO32- data locations. Δ14C and CO32- data are compiled from published estimates. For δ13C we take the compilation of Oliver et al. (2010). MIS and model box-averaged data and their respective standard deviations are shown in Tables S4–S7.

Oliver et al. (2010) compiled a global dataset of 240 cores of marine δ13C data encompassing benthic and planktonic species for the last  150 kyr. Oliver et al. (2010) observed considerable uncertainties associated with the broad range of species included, particularly for the planktonic foraminifera. By comparison, Peterson et al. (2014) aggregated marine δ13C for the LGM and late Holocene periods, as time period averages, exclusively sampling benthic C. wuellerstorfi data, which are a more reliable indicator of marine δ13C (Oliver et al.2010; Peterson et al.2014). To narrow the range of uncertainty, we constrain our use of marine δ13C data to the deep and abyssal (> 2500 m) benthic Cibicides species foraminifera samples in the Oliver et al. (2010) dataset, supplemented with Cibicides species δ13C proxy data from Govin et al. (2009) and Piotrowski et al. (2009) (Table 2). Figure 3 shows the δ13C data locations from Oliver et al. (2010), which are concentrated in the Atlantic Ocean. We map and average the carbon isotope data into SCP-M's boxes on depth and latitude coordinates (Fig. 1), averaged for each MIS time slice.

2.3.2 Carbonate ion proxy

We aggregate ocean carbonate ion proxy data (as deduced from B∕Ca) from the sources shown in Table 2 and locations in Fig. 3, map into SCP-M box coordinates, and average the data across MIS. The data coverage for CO32- is relatively sparse, with < 20 individual site locations across the global ocean. However, the depth and lateral coverage of SCP-M's boxes is large, particularly in the case of the deep-ocean boxes, which cover the full lateral extent of the Pacific–Indian and Atlantic oceans, and depth ranges of 100–2500 m (Pacific–Indian) and 250–2500 m (Atlantic). CO32- can vary by more than 100 µmol kg−1 across the depth range 100–2500 m and can vary by up to  200 µmol kg−1 in the shallow ocean (e.g. Sarmiento and Gruber2006; Yu et al.2014b, a). Some boxes contain only one core, creating an exceptionally low standard deviation range relative to the other ocean proxies. In other cases, such as the deep Atlantic Ocean, the data points are clustered within the 2000–2500 m depth range, the bottom third of the corresponding SCP-M box. This clustering becomes a problem for the SCP-M box model, which outputs average concentrations over the complete depth range of each box – a drawback of using a large-resolution box model to analyse proxy data at a global ocean level. Furthermore, the very low standard deviations associated with the CO32- data (shown in Table S6) cause it to assume a disproportionate weighting in the model–data optimization, which uses standard deviation for weighting of proxies, relative to ocean δ13C and Δ14C. The latter proxies often have box standard deviations up to 100 % of their mean value, when averaged across a box. This issue is also an artefact of our procedure necessary to normalize the different proxies (each in unique units) in a multi-proxy model–data optimization, by using the standard deviation as a weighting. To deal with this, we assign an arbitrary standard deviation (weighting) of 20 µmol kg−1 to CO32- data observations in our model–data optimizations, which acts as a feasible weighting for the processing of CO32- relative to the other ocean proxy data. This value is a small fraction of the variation in CO32- concentrations observed over the depth range 100–2500 m in the modern ocean (e.g. Key et al.2004; Yu et al.2014b).

3 Data analysis

In this section we describe the proxy data used to constrain the glacial–interglacial model–data experiments. We depict the major changes in atmospheric CO2, δ13C and Δ14C and ocean δ13C, Δ14C and CO32- proxy data across the model box locations and MIS in the last glacial–interglacial cycle. We mainly refer to changes in the MIS-averaged proxy data.

Figure 4MIS atmosphere data for (a) atmospheric CO2 (Monnin et al.2004; MacFarling Meure et al.2006; Bereiter et al.2012; Rubino et al.2013; Schneider et al.2013; Ahn and Brook2014; Marcott et al.2014), (b) δ13C (Elsig et al.2009; Schmitt et al.2012; Schneider et al.2013; Eggleston et al.2016) and (c) Δ14C (Reimer et al.2009). Data are shown as fine lines, with bold horizontal lines for MIS-sliced data. Natural observations for Δ14C do not exist beyond  50 ka due to the radioactive decay of 14C. Data behind the figure are shown in Table S4.


Figure 4 shows the atmospheric data used to constrain the model–data experiments, averaged into MIS time slices. There are many fluctuations and transient changes throughout the last glacial–interglacial cycle, but there are three major sustained reductions in atmospheric CO2 concentration in the lead-up to the LGM (Fig. 4a): an average drop of  25 ppm during MIS 5d (115–100 ka), a further average drop of  30 ppm during MIS 4 (72–65 ka) and finally a fall of  20 ppm in the period leading up to the LGM (during MIS 3 and 2, 40–18 ka). These are the three major CO2 events described in Kohfeld and Chase (2017) (although MIS-averaged in our analysis) and, combined with additional reductions of −10 ppm throughout the period, yield a total drop of −85 ppm from the last interglacial to the LGM. Transient changes in atmospheric CO2 concentration occur throughout the glacial cycle, including during MIS 5c–5a, MIS 4 and throughout MIS 3. As discussed in the Introduction, this sequence of CO2 reductions is likely the result of oceanic drivers with biogeochemical and terrestrial feedbacks (e.g. Ganopolski et al.2010; Menviel et al.2012; Brovkin et al.2012; Ganopolski and Brovkin2017; Kohfeld and Chase2017). Atmospheric CO2 concentration increases by  85 ppm in the glacial termination and Holocene periods, a transition in the carbon cycle which has occupied substantial research effort in the last 4 decades but with a growing consensus of multiple physical and biogeochemical drivers and feedbacks. Kohfeld and Ridgwell (2009) and Köhler et al. (2010) provide summaries of the potential candidate mechanisms to explain the glacial–interglacial changes in atmospheric CO2, while recent model–data studies have attempted to explain the specific physical and biogeochemical drivers of the LGM–Holocene change in atmospheric CO2 (Tagliabue et al.2009; Menviel et al.2016; Muglia et al.2018; O'Neill et al.2019).

Figure 4b shows atmospheric δ13C over the last glacial–interglacial cycle. Eggleston et al. (2016) explained the glacial–interglacial atmospheric δ13C pattern in terms of ongoing changes in SST, AMOC, Southern Ocean upwelling, dust-driven Southern Ocean biological export productivity and the terrestrial biosphere. Atmospheric δ13C (Fig. 4b) was  0.4‰ higher in the Holocene (MIS 1) and LGM (MIS 2) periods than in the last interglacial (MIS 5e) and penultimate glacial periods (MIS 6, not shown in Fig. 4b), as described in Schneider et al. (2013) and Eggleston et al. (2016). There were temporary falls in atmospheric δ13C between MIS 5e and 5d (between 120 and 110 ka), during MIS 4 (between 69 and 58 ka), during MIS 3 (between 50 and 35 ka), and in the last glacial termination between MIS 2 and 1 (between 19 and 16 ka). The cause of the observed increase in atmospheric δ13C across the last glacial–interglacial cycle may be the effect of accumulation and freezing or burial in glacial sediments of peat and other soil organic matter at the high latitudes (e.g. Tarnocai et al.2009; Ciais et al.2012; Schneider et al.2013; Eggleston et al.2016; Ganopolski and Brovkin2017; Treat et al.2019). According to Treat et al. (2019), peatlands and other vegetation accumulated carbon in the relatively warm periods, and these carbon stocks were then frozen and/or buried in glacial and other sediments during the cooler periods, throughout the last glacial–interglacial cycle. This buried or frozen stock of carbon mostly persists to the present day (Tarnocai et al.2009; Ciais et al.2012). Schneider et al. (2013) evaluated several possible candidates for the rising atmospheric δ13C pattern across the last glacial–interglacial cycle and could not discount any (1) changes in the carbon isotope fluxes of carbonate weathering and sedimentation on the seafloor, (2) variations in volcanic outgassing, or (3) peat and permafrost build-up throughout the last glacial–interglacial cycle.

The large drop in atmospheric δ13C observed during MIS 4 reverses in MIS 3 (Fig. 4b). This excursion in the δ13C pattern likely resulted from sequential changes in SST (cooling), AMOC, Southern Ocean upwelling and marine biological productivity (Eggleston et al.2016). Eggleston et al. (2016) parsed the atmospheric δ13C signal into its component drivers across MIS 5a–3 using a stack of proxy indicators. Eggleston et al. (2016) highlighted the sequence of events between the end of MIS 5a and beginning of MIS 3 and their cumulative effects to deliver the full change in atmospheric δ13C. Our MIS-averaging approach as shown in Fig. 4b fails to capture the full amplitude of the changes in atmospheric δ13C during MIS 4 and MIS 3 and only captures the changes in the mean-MIS value, serving to understate the full extent of transient changes in responsible processes. In addition, the MIS-averaging approach misses the sequential timing of changes in processes within each MIS. These are limitations of our steady-state, MIS-averaging approach. The reduction in atmospheric δ13C at the last glacial termination, between the LGM and Holocene (Fig. 4b), coincident with a large atmospheric CO2 increase, is attributed to the release of deep-ocean carbon to the atmosphere resulting from increased ocean circulation and Southern Ocean upwelling (Schmitt et al.2012). The subsequent rebound of δ13C in the termination period and the Holocene is believed to result from terrestrial biosphere regrowth, in response to increased CO2 and carbon fertilization (Schmitt et al.2012; Hoogakker et al.2016).

Figure 4c shows atmospheric Δ14C over the last 50 kyr (Reimer et al.2009). During this period Δ14C is heavily influenced by declining atmospheric 14C production (Broecker and Barker2007; Muscheler et al.2014). In addition, an acceleration in atmospheric Δ14C decline at the last glacial termination is attributed to the release of old, 14C-depleted waters from the deep ocean, due mainly to increased Southern Ocean upwelling of Δ14C-depleted deep source waters (e.g. Marchitto et al.2007; Skinner et al.2010; Burke and Robinson2012; Siani et al.2013). Broecker and Barker (2007) characterized the drop in atmospheric Δ14C at the last glacial termination as “the mystery interval” and questioned whether there existed a Δ14C-depleted ocean reservoir source of sufficient size to contribute to the drop.

Figure 5 shows deep and abyssal ocean δ13C data mapped into SCP-M box model space and averaged across MIS. The visual offset between deep and abyssal proxy data values is regularly interpreted as an indicator of the strength of deep-ocean circulation and/or mixing, or biological productivity, during the LGM and the Holocene (e.g. Sikes et al.2000; Curry and Oppo2005; Marchitto et al.2007; Oliver et al.2010; Skinner et al.2010; Burke and Robinson2012; Siani et al.2013; Yu et al.2013, 2014a; Skinner et al.2015, 2017). The deep–abyssal Atlantic δ13C time series (Fig. 5a) exhibits modest widening in the MIS-average deep and abyssal offset between MIS 5e and 5d and again during MIS 5b and then a more substantial widening during MIS 4 and during MIS 2 (the LGM). The widening of the offset during MIS 4 and MIS 2 is caused primarily by more negative abyssal δ13C values. The offset is almost closed in MIS 1 (the Holocene). The deep Atlantic δ13C range itself also widens considerably from MIS 4 and narrows after the LGM. Oliver et al. (2010) and Kohfeld and Chase (2017) interpreted these patterns as the result of weakened deep Atlantic Ocean circulation during MIS 4 and during the LGM, strengthening in the post glacial period.

The Pacific–Indian δ13C data (Fig. 5b) show a drop in abyssal δ13C and widening in the MIS-average deep–abyssal offset between MIS 5e and 5d (Govin et al.2009) which continued throughout the last glacial build-up. Importantly, the more negative abyssal δ13C values during MIS 5d–5a seen in Fig. 5b occur at the same time that deep-ocean and atmospheric δ13C becomes more positive (Fig. 4b), suggesting that the abyssal Pacific–Indian oceans became more isolated from the deep ocean and atmosphere during this period. This is qualitative evidence for slowing ocean circulation or increased biological export productivity in the Pacific–Indian oceans, at that time (Govin et al.2009). This also corresponds with a  50 ppm fall in CO2 across the period spanning MIS 5e to 5b (Fig. 4a). Abyssal Pacific–Indian δ13C drops further and most noticeably during MIS 4 and again during the LGM and then rebounds from the LGM into the Holocene period, as also observed in the Atlantic Ocean δ13C data. Statistical analysis of the δ13C data provided in Fig. S1 and Table S8 supports our qualitative interpretation of the Atlantic and Pacific–Indian δ13C proxy data.

Figure 5MIS ocean data mapped into SCP-M box model dimensions for δ13C (Govin et al.2009; Piotrowski et al.2009; Oliver et al.2010). Data (round circles) are mapped into deep (2500 m average depth) and abyssal (3700 (Atlantic) – 4000 m (Pacific–Indian) average depth) model boxes and averaged across MIS slices (bold lines). Data behind the figure are shown in Table S5.


Ocean Δ14C data cover the period MIS 1–3 and the LGM and Holocene in most detail (Fig. 6). We show ocean ΔΔ14C, which is ocean minus atmospheric Δ14C. This calculation is made in attempt to normalize the effects of varying atmospheric 14C production through the glacial–interglacial cycle (Broecker and Barker2007; Muscheler et al.2014), which imparts a dominant influence on the ocean Δ14C trajectory. Given the sparse data coverage for MIS 3 we focus our analysis on MIS 1 and 2. The ΔΔ14C time series exhibits two key features across the MIS 2 (LGM) and MIS 1 (Holocene) periods. First, there is a narrowing in the spread of values between the shallow and abyssal ocean from the LGM to the Holocene, in both the Atlantic (Fig. 6a) and Pacific–Indian (Fig. 6b) basins. Second, all ocean boxes display an increase in ΔΔ14C from the LGM to the Holocene, towards equilibrium with the atmosphere. These patterns are believed to represent increased overturning circulation and Southern Ocean upwelling in the Atlantic and Pacific–Indian basins across the LGM–Holocene. Increased ocean overturning brought old, Δ14C-negative water up from the deep and abyssal oceans, resulting in mixing with shallow and intermediate waters and eventually into the surface Southern Ocean and contact with the atmosphere (where 14C is produced) – known as “increased ventilation” (e.g. Sikes et al.2000; Marchitto et al.2007; Bryan et al.2010; Skinner et al.2010; Burke and Robinson2012; Siani et al.2013; Davies-Walczak et al.2014; Skinner et al.2014; Hines et al.2015; Freeman et al.2016; Sikes et al.2016; Skinner et al.2017).

Figure 6MIS stage ocean data mapped into box model dimensions for ΔΔ14C. Data (round circles) are mapped into deep (2500 m average depth) and abyssal (3700 (Atlantic)–4000 m (Pacific–Indian) average depth) model boxes and averaged across MIS slices (bold lines). Natural observations do not exist beyond  50 ka due to the radioactive decay of 14C. ΔΔ14C is ocean minus atmosphere Δ14C. Note that this calculation is not done with the average ocean box and atmosphere values for each MIS; rather ΔΔ14C represents the difference between each ocean data point and the contemporary atmospheric Δ14C value. Data behind the figure are shown in Table S7.


The Atlantic Ocean CO32- time series shows a similar pattern to ΔΔ14C and δ13C, with a wide dispersion of shallow–abyssal and deep–abyssal concentrations at the LGM that narrows by the Holocene (Fig. 7). This pattern has been interpreted as varying strength and/or depth of AMOC and biological productivity in the Atlantic Ocean (e.g. Yu et al.2013, 2014b, a, 2016). The abyssal Atlantic CO32- pattern, which spans the last glacial–interglacial cycle, is punctuated by two downward excursions (Fig. 7). These occur during MIS 4 and MIS 2, corresponding to the second and third major atmospheric CO2 drops in the last glacial–interglacial cycle (Kohfeld and Chase2017), respectively (Fig. 4a). The lower deep Atlantic Ocean CO32- values during MIS 4 were interpreted by Yu et al. (2016) as shoaling of AMOC and increased carbon storage in the deep–abyssal Atlantic Ocean. This signal is repeated at the LGM, where further shoaling and slowing AMOC contributed to deep oceanic drawdown of CO2 from the atmosphere (Yu et al.2013, 2014b, a). There is also a modest drop in abyssal Atlantic Ocean CO32- during MIS 5b (−13µmol kg−1 relative to MIS 5c), which coincides with a minor drop in abyssal Atlantic Ocean δ13C (−0.19 ‰) and atmospheric CO2 (−14 ppm), indicating a common link. Menviel et al. (2012) modelled a transient slowdown in AMOC for this period, which could explain these features.

The Pacific Ocean is thought to partially buffer the effects of ocean circulation on CO32- concentrations (Fig. 7) via changes in shallow (reef) and deep carbonate production and dissolution and therefore displays less variation across the MIS (Yu et al.2014b; Qin et al.2017, 2018). The deep and abyssal Pacific–Indian ocean data show a gradual trend of increasing CO32- through the glacial–interglacial cycle (Fig. 7), suggesting that it is influenced more by variations in shallow- or deep-sea carbonate production or dissolution and less by deep-ocean circulation (Yu et al.2014b; Qin et al.2017, 2018). Notable exceptions are during MIS 5d and MIS 4. Between MIS 5e and 5d, both deep and abyssal Pacific–Indian ocean CO32- drops (Fig. 7), aligning with the contemporary drop in abyssal ocean δ13C and atmospheric CO2 (Figs. 5b and 4a), suggesting a possible common driver, and providing additional qualitative evidence for changes in either Pacific–Indian ocean circulation or biology, at this time. During MIS 4, there is a drop in deep and abyssal Pacific–Indian CO32- and a modest widening in the average deep–abyssal offset from MIS 5b and 5a, also suggestive of the influence of deep-ocean circulation and/or biological export productivity (Fig. 7). The widest Pacific–Indian deep–abyssal offset CO32- is observed during MIS 3, also seen in the ΔΔ14C data (Figs. 57), indicating it is a persistent feature of the proxy records. This suggests MIS 3 may be the nadir of Pacific–Indian ocean circulation and/or the peak in biological activity in the last glacial–interglacial cycle or at least that important changes in this part of the ocean took place in MIS 3, prior to the LGM.

Figure 7MIS stage ocean data mapped into box model dimensions for carbonate ion proxy. Data (round circles) are mapped into deep (2500 m average depth) and abyssal (3700 (Atlantic)–4000 m (Pacific–Indian) average depth) model boxes and averaged across MIS slices (bold lines). Data behind the figure are shown in Table S6.


4 Results

Figure 8 shows the data-optimized MIS-average values returned from the model–data experiments for GOC, AMOC and Atlantic Southern Ocean biological productivity parameters, in each MIS (“X” symbols). The optimized values take account of data variance, due to the weighting of proxy data points by their standard deviation in the model–data optimization equation (Eq. 3). The full range of model–data experiment results are shown in Table S9. The GOC parameter (Ψ1) value falls from 29 to 22 Sv between MIS 5e and 5d, with gradual declines during MIS 5c–5a and a slight acceleration in the rate of decline during MIS 5a–3. GOC reaches its minimum glacial value (16 Sv) in MIS 3, then increases from 16 to 29 Sv between MIS 2 (the LGM) and the Holocene. AMOC (Ψ2) weakens modestly in MIS 5d (−2 Sv), with a further drop during MIS 5b (−2 Sv) that is partially reversed in MIS 5a. AMOC weakens further in MIS 4, achieving its glacial nadir (13 Sv), which is maintained until the LGM before increasing to 18 Sv in MIS 1. Importantly, Ψ2 closely follows the abyssal Atlantic (>2500 m, single box covering North and South Atlantic) δ13C and CO32- data patterns across the glacial–interglacial cycle and ΔΔ14C from the LGM to the Holocene (Figs. 57). Ψ2 remains near its modelled last interglacial value (MIS 5e, 18 Sv), during MIS 5d and 5c before dropping in MIS 5b (abyssal Atlantic δ13C and CO32- and atmospheric CO2 also drop at this point) and partly rebounding during MIS 5a and then falling synchronously with abyssal Atlantic δ13C and CO32- concentrations during MIS 4. Southern Ocean biological export productivity (Z) fluctuates around its last interglacial (MIS 5e) value during the time period spanning MIS 5d–5b and then increases during MIS 4. Atlantic (Pacific–Indian) Southern Ocean Z spikes to 4.7 (3.3) mol C m−2 yr−1 in the LGM and then falls to 3.4 (2.4) mol C m−2 yr−1 in MIS 1.

Figure 8Model–data experiment results for global overturning circulation (a), Atlantic meridional overturning circulation (b) and Atlantic Southern Ocean biological export productivity (c). “X” symbols mark the optimal parameter values returned from the model–data experiments. The optimized values take account of data variance, due to the weighting of proxy data points by their standard deviation in the model–data optimization equation (Eq. 3). Data for optimized parameter values shown in the figure are contained in Table S9.


Figure 9 shows the optimized model–data output for atmospheric CO2 and ocean CO32- concentrations compared with the proxy data observations, in each MIS. This shows how well the model is constrained by the proxy data and also how well the model–data output of parameter values can explain the proxy data patterns as described in the data analysis section (Sect. 3). The model–data results fall within 1 standard deviation of atmospheric CO2 and deep and abyssal CO32- data and mostly on the MIS means, across the MIS periods (Fig. 9). The modelled abyssal Pacific–Indian CO32- falls close to the MIS proxy data means across the glacial–interglacial cycle but misses some of the variations in the data – particularly between MIS 4 and MIS 3 (Fig. 9). This is a result of the abyssal ocean box carbonate dissolution equations in SCP-M, which effectively buffer changes in the relative balance of DIC and alkalinity from ocean physical and biological changes, and possibly the large box sizes in SCP-M which miss some detail for sparse CO32- data.

Figure 9Values returned from the model–data experiment for (a) atmospheric CO2 and carbonate ion proxy for (b) deep Atlantic (2500 m average depth), (c) abyssal Atlantic (3700 m average depth), (d) deep Pacific–Indian (2500 m average depth) and (e) abyssal Pacific–Indian (4000 m average depth). Model–data experiment results are shown as dots, with mean proxy data shown as solid lines and 1 standard deviation range by dashed lines, in each MIS. A default standard deviation of 20 µmol kg−1 is used as discussed in the text. CO32- data for the SCP-M deep Atlantic box in (b) do not extend beyond 50 ka. Model results for each box in each MIS are shown in Tables S10 and S12.


The model–data results show good agreement with atmospheric, deep and abyssal δ13C data throughout the MIS (Fig. 10). The results mostly fall on the mean, and all are within the standard deviation for atmospheric δ13C data in the MIS. Nearly all results fall within standard deviation for the deep and abyssal Atlantic and Pacific–Indian oceans. The modelled abyssal Pacific–Indian box δ13C underestimates mean MIS δ13C in most MIS time slices, which may reflect a discrepancy between the average depth of the δ13C proxy data and SCP-M abyssal ocean box or a bias in the model's equations.

Figure 10Values returned from the model–data experiment for δ13C for (a) atmosphere, (b) deep Atlantic (2500 m average depth), (c) abyssal Atlantic (3700 m average depth), (d) deep Pacific–Indian (2500 m average depth) and (e) abyssal Pacific–Indian (4000 m average depth). Model–data experiment results are shown as dots, with proxy data mean (solid lines) and 1 standard deviation (dashed lines) in each MIS. Model results for each box in each MIS are shown in Tables S10 and S11.


Figure 11 shows model–data results for atmospheric Δ14C and ocean ΔΔ14C compared with data, for MIS 1–3. Model–data results fall within 1 standard deviation of the data for all observations that were modelled and replicate the dramatic compression in deep–abyssal ΔΔ14C and ocean–atmosphere offsets between MIS 2 (LGM) and MIS 1 (the Holocene) as shown in the data (Fig. 11).

Figure 11Values returned from the model–data experiment for (a) atmospheric Δ14C and ΔΔ14C for (b) deep Atlantic (2500 m average depth), (c) abyssal Atlantic (3700 m average depth), (d) deep Pacific–Indian (2500 m average depth) and (e) abyssal Pacific–Indian (4000 m average depth). ΔΔ14C is ocean minus atmospheric Δ14C, calculated to correct for the varying atmospheric Δ14C signal. Model–data experiment results are shown as dots, with proxy data mean (solid lines) and 1 standard deviation (dashed lines) in each MIS. Model–data experiment results prior to MIS 4 are omitted, due to the radioactive decay of 14C, which precludes natural observations prior to  50 ka. Model results for each box in each MIS are shown in Tables S10 and S13.


Figure 12 shows model–data output for the terrestrial biosphere NPP and carbon stock during the last glacial–interglacial cycle. The NPP and carbon stock follow atmospheric CO2 downwards in the lead-up to the LGM and rebound from the LGM to the Holocene. In our model this is driven by carbon fertilization from atmospheric CO2 (Kaplan et al.2002; Otto et al.2002; Harman et al.2011; Hoogakker et al.2016). However, other studies emphasize the important role of temperature and precipitation in influencing NPP (François et al.1999; van der Sleen et al.2015). Notably, there is a distinct drop in NPP during MIS 4, a period where atmospheric CO2 falls by  30 ppm (Fig. 4a). Hoogakker et al. (2016) provided a reconstruction of NPP through the last glacial–interglacial cycle using pollen data and climate models, shown for comparison in Fig. 12a. Our model–data results for NPP typically fall in the upper and lower end of the range of NPP values from Hoogakker et al. (2016). However, our model–data estimates of NPP for MIS 5d and 5e underestimate the NPP calculated by Hoogakker et al. (2016) (which extends only to 120 ka). We model the terrestrial biosphere carbon stock as falling by 385 Pg C from the last interglacial to the LGM and increasing by  600 Pg C from the LGM to the Holocene (Fig. 12b).

Figure 12(a) Model–data output for the terrestrial biosphere net primary productivity (NPP) in each MIS time slice (black dashed lines) compared with the range of estimates provided by Hoogakker et al. (2016) (grey area). (b) Model–data output for the terrestrial biosphere carbon stock for each MIS time slice.


5 Discussion

5.1 Last glacial–interglacial cycle

This study applies a carbon cycle box model to diagnose the values for ocean circulation and Southern Ocean biological export productivity during the last glacial–interglacial cycle, optimized for ocean and atmospheric proxy data. This study continues efforts to simulate the last glacial–interglacial cycle of atmospheric CO2 (e.g. Ganopolski et al.2010; Brovkin et al.2012; Menviel et al.2012; Ganopolski and Brovkin2017) but with a simpler box model and using a non-transient model–data optimization to estimate parameter values.

There were three major episodes in which atmospheric CO2 concentration fell during the last glacial–interglacial cycle (Fig. 4a), accompanied by changes in atmospheric δ13C (Fig. 4b), Δ14C (Fig. 4c) and ocean δ13C, Δ14C and CO32- (Figs. 57). Our model–data results show that glacial–interglacial atmospheric CO2 and the other proxy patterns are delivered by a host of physical and biogeochemical changes. These changes include weakened GOC, AMOC and strengthened Southern Ocean biological export productivity (Figs. 8, 9, 10, 11) and changes in SST, salinity, ocean volume, the terrestrial biosphere, reef carbonates and atmospheric 14C production (Figs. 2 and 12).

Our model–data results show that an initial fall in GOC took place during MIS 5d (Fig. 8), as MIS-average atmospheric CO2 concentration fell by  25 ppm. This was also a time of substantial cooling in SST (Fig. 2a). GOC drifted lower until achieving its glacial minimum level in MIS 3 and MIS 2. AMOC weakened in MIS 4, at the same time that North Atlantic SST cooled dramatically (Fig. 2a) and MIS-average atmospheric CO2 fell  30 ppm. GOC and AMOC were both equal to their glacial lows at the LGM and accompanied by increased Southern Ocean biological export productivity, yielding the LGM minima in atmospheric CO2 and the final major fall in CO2 during the glacial cycle. We model elevated Southern Ocean biological productivity during MIS 4 and MIS 2, relative to interglacial values (MIS 5e and 1). Importantly, the transition from MIS 3 to MIS 2, which incorporates the LGM and increased Southern Ocean biological productivity, only accounted for an average 15 ppm reduction in CO2 (Figs. 4, 9). Therefore, our results suggest that an increase in Southern Ocean biological productivity during this period was an additional “kicker” to achieve the LGM atmospheric CO2 minima, following prior reductions of  70 ppm in the lead-up, which were delivered mainly by ocean physical processes and SST. The finding of increased biological productivity, while mostly constrained in our model to MIS 4 and 2, and a modest yet essential contributor to the overall glacial CO2 drawdown corroborates proxy data (e.g. Martinez-Garcia et al.2014; Lambert et al.2015; Kohfeld and Chase2017; Shaffer and Lambert2018) and recent model–data exercises (e.g. Menviel et al.2016; Muglia et al.2018).

For the Holocene, we model GOC and AMOC returning to values similar to the modern ocean estimates of Talley (2013). Our Holocene result for Atlantic (Pacific–Indian) Southern Ocean biological export productivity, of 3.4 (2.4) mol C m−2 yr−1 (Fig. 8), falls within modern observations for the Southern Ocean of 0.5–6 mol C m−2 yr−1 (e.g. Lourey and Trull2001; Weeding and Trull2004; Ebersbach et al.2011; Jacquet et al.2011; Cassar et al.2015; Arteaga et al.2019). Our model–data experiment results also reproduce values that fall within 1 standard deviation of the mean value in nearly all model boxes, for all of the atmosphere and ocean proxies in each MIS (Figs. 911).

Kohfeld and Chase (2017) suggested that sequential falls in atmospheric CO2 concentration were first the result of temperature, sea-ice cover, and potentially sea-ice-cover-induced Atlantic Southern Ocean “barrier mechanisms” or shallow stratification during MIS 5d and, second, followed by falls in deep Atlantic Ocean circulation and potentially dust-driven Southern Ocean biological productivity during MIS 4. Finally, a synthesis of those factors including enhanced Southern Ocean biology delivered the LGM atmospheric CO2 minimum. Our model–data results mostly agree with the Kohfeld and Chase (2017) hypothesis for glacial–interglacial CO2, particularly with regard to lower SST early in the glacial inception followed by weaker deep Atlantic Ocean circulation and stronger Southern Ocean biological export productivity later in the glacial cycle. However, we also posit a role for slowing GOC and no direct role for increased sea-ice cover in delivering lower atmospheric CO2 at the last glacial inception. Stephens and Keeling (2000) proposed that expanded glacial sea-ice cover around Antarctica could deliver LGM CO2 changes on its own, as a result of reduced air–sea gas exchange or in combination with ice-driven ocean stratification. However, Köhler et al. (2010) demonstrated with a carbon cycle box model that increased sea-ice cover leads to increased atmospheric CO2, due to less in-gassing of CO2 into the cold waters surrounding Antarctica. Kohfeld and Ridgwell (2009) reviewed estimates of the effects of decreased sea-ice cover at the last glacial termination and found a best estimate of −5 ppm within a range of −14–0 ppm, which is in the opposite direction to that envisaged by Stephens and Keeling (2000) and Kohfeld and Chase (2017). The modelling work by Stephens and Keeling (2000) was discounted by Kohfeld and Ridgwell (2009) because it assumed nearly all ocean degassing of CO2 was confined to the polar Antarctic region, when modern observations suggest the locus of outgassing is in the equatorial ocean (Takahashi et al.2003). In SCP-M, the effects of polar Southern Ocean sea-ice cover, modelled as a slowing down in air–sea gas exchange in the polar Southern Ocean surface box, are modest. This modelling result reflects the offsetting effects of upwelled nutrient- (and carbon) rich waters (degassing and higher CO2) against the effects of lower temperatures and enhanced biological export productivity (in-gassing and lower CO2). This finding may reflect our approach to treat Southern Ocean sea-ice cover simply as a regulator of the rate of air–sea gas exchange. Our approach may neglect other effects of sea-ice cover including as a contributor to changes in Southern Ocean brine formation, buoyancy forcing, upwelling, mixing, deep-ocean stratification and NADW formation rates (Morrison et al.2011; Brovkin et al.2012; Ferrari et al.2014; Kohfeld and Chase2017; Jansen2017; Marzocchi and Jansen2017). For example, Brovkin et al. (2012) found that in the CLIMBER-2 model, atmospheric CO2 was more sensitive to sea-ice cover when it was linked to weakened vertical diffusivity in the Southern Ocean of tracers such as DIC, thereby reducing outgassing of CO2. The synergistic effects of increased Antarctic Southern Ocean sea-ice cover discussed by Kohfeld and Chase (2017), in terms of reduced ocean vertical mixing rates to deliver reductions in atmospheric CO2, could be tested with a more complex model than SCP-M.

In addition to lower SST, increased-sea-ice cover and the other model forcings (Fig. 2), SCP-M requires additional changes in the ocean to deliver the  25 ppm fall in average CO2 concentration during MIS 5d and satisfy the other atmospheric and ocean proxy data. We model a weakening in GOC of  7 Sv during MIS 5d and further weakening until the LGM, a substantial change in the global ocean and not just the Atlantic basin. This underscores the importance of the global ocean in any hypothesis for the last glacial–interglacial cycle or LGM–Holocene (Fig. 8). We note that our simplified representation of GOC, as per Talley (2013), includes features that may be separated out or characterized differently in other models or hypotheses, such as AABW formation rate, Southern Ocean upwelling or shallow mixing or stratification, Pacific and Indian deepwater formation (PDW/IDW), or northward extension of AABW versus NADW formation of abyssal waters in the Atlantic Ocean (e.g. Menviel et al.2016; Kohfeld and Chase2017).

The period MIS 5e–5d does not feature in some oceanographic theories of glacial inception atmospheric CO2 decline, largely due to a focus on Atlantic Ocean data and a lack of any obvious changes in the Atlantic shallow–deep–abyssal proxy offsets at that period, as observed clearly during MIS 4 and the LGM (e.g. Oliver et al.2010; Yu et al.2016; Kohfeld and Chase2017). However, Govin et al. (2009) proposed an expansion of AABW across the Southern Ocean and weakening of circumpolar deep-water upwelling during MIS 5d, based on deep-ocean δ13C from the Atlantic and Indian basins. The proxy evidence of Govin et al. (2009) supports the concept of De Boer and Hogg (2014) that the glacial ocean could have exhibited slower and at the same time more expansive formation of AABW. Ganopolski et al. (2010) and Brovkin et al. (2012) modelled cooling SST and substitution of NADW by denser waters of Antarctic origin in the abyssal ocean as the main drivers of falling atmospheric CO2 at the last glacial inception. Menviel et al. (2012) modelled a transient slowdown in the rate of overturning circulation in the North Atlantic across MIS 5e–5d. Despite these findings, changes in ocean circulation at the last glacial inception are not obvious in Atlantic Ocean δ13C proxy data (Oliver et al.2010; Kohfeld and Chase2017).

To illustrate the plausibility of a slowdown in GOC during the last glacial inception in the context of deep-ocean δ13C proxy data, we show a model experiment testing the sensitivity of atmospheric CO2 and abyssal ocean δ13C to slowed GOC under MIS 5e and MIS 5d conditions (Fig. 13). Shown for comparison are the standard deviation of data values for abyssal ocean δ13C for MIS 5e (Fig. 13b). The experiment shows that slowing GOC from the MIS 5e model–data optimized value of 29 Sv (e.g. Fig. 8) delivers lower values for atmospheric CO2 (Fig. 13a) and more negative abyssal Pacific–Indian δ13C (Fig. 13b). However, in the experiment of decreasing GOC, modelled atmospheric CO2 crosses the  25 ppm change of the MIS 5e–5d transition well before the model's abyssal Pacific–Indian box δ13C breaches 1 standard deviation of the abyssal Pacific–Indian δ13C data for MIS 5e (Fig. 13b). Changes in the deep–abyssal δ13C offsets are also muted (Fig. 13c) relative to atmospheric CO2 and particularly for the Atlantic Ocean. The observation is even more obvious when including other ocean changes for the MIS 5e–5d transition, such as SST. When these changes are incorporated (shown as the “x” symbols in Fig. 13a and b), the atmospheric CO2 change across MIS 5e–5d is even more quickly satisfied by the modelled reduction in GOC, while abyssal ocean δ13C remains near its MIS 5e box average and well within 1 standard deviation. Despite a range of GOC variation that surpasses the MIS 5e–5d atmospheric CO2 reduction, the abyssal Atlantic δ13C result hardly varies, a particularly interesting finding. In SCP-M this can be explained by a reduced rate of AABW formation as a part of slowing GOC, leading to relatively greater influence of other Atlantic Ocean processes such as the deep–abyssal mixing and AMOC, which mixes deep water with a more positive δ13C into the abyssal Atlantic and offsets the effects of slowing GOC. Slowing GOC by itself leads to a more negative abyssal δ13C, as per the Pacific–Indian basin results. This type of dynamic could help explain why hypothesized or modelled changes in the ocean during the last glacial inception (e.g. Govin et al.2009; Menviel et al.2012; Brovkin et al.2012) do not show up more obviously in the deep and abyssal Atlantic Ocean δ13C proxy data (Oliver et al.2010; Kohfeld and Chase2017).

These observations from Fig. 13 could be exaggerated in SCP-M due to the large size of its ocean boxes and therefore relatively large spread of δ13C values and standard deviations for each box. In addition, this experiment may reflect idiosyncrasies in the SCP-M model design and its simple parameterization of ocean circulation and mixing. A finer-resolution model may show a greater sensitivity of the ocean box δ13C to variations in ocean circulation. Menviel et al. (2015) analysed the sensitivity of ocean and atmospheric δ13C to variations in NADW, AABW and North Pacific Deep Water (NPDW) formation rates in the context of past changes in atmospheric δ13C and CO2. Their modelling, using the more spatially detailed LOVECLIM and Bern3D models, showed modest but location-dependent sensitivities of ocean δ13C to slowing ocean circulation and particular sensitivity to AABW. These models are higher resolution and show greater sensitivity of δ13C to ocean circulation over depth intervals not differentiated in the SCP-M boxes. However, our simple experiment illustrated in Fig. 13 does highlight the potential for important changes in the ocean during glacial–interglacial periods to go unnoticed when focussed on one set of ocean proxy data and without validation by modelling.

As shown in Fig. 13, analysing Atlantic Ocean data in isolation and only qualitatively assessing ocean proxy data offsets (e.g. solely relying on standard deviations) may obscure features that could have contributed meaningfully to glacial falls in atmospheric CO2 (e.g. GOC). According to Talley (2013), GOC is a key part of the global ocean carbon cycle, operating in the Atlantic, Pacific and Indian ocean basins. Given it is a global feature, spread across all basins, its global changes may not show up as dramatic changes in proxy data offsets in any particular basin, despite it exerting a strong influence on atmospheric CO2. A number of authors highlight changes in Δ14C distributions in the Pacific Ocean during the LGM and Holocene, providing qualitative evidence of changes in ocean circulation in this basin and of it being a potential driver for post-glacial increase in atmospheric CO2 concentration (e.g. Sikes et al.2000; Marchitto et al.2007; Stott et al.2009; Cook and Keigwin2015; Skinner et al.2015; Ronge et al.2016; Skinner et al.2017). Ocean Δ14C values are particularly sensitive to ocean circulation rates (Broecker et al.1980). However, Δ14C proxy records in periods prior to the LGM and Holocene are sparse because they can only extend to  50 ka due to their radioactive decay in nature; therefore they cannot be applied to the glacial inception period.

Figure 13Sensitivity of atmospheric CO2 concentration and ocean δ13C to a downward variation in global ocean circulation parameter Ψ1 in MIS 5e in SCP-M. The x axis shows the range of variation in Ψ1 in Sv and the y axes show the model results for (a) atmospheric CO2 and (b) abyssal ocean δ13C in each basin. Shaded areas are the ± standard deviations for abyssal δ13C in MIS 5e. Panel (c) shows the deep–abyssal δ13C offset for each basin. Atmospheric CO2 in MIS 5e and 5d is shown for reference. The “x” symbols in (a) and (b) show the same experiment including other changes in the ocean across MIS 5e–5d: SST, salinity, Antarctic sea-ice cover, ocean volume and coral reef carbonate production. Southern Ocean biological export productivity is not varied in this experiment.


There is qualitative multi-proxy evidence for a slowdown or shoaling of AMOC during MIS 4. Kohfeld and Chase (2017) evaluated Atlantic basin δ13C data and surmised that Atlantic deep-ocean circulation slowed or shoaled during MIS 4. Yu et al. (2016) and Chalk et al. (2019) came to similar conclusions from the analysis of carbonate proxy records. Piotrowski et al. (2009) further suggested a reduced proportion of AMOC-sourced waters in the deep Indian Ocean during MIS 4, as deduced from Indian Ocean δ13C data. Our model–data results corroborate these findings, with a pronounced weakening in AMOC during MIS 4. SCP-M does not take explicit account of AMOC shoaling due to its rigid box boundaries, and therefore the change in proxy data between MIS 5a and 4 is resolved as weakening AMOC, which could understate the importance of this event. We also model a drop in AMOC during MIS 5b, which replicates abyssal Atlantic δ13C and CO32- observations (Figs. 5 and 7) and also accompanies a transient fall in atmospheric CO2 of 14 ppm during that period (Fig. 4). Menviel et al. (2012) modelled a transient but more dramatic decline in AMOC rate during MIS 5b and a more protracted but similarly large decline during MIS 4 (also modelled by Ganopolski et al.2010), in addition to a deepening in the remineralization depth of organic carbon.

Our model–data results indicate a role for increased Southern Ocean biological export productivity in achieving glacial troughs in atmospheric CO2 concentration during MIS 4 and MIS 2. Our finding of increased biological productivity, while constrained to MIS 4 and MIS 2, and a modest contributor to the overall glacial atmospheric CO2 drawdown aligns with proxy data for increased iron-rich continental dust supply to the Southern Ocean in these periods (e.g. Martinez-Garcia et al.2014; Lambert et al.2015; Kohfeld and Chase2017) and recent model–data exercises (e.g. Menviel et al.2016; Muglia et al.2018; Khatiwala et al.2019). Martin (1990) pioneered the “iron hypothesis”, which invoked the increased supply of continent-borne dusts to the Southern Ocean in glacial periods. Increased dust supply stimulated more plankton productivity where plankton was bio-limited in nutrients supplied in the dust, such as iron (Martin1990). Since then, the iron hypothesis has retained an important place in the debate over glacial–interglacial cycles of CO2. Watson et al. (2000) took experimental data on the effects of iron supply on plankton productivity in the Southern Ocean (Boyd2000) and applied this to a carbon cycle model across glacial–interglacial cycles. Their modelling, informed by the ocean experiment data, suggested that variations in the Southern Ocean iron supply and plankton productivity could account for large ( 40 ppm) swings in atmospheric CO2, with peak activity in the last glacial cycle during MIS 4 and MIS 2. Debate has continued over the magnitude of the contribution of Southern Ocean biological productivity to the glacial CO2 drawdown. According to Kohfeld et al. (2005), based on sediment data, enhanced Southern Ocean biological productivity could account for no more than half of the glacial CO2 drawdown. Others emphasize that Southern Ocean biological export productivity fluxes may have been weaker in the LGM in absolute terms but that with weaker Southern Ocean upwelling, the iron-enhanced productivity contributed to a stronger biological pump of carbon and was a major contributor to the LGM CO2 drawdown (Jaccard et al.2013; Martinez-Garcia et al.2014; Yamamoto et al.2019). Importantly, our finding for increased biological export productivity during MIS 4 and 2 is delivered without any model-simulated iron dust fertilization of the Southern Ocean and entirely on account of model results best-fit to the atmospheric and ocean proxy data used. Therefore the finding is a robust independently derived support for increased biological export productivity during MIS 4 and in particular MIS 2. It is important to note our model–data experiments assume unchanged biological export productivity in surface boxes outside of the Atlantic and Pacific–Indian subpolar Southern Ocean boxes across the last glacial–interglacial period. Some authors posit that low-latitude biological export productivity may have been stronger at the LGM due to increased shelf-sourced phosphorus (Broecker1981, 1982; Filippelli et al.2007; Tamburini and Föllmi2009; Menviel et al.2012) or increased biological matter remineralization depth (Matsumoto2007; Menviel et al.2012). Others argue that low-latitude biological export productivity was weaker at the LGM due to less upwelling of thermocline waters and lower shallow-ocean nutrient levels (Calvo et al.2011; Hayes et al.2011; Winckler et al.2016). Weaker (stronger) glacial biological export productivity in the low-latitude surface boxes would reduce (increase) the sensitivity of atmospheric CO2 to ocean circulation in our model–data experiments.

5.2 Contribution and attribution analysis

Table 3 shows a contribution analysis for the data observations in each MIS model–data optimization of ocean parameter values. The ranking is based on the relative standard deviation (RSD) for each proxy data observation (or set of data observations) in each MIS, with the highest ranking (e.g. 1) given to the data observation with the lowest RSD in each model box or MIS. The contribution analysis shows that atmospheric δ13C and CO2 exert the greatest influence on the optimization results throughout the MIS experiments. This reflects that each of these atmospheric data time series is derived from a single source and does not require locational averaging as in the ocean boxes. For the atmosphere data, only MIS averaging (not model box dimension) takes place, and therefore there is a lower standard deviation of the data in most MIS time slices. For the ocean boxes, averaging on depth and latitude takes place as well as MIS averaging to derive a box or MIS mean data value. Using a box model with large boxes such as SCP-M means that large parts of the ocean are averaged into the ocean box mean value, and therefore there is an increased spread of data values around the mean for those boxes. Therefore, the model–data results show a precise fit to the atmospheric δ13C and CO2 data as shown in Figs. 911. The results for oceanic variables are typically less precise but also fall within the standard deviations of the data observations for each box and MIS (Figs. 911). Others have attempted glacial–interglacial model–data studies focusing only on the ocean data without matching atmospheric data (e.g. LeGrand and Wunsch1995; Gebbie and Huybers2006; Hesse et al.2011; Zhao et al.2017; Kurahashi-Nakamura et al.2017). While these studies could potentially elicit more detail on oceanic processes, they are also potentially fraught due to the high spread of data values for the oceanic data and could return results that are not consistent with the relatively well-constrained glacial–interglacial atmosphere data. For our study, the express purpose is to identify causes of changes in atmospheric CO2 concentration, so it is appropriate that atmospheric data observations make an important contribution to the model results. However, as shown in Figs. 911 this is not at the expense of providing plausible results for the ocean variables. Additional parameter sensitivity analysis of the model–data experiments is shown in Fig. S2.

Table 3Contribution of proxy data observations to the model–data experiment results for ocean parameter values in each MIS. Each proxy data observation from each model box is ranked from 1 to 6 in each MIS based on the relative standard deviation (RSD) of its data points. A ranking of 1 is given to the data observation with the smallest RSD in each MIS. A smaller RSD gives the data observation a higher weighting in the model–data optimization and therefore a greater contribution to the model results. Δ14C proxy data do not exist for periods before MIS 3.

NA – not available.

Download Print Version | Download XLSX

Figure 14 shows the contribution to the glacial drawdown in atmospheric CO2 by each mechanism we modelled, relative to the last interglacial period (MIS 5e), in SCP-M. Our model–data study finds that approximately half of the glacial atmospheric CO2 drawdown is contributed by weakened ocean circulation (GOC and AMOC), with the other half contributed by a combination of lower SST, increased Southern Ocean biological export productivity, varying coral reef carbonate production and dissolution, and increased polar Southern Ocean sea-ice cover. Weakened GOC delivers the highest contribution to falling CO2, followed by lower SST, weakened AMOC and stronger Southern Ocean biological export productivity. Lower SST leads to modest reductions in atmospheric CO2 concentration early in the glacial cycle, increasing as the ocean cools further during MIS 4, and is an important contributor to decreased CO2 in the LGM (Kohfeld and Chase2017). Some studies observed that early versions of box models tended to overstate the effects of SST and other processes at high latitudes on atmospheric CO2 relative to general circulation models (GCMs) (Broecker et al.1999; Archer et al.2000; Ridgwell2001; Kohfeld and Ridgwell2009). However, our modelled estimate of 28 ppm for the contribution of SST to the glacial–interglacial atmospheric CO2 change (Fig. 14) falls within the range of GCM-derived estimates of 21–30 ppm (mean value 26 ppm) compiled by Kohfeld and Ridgwell (2009), is similar to that of Menviel et al. (2012) (27.5 ppm) and substantially less than another recent GCM-derived estimate of 44 ppm (Khatiwala et al.2019). Southern Ocean biological export productivity strengthens during MIS 4 and contributes a peak of −13 ppm by MIS 2 (LGM).

The smaller glacial terrestrial biosphere contributes 13 ppm CO2 to the atmosphere during the LGM (MIS 2), consistent with other modelled estimates (Köhler et al.2010; Menviel et al.2012; Ganopolski and Brovkin2017). Other parameters contribute lesser increases in CO2 (salinity, ocean volume) and decreases (Antarctic sea-ice, coral reefs) during the glacial cycle. Our estimate for coral reefs of −9 ppm CO2 is at the lower end of the range of 6–20 ppm summarized in Kohfeld and Ridgwell (2009), suggesting that our simple parameterization of the coral reef carbon and alkalinity fluxes could underestimate its effect, likely due to the assumed fast mixing rates of reef carbon and alkalinity into the surface boxes in SCP-M. Ridgwell et al. (2003) modelled +20 ppm atmospheric CO2 from coral reef carbonate accumulation in the Holocene period, noting a high sensitivity of their model to coral reef accumulation rates. It is likely that our model–data results underestimate the contribution of AMOC because our model does not explicitly resolve AMOC shoaling (e.g. Menviel et al.2012; Brovkin et al.2012; Yu et al.2016; Eggleston et al.2016; Kohfeld and Chase2017; Menviel et al.2020), other than a linear–positive linkage between the AMOC parameter and a deep–abyssal Atlantic box mixing term (less mixing between the deep and abyssal Atlantic boxes as AMOC slows). Therefore, the analysis could miss additional features of the AMOC mechanism which could contribute to greater atmospheric CO2 drawdown in Fig. 14. The contribution of the model parameters to the glacial atmospheric CO2 drawdown shown in Fig. 14 incorporates the effects of various feedbacks in the model such as continental weathering and calcium carbonate compensation.

Figure 14Impacts on atmospheric CO2 concentration of model parameters in the model–data experiment results, from the last interglacial period (MIS 5e) to the Last Glacial Maximum (MIS 2). SST: sea surface temperature; ReefC: shallow carbonate production or dissolution; GOC: global ocean circulation; AMOC: Atlantic meridional overturning circulation; SO Bio Export: Southern Ocean biological export productivity.


5.3 The LGM and Holocene

Within the context of LGM–Holocene studies, our findings corroborate the hypothesis that a number of mechanisms, not one singular factor, delivered the  85 ppm increase in atmospheric CO2 concentration from the LGM to the Holocene (e.g. Kohfeld and Ridgwell2009; Köhler et al.2010; Sigman et al.2010; Hain et al.2010; Menviel et al.2012; Brovkin et al.2012; Ferrari et al.2014; Menviel et al.2016; Ganopolski and Brovkin2017; Kohfeld and Chase2017; Muglia et al.2018). This finding is more obvious when the sequential nature of changes is observed over the full glacial–interglacial cycle, as distinct from analysing the LGM and Holocene in isolation. Our model–data results agree with those of Menviel et al. (2016), who showed that the LGM oceanic δ13C and Δ14C records were most consistent with a weak GOC and AMOC. Menviel et al. (2016) further showed that this weak oceanic circulation would significantly increase the deep-ocean carbon content (and thus significantly contribute to the pCO2 decrease). The longer timescale of our analysis highlights that changes in GOC and AMOC took place earlier in the glacial cycle than the LGM, and were at or near their glacial minima prior to the LGM. However, some caution is required as our model–data results reflect the mean MIS state. For example, Menviel et al. (2014) modelled substantial variability in AMOC during Dansgaard–Oeschger events during MIS 3. Such variability is averaged out in our MIS state experiments. Our model–data results also support increased Southern Ocean biological export productivity in the LGM as a contributor to the lower LGM atmospheric CO2 concentration (and in MIS 4) as well as lower SST and to a lesser extent coral reef carbonate dissolution.

5.4 The terrestrial biosphere

Our modelled increase in the terrestrial biosphere carbon stock from the LGM to Holocene, of  600 Pg C (Fig. 12), falls within but towards the upper end of recent estimates of this change of 300–850 Pg C (e.g. Joos et al.2004; Brovkin et al.2007; Köhler et al.2010; Prentice et al.2011; Brovkin et al.2012; Ciais et al.2012; Peterson et al.2014; Menviel et al.2016; Jeltsch-Thommes et al.2019). Brovkin et al. (2007), Brovkin et al. (2012) and Köhler et al. (2010) all modelled a  500–550 Pg C increase in the terrestrial biosphere between the LGM and Holocene (Prentice et al. (2011) estimated 550–694 Pg C). According to François et al. (1999), palynological and sediment data infer that the terrestrial biosphere carbon stock was 700–1350 Pg C smaller in the LGM than the present. Ciais et al. (2012) pointed to a growth of a large carbon pool in steppes and tundra during the LGM as an offsetting feature to the declining tropical biosphere, leading to a smaller estimate of  330 Pg C (Ganopolski and Brovkin (2017) modelled a similar estimate of 350 Pg C). Jeltsch-Thommes et al. (2019) estimated a glacial–interglacial change in terrestrial biosphere of 850 Pg C (median estimate; range 450 to 1250 Pg C), a similar estimate to that of Joos et al. (2004) of 820–850 Pg C. Jeltsch-Thommes et al. (2019) demonstrated the importance of including ocean-sediment and weathering fluxes in their modelling estimates, and suggested other studies may underestimate the full deglacial change in the terrestrial biosphere carbon stock. While our model results ( 600 Pg C) are higher than some estimates of the LGM–Holocene change in the terrestrial biosphere (e.g. Ciais et al.2012; Menviel et al.2016; Ganopolski and Brovkin2017), they are mostly in good agreement (e.g. Joos et al.2004; Brovkin et al.2007; Köhler et al.2010; Prentice et al.2011; Brovkin et al.2012; Peterson et al.2014; Jeltsch-Thommes et al.2019), and our NPP estimates mostly align with the glacial–interglacial cycle NPP reconstruction of Hoogakker et al. (2016) as shown in Fig. 12. The driver for NPP in the simple terrestrial biosphere module in SCP-M is atmospheric CO2 concentration via carbon fertilization (e.g. Otto et al.2002; Kaplan et al.2002; Joos et al.2004; Hoogakker et al.2016). Temperature and precipitation also exert important controls on NPP (e.g. François et al.1999; van der Sleen et al.2015), which are not accounted for in our model.

The isotopic fractionation behaviour of the terrestrial biosphere may also vary in glacial–interglacial time frames. This has been studied for the LGM, Holocene and the present day (e.g. Collatz et al.1998; François et al.1999; Kaplan et al.2002; Köhler and Fischer2004; Joos et al.2004; Kohn2016). The variation in isotopic fractionation within the terrestrial biosphere reflects changes in the relative proportions of plants with the C3 and C4 photosynthetic pathways but also strong variations within the same photosynthetic pathways themselves (François et al.1999; Kohn2010; Schubert and Jahren2012; Kohn2016). The drivers for these changes include relative sea level and exposed land surface area (François et al.1999), global tree-line extent (Köhler and Fischer2004), atmospheric temperature and CO2 (Collatz et al.1998; François et al.1999; Köhler and Fischer2004; Kohn2010; Schubert and Jahren2012), global and localized precipitation and humidity (Huang et al.2001; Kohn2010; Schubert and Jahren2012; Kohn2016), and also changes in the intercellular CO2 pressure in the leaves of C3 plants (François et al.1999). Estimated changes in average terrestrial biosphere δ13C signature between the LGM and the Holocene fall in the range -0.3-1.8 (less negative δ13C signature in the LGM), with further changes estimated from the onset of the Holocene to the pre-industrial and even greater changes to the present day (due to rising atmospheric CO2). This feature has been covered in detail in studies that focussed on the terrestrial biosphere between the LGM and Holocene but less so in modelling and model–data studies of the last glacial–interglacial cycle. Menviel et al. (2016) provided a sensitivity of −0.7 ‰ and +0.5 ‰ around an average LGM terrestrial biosphere δ13C value of −23.3‰, based on previous modelling of the LGM–Holocene time frame by Joos et al. (2004). Another modelling study (Menviel and Joos2012) assessed the variation in LGM–Holocene δ13C of the terrestrial biosphere to be a minor factor and it was not considered. Köhler and Fischer (2004) assessed the changing δ13C signature of plants between the LGM and Holocene to be a minor factor in setting δ13C of marine DIC, compared to changes in the absolute size of the terrestrial biosphere across this period. Given the uncertainty and ranges of starting estimates of terrestrial biosphere δ13C (for example, the very large range in present-day estimates of C3 plant δ13C; Kohn2010, 2016), the uncertain LGM–Holocene changes, the large number of potential drivers of relative C3 and C4, and the further uncertainty in extrapolating the posited LGM–Holocene changes back for the preceding 100 kyr and finally the modest changes relative to the average δ13C signature, we omit this feature with the caveat that there is added uncertainty in our terrestrial biosphere results with respect to the δ13C signature applied. Our choice of a constant terrestrial biosphere δ13C signature of −23‰ is similar to values assumed by Menviel et al. (2016) and Jeltsch-Thommes et al. (2019) (−23.3‰, −24‰ respectively) but more negative than assumed in Brovkin et al. (2002), Köhler and Fischer (2004), and Joos et al. (2004) (−16‰, −17‰). In summary, our aim is not to contribute new findings of the terrestrial biosphere but to ensure that the simple representation of the terrestrial biosphere in SCP-M provides the appropriate feedbacks to our (exhaustive) glacial–interglacial cycle model–data optimization experiments that are in line with the published estimates discussed above.

5.5 Advantages and limitations of this study

The use of a simple box model for this model–data study, SCP-M, enabled a range of proxies to be incorporated into MIS data reconstructions and a large number of simulations ( 9000 in each MIS) to explore possible parameter combinations in each MIS. However, the use of a simple box model means that some details are lost in the analysis. Given the large spatial coverage of the SCP-M boxes, data for large areas of the ocean are averaged. In the case of a carbonate ion proxy, we apply a default estimate of standard deviation to account for the large volume of ocean covered by SCP-M's boxes relative to the proxy data locations and to enable the normalization of the carbonate ion proxy data in a procedure that uses the data standard deviation as a weighting. Despite this caveat, we believe that the model–data experiment results provide a good match to the data across the various atmospheric and ocean proxies as shown in Figs. 911.

Most major processes in the SCP-M model are simply parameterized, allowing them to be free-floated in model–data experiments. However, the driving factors behind parameter value changes can only be speculated on. For example, slowdown in GOC may be the result of changing wind patterns or buoyancy fluxes around Antarctica (Morrison and Hogg2013), Antarctic sea-ice cover (Ferrari et al.2014) or shoaling AMOC leading to extensive filling of the abyssal ocean by waters sourced from GOC (Curry and Oppo2005; De Boer and Hogg2014; Jansen2017). Probing the root cause of our model–data findings would require a more detailed physical and/or biogeochemical model. Furthermore, we apply a simple representation of the terrestrial biosphere in our model–data experiments, relying primarily on atmospheric CO2 as the driver for NPP. This approach provided reasonable results for the terrestrial biosphere carbon stock and NPP, on the whole, but may miss some detail in the terrestrial biosphere during the last glacial–interglacial cycle. Our MIS time slicing obscures details in the proxy records within MIS. For example, Yu et al. (2013) observed a transient drop in carbonate ion concentrations in the deep Pacific Ocean during MIS 4, and there are large transient changes in atmospheric δ13C during MIS 4 and MIS 3. Ganopolski et al. (2010) and Menviel et al. (2012) modelled transient collapses and rebounds in AMOC during MIS 4 (and other short-term changes in atmospheric dust supply and depth of biological nutrient remineralization), which could have contributed to the full observed magnitude of changes in atmospheric δ13C across this period (e.g. Eggleston et al.2016) – not captured with our MIS-averaging approach. We omitted the last glacial termination from our analysis, a period in which atmospheric CO2 concentration increased by  85 ppm in 8 kyr. Future model–data optimization work could probe this period at 1 kyr intervals or with transient, data-optimized simulations to profile the unwinding of processes that led to the last glacial cycle CO2 drawdown.

In summary, while the carbon cycle box model we applied is high level in nature and there are caveats, the modelling itself is heavily constrained by natural observations and proxy data from the carbon cycle. Therefore, this work presents a plausible set of modelled outcomes for the last glacial–interglacial cycle.

6 Conclusions

Multiple processes drove changes in atmospheric CO2 concentration during the last glacial–interglacial cycle. Against a backdrop of varied SST, salinity, sea-ice cover, ocean volume and reef carbonates, we modelled sequentially weaker GOC (first) and AMOC (second) to reduce atmospheric CO2 concentration in the lead-up to the LGM. During the LGM, increased Southern Ocean biological export productivity delivered an incremental fall in atmospheric CO2 concentration, resulting in the glacial cycle CO2 minimum. GOC, AMOC, Southern Ocean biology and SST rebounded to modern values between the LGM and Holocene, contributing to the sharp post-glacial increase in atmospheric CO2 concentration. The terrestrial biosphere played an important negative feedback role during the glacial–interglacial cycle, releasing δ13C-depleted CO2 to the atmosphere at times during the glaciation and taking up CO2 during the termination and Holocene. These model–data results were achieved with a simple carbon cycle box optimized for proxy data for atmospheric CO2, atmospheric and ocean δ13C and Δ14C, and ocean CO32-. Our results agree with hypotheses for glacial–interglacial cycle CO2 that include varying ocean circulation, Southern Ocean biological export productivity and other physical and biogeochemical changes in the marine and terrestrial carbon cycle (e.g. Kohfeld and Ridgwell2009; Sigman et al.2010; Ganopolski et al.2010; Brovkin et al.2012; Menviel et al.2012; Ferrari et al.2014; Menviel et al.2016; Kohfeld and Chase2017; Ganopolski and Brovkin2017). We emphasize the need to include the Pacific and Indian oceans in an evaluation of the oceanic carbon cycle, particularly in relation to the last glacial–interglacial cycle and the LGM–Holocene transition.

Many uncertainties exist in the data and the prescribed nature of carbon cycle processes in a box model. However, such uncertainty is largely inescapable when dealing with models and proxy data. We propose these model–data results as one set of plausible results for the last glacial carbon cycle, in agreement with available proxy data, and see them as encouraging for the use of models and data to help constrain hypotheses for past changes in the Earth's carbon cycle.

Code and data availability

The model code, processed data files, model–data experiment results and any (published) raw proxy data gathered in the course of this work are located at (O'Neill et al.2021). No original data were created in this work nor were unpublished data used. Figure S3 contains an overview of the files contained in the repository. For more detail on the SCP-M equations, see O'Neill et al. (2019).


The supplement related to this article is available online at:

Author contributions

CO undertook model development work, data gathering, modelling and model–data experiments. AH provided the oceanographic interpretation and guided modelling and data analysis. ME designed model–data experiments, provided input into data analysis, and guided the modelling of the marine biology and carbon isotopes. BO contributed glacial–interglacial cycle model forcings and input to modelling of the reef carbonates. SE contributed to the modelling of the marine biology and carbonate pump. All authors contributed to drafting and reviewing the document.

Competing interests

The authors declare that they have no conflict of interest.


Stewart Fallon provided input to the processing of radiocarbon data. Malcolm Sambridge provided input on model–data optimization and inversions. Jimin Yu provided helpful discussions and published carbonate ion proxy data.

Review statement

This paper was edited by Laurie Menviel and reviewed by Pearse Buchanan and three anonymous referees.


Adkins, J., McIntyre, K., and Schrag, D.: The Salinity, Temperature, and δ18O of the Glacial Deep Ocean, Science, 298, 1769–1773, 2002. a, b, c, d, e, f, g, h

Ahn, J. and Brook, E.: Siple Dome ice reveals two modes of millennial CO2 change during the last ice age, Nat. Commun., 5, 3723,, 2014. a, b

Anderson, R., Chase, Z., Fleisher, M., and Sachs, J.: The Southern Ocean's biological pump during the Last Glacial Maximum, Deep Sea Res. Part II, 49, 1909–1938, 2002. a

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

Archer, D., Eshel, G., Winguth, A., Broecker, W., Pierrehumbert, R., Tobis, M., and Jacob, R.: Atmospheric CO2: sensitivity to the biologicalpump in the ocean, Global Biogeochem. Cy., 14, 1219–1230, 2000. a

Arteaga, L., Pahlow, M., Bushinsky, S., and Sarmiento, J.: Nutrient Controls on Export Production in the Southern Ocean, Global Biogeochem. Cy., 33, 942–956, 2019. a

Barker, S., Knorr, G., Vautravers, M., Diz, P., and Skinner, L.: Extreme deepening of the Atlantic overturning circulation during deglaciation, Nat. Geosci., 3, 567–571, 2010. a

Bereiter, B., Lüthi, D., Siegrist, M., Schüpbach, S., Stocker, T. F., and Fischer, H.: Mode change of millennial CO2 variability during the last glacial cycle associated with a bipolar marine carbon seesaw, P. Natl. Acad. Sci. USA, 109, 9755–9760, 2012. a, b

Bereiter, B., Eggleston, S., Schmitt, J., Nehrbass-Ahles, C., Stocker, T., Fischer, H., Kipfstuhl, S., and Chappellaz, J.: Revision of the EPICA Dome C CO2 record from 800 to 600kyr before present, Geophys. Res. Lett, 42, 542–549, 2015. a

Berger, W.: Increase of carbon dioxide in the atmosphere during deglaciation: The coral reef hypothesis, Naturwissenschaften, 69, 87–88, 1982. a

Boyd, P., Watson, A., Law, C., Abraham, E., Trull, T., Murdoch, R., Bakker, D., Bowie, A., Buesseler, K., Chang, H., Charette, M., Croot, P., Downing, K., Frew, R., Gall, M., Hadfield, M., Hall, J., Harvey, M., Jameson, G., LaRoche, J., Liddicoat, M., Ling, R., Maldonado, M., McKay, R., Nodder, S., Pickmere, S., Pridmore, R., Rintoul, S., Safi, K., Sutton, P., Strzepek, R., Tanneberger, K., Turner, S., Waite, A., and Zeldis, J.: A mesoscale phytoplankton bloom in the polar Southern Ocean stimulated by iron fertilization, Nature, 407, 695–702, 2000. a

Broecker, W.: Glacial to Interglacial Changes in Ocean and Atmosphere Chemistry, in: Climatic Variations and Variability: Facts and Theories, edited by: Berger, A., NATO Advanced Study Institutes Series (Series C – Mathematical and Physical Sciences), vol. 72, pp. 111–121, Springer, Dordrecht, 1981. a

Broecker, W., Yu, J., and Putnam, A.: Two contributors to the glacial CO2 decline, Earth Planet. Sci. Lett., 429, 191–196, 2015. a

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

Broecker, W. S. and Barker, S.: A 190 ‰ drop in atmosphere's Δ14C during the “Mystery Interval” (17.5 to 14.5 kyr), Earth Planet. Sci. Lett., 256, 90–99, 2007. a, b, c

Broecker, W. S., Peng, T. H., and Engh, R.: Modeling the carbon system, Radiocarbon, 22, 565–598, 1980. a

Broecker, W. S., Lynch-Stieglitz, J., Archer, D., Hofmann, M., Maier-Reimer, E., Marchal, O., Stocker, T., and Gruber, N.: How strong is the Harvardton-Bear constraint?, Global Biogeochem. Cy., 13, 817–820, 1999. a

Brovkin, V., Claussen, J. B. M., Ganopolski, A., Kubatzki, C., Petoukhov, V., and Andreev, A.: Carbon cycle, vegetation, and climate dynamics in the Holocene: Experiments with the CLIMBER-2 model, Global Biogeochem. Cy., 16, 1139,, 2002. 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, b, c

Brovkin, V., Ganopolski, A., Archer, D., and Munhoven, G.: Glacial CO2 cycle as a succession of key physical and biogeochemical processes, Clim. Past, 8, 251–264,, 2012. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r

Bryan, S., Marchitto, T., and Lehman, S.: The release of 14C-depleted carbon from the deep ocean during the last deglaciation: Evidence from the Arabian Sea, Earth Planet. Sci. Lett., 298, 244–254, 2010. a, b

Burke, A. and Robinson, L.: The Southern Ocean's Role in Carbon Exchange During the Last Deglaciation, Science, 335, 557–561, 2012. a, b, c, d, e, f

Calvo, E., Pelejero, C., Pena, L., Cacho, I., and Logan, G.: Eastern equatorial pacific productivity and related-CO2 changes since the last glacial period, P. Natl. Acad. Sci., 108, 5537–5541, 2011. a

Cassar, N., Wright, S., Thomson, P., Trull, T., Westwood, K., de Salas, M., Davidson, A., Pearce, I., Davies, D., and Matear, R.: The relation of mixed-layer carbon export production to plankton community in the Southern Ocean, Global Biogeochem. Cy., 29, 446–462, 2015. a

Chalk, T., Foster, G., and Wilson, P.: Dynamic storage of glacial CO2 in the Atlantic Ocean revealed by boron [CO32-] and pH records, Earth Planet. Sci. Lett., 510, 1–11, 2019. a, b

Chen, T., Robinson, L., Burke, A., Southon, J., Spooner, P., Morris, P., and Ng, H.: Synchronous centennial abrupt events in the ocean and atmosphere during the last deglaciation, Science, 349, 1537–1541, 2015. a

Ciais, P., Tagliabue, A., Cuntz, M., Bopp, L., Scholze, M., Hoffmann, G., Lourantou, A., Harrison, S. P., Prentice, I. C., Kelley, D. I., 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, b, c, d, e, f, g

Collatz, G., Berry, J., and Clark, J.: Effects of climate and atmospheric CO2 partial pressure on the global distribution of C4 grasses: present, past, and future, Oecologia, 114, 441–454, 1998. a, b

Cook, M. and Keigwin, L.: Radiocarbon profiles of the NW Pacific from the LGM and deglaciation: Evaluating ventilation metrics and the effect of uncertain surface reservoir ages, Paleoceanography, 30, 174–195, 2015. a

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

Davies-Walczak, M., Mix, A., Stoner, J., Southon, J., Cheseby, M., and Xuan, C.: Late Glacial to Holocene radiocarbon constraints on North Pacific Intermediate Water ventilation and deglacial atmospheric CO2 sources, Earth Planet. Sci. Lett., 397, 57–66, 2014. a, b

De Boer, A. M. and Hogg, A. M. C.: Control of the glacial carbon budget by topographically induced mixing, Geophys. Res. Lett, 41, 4277–4284, 2014. a, b, c

DeVries, T. and Weber, T.: The export and fate of organic matter in the ocean: New constraints from combining satellite and oceanographic tracer observations, Paleoceanography, 31, 535–555, 2017. a, b, c

Dunne, J. P., Armstrong, R. A., Gnanadesikan, A., and Sarmiento, J. L.: Empirical and mechanistic models for the particle export ratio, Global Biogeochem. Cy., 19, GB4026,, 2005. a, b, c

Ebersbach, F., Trull, W., Davies, D., and Bray, S.: Controls on mesopelagic particle fluxes in the Sub-Antarctic and Polar Frontal Zones in the Southern Ocean south of Australia in summer – Perspectives from free-drifting sediment traps, Deep Sea Res. Part II, 58, 2260–2276, 2011. a

Eggleston, S., Schmitt, J., Bereiter, B., Schneider, R., and Fischer, H.: Evolution of the stable carbon isotope composition of atmospheric CO2 over the last glacial cycle, Paleoceanography, 31, 434–452, 2016. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Elsig, J., Schmitt, J., Leuenberger, D., Schneider, R., Eyer, M., Leuenberger, M., Joos, F., Fischer, H., and Stocker, T. F.: Stable isotope constraints on Holocene carbon cycle changes from an Antarctic ice core,, Nature, 461, 507–510, 2009. a, b

Ferrari, R., Jansen, M., Adkins, J., Burke, A., Stewart, A. L., and Thompson, A.: Antarctic sea ice control on ocean circulation in present and glacial climates, P. Natl. Acad. Sci., 111, 8753–8758, 2014. a, b, c, d, e, f

Filippelli, G., Latimer, J., Murray, R., and Flores, J.-A.: Productivity records from the Southern Ocean and the equatorial Pacific Ocean: Testing the glacial Shelf-Nutrient Hypothesis, Deep Sea Res. Part II, 54, 2443–2452, 2007. a

Follows, M. J., Ito, T., and Dutkiewicz, S.: On the solution of the carbonate chemistry system in ocean biogeochemistry models, Ocean Modell., 12, 290–30, 2006. a

François, L., Godderis, Y., Warnant, P., Ramstein, G., de Noblet, N., and Lorenz, S.: Carbon stocks and isotopic budgets of the terrestrial biosphere at mid-Holocene and last glacial maximum times, Chem. Geol., 159, 163–199, 1999. a, b, c, d, e, f, g, h, i

Freeman, E., Skinner, L., Waelbroeck, C., and Hodell, D.: Radiocarbon evidence for enhanced respired carbon storage in the Atlantic at the Last Glacial Maximum, Nat. Commun., 7, 1–8,, 2016. a, b

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, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s

Ganopolski, A., Calov, R., and Claussen, M.: Simulation of the last glacial cycle with a coupled climate ice-sheet model of intermediate complexity, Clim. Past, 6, 229–244,, 2010. a, b, c, d, e, f, g, h, i, j

Gebbie, G. and Huybers, P.: Meridional circulation during the Last Glacial Maximum explored through a combination of South Atlantic δ18O observations and a geostrophic inverse model, Geochem. Geophys. Geosys, 7, Q11N07,, 2006. a

Govin, A., Michel, E., Labeyrie, L., Waelbroeck, C., Dewilde, F., and Jansen, E.: Evidence for northward expansion of Antarctic Bottom Water mass in the Southern Ocean during the last glacial inception, Paleoceanography, 24,, 2009. a, b, c, d, e, f, g, h

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, b, c

Harman, I., Trudinger, C., and Raupach, M.: SCCM – the Simple Carbon-Climate Model: Technical Documentation, CAWCR Technical Report 047, CSIRO Centre for Australian Weather and Climate Research, CSIRO Marine and Atmospheric Research, FC Pye Laboratory, GPO Box 3023, Canberra, ACT, 2601, Australia, 2011. a, b, c

Harrison, K. G.: Role of increased marine silica input on paleo-pCO2 levels, Paleoceanography, 15, 292–298, 2000. a

Hayes, C., Anderson, R. F., and Fleisher, M. Q.: Opal accumulation rates in the equatorial Pacific and mechanisms of deglaciation, Paleoceanography, 26, PA1207,, 2011. a

Henson, S. A., Sanders, R., Madsen, E., Morris, P. J., Moigne, F. L., and Quartly, G. D.: A reduced estimate of the strength of the ocean's biological carbon pump, Geophys. Res. Lett., 38, L04606,, 2011. a, b, c

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

Hines, S., Southon, J., and Adkins, J.: A high-resolution record of Southern Ocean intermediate water radiocarbon over the past 30,000 years, Earth Planet. Sci. Lett., 432, 46–48, 2015. a, b

Hogg, A. M.: Glacial cycles and carbon dioxide: A conceptual model, Geophys. Res. Lett., 35, L01701,, 2008. a, b, c

Hoogakker, B. A. A., Smith, R. S., Singarayer, J. S., Marchant, R., Prentice, I. C., Allen, J. R. M., Anderson, R. S., Bhagwat, S. A., Behling, H., Borisova, O., Bush, M., Correa-Metrio, A., de Vernal, A., Finch, J. M., Fréchette, B., Lozano-Garcia, S., Gosling, W. D., Granoszewski, W., Grimm, E. C., Grüger, E., Hanselman, J., Harrison, S. P., Hill, T. R., Huntley, B., Jiménez-Moreno, G., Kershaw, P., Ledru, M.-P., Magri, D., McKenzie, M., Müller, U., Nakagawa, T., Novenko, E., Penny, D., Sadori, L., Scott, L., Stevenson, J., Valdes, P. J., Vandergoes, M., Velichko, A., Whitlock, C., and Tzedakis, C.: Terrestrial biosphere changes over the last 120 kyr, Clim. Past, 12, 51–73,, 2016. a, b, c, d, e, f, g, h

Huang, Y., Street-Perrott, F., Metcalfe, S., Brenner, M., Moreland, M., and Freeman, K.: Climate change as the dominant control on glacial-interglacial variations in C3 and C4 plant abundance, Science, 293, 1647–1651, 2001. a

Jaccard, S., Hayes, C., Martínez-García, A., an R.F. Anderson, D. H., Sigman, D., and Haug, G.: Two Modes of Change in Southern Ocean Productivity Over the Past Million Years, Science, 339, 1419–1423, 2013. a

Jacquet, S., Lam, P., Trull, T., and Dehairs, F.: Carbon export production in the Polar Front Zone and Subantarctic Zone south of Tasmania, Deep Sea Res. Part II, 58, 2277–2292, 2011. a

Jansen, M.: Glacial ocean circulation and stratification explained by reduced atmospheric temperature, P. Natl. Acad. Sci., 114, 45–50, 2017. a, b, c

Jeltsch-Thömmes, A., Battaglia, G., Cartapanis, O., Jaccard, S. L., and Joos, F.: Low terrestrial carbon storage at the Last Glacial Maximum: constraints from multi-proxy data, Clim. Past, 15, 849–879,, 2019. a, b, c, d, e

Joos, F., Gerber, S., Prentice, I. C., Otto-Bliesner, B., and Valdes, P.: Transient simulations of Holocene atmospheric carbon dioxide and terrestrial carbon since the Last Glacial Maximum, Global Biogeochem. Cy., 18, GB2002,, 2004. a, b, c, d, e, f, g, h

Kaplan, J., Prentice, I., Knorr, W., and Valdes, P.: Modeling the dynamics of terrestrial carbon storage since the Last Glacial Maximum, Geophys. Res. Lett., 22, 2074,, 2002. a, b, c

Key, R.: Ocean process tracers: Radiocarbon, in: Encyclopedia of Ocean Sciences, pp. 2338–2353, Academic Press, London, 2001. a, b, c

Key, R., Kozyr, A., Sabine, C. L., Lee, K., Wanninkhof, R., Bullister, J. L., Feely, R. A., Millero, F. J., Mordy, C., and Peng, T.-H.: A global ocean carbon climatology: Results from Global Data Analysis Project (GLODAP), Global Biogeochem. Cy., 18, GB4031,, 2004. a

Khatiwala, S., Schmittner, A., and Muglia, J.: Air-sea disequilibrium enhances ocean carbon storage during glacial periods, Sci. Adv., 5, 1–10,, 2019. a, b, c, d

Kleypas, J.: Modeled estimates of global reef habitat and carbonate production since the Last Glacial Maximum, Paleoceanogr. Paleoclim., 12, 533–545, 1997. a

Knox, F. and McElroy, M.: Changes in Atmospheric CO2: Influence of the Marine Biota at High Latitude, J. Geophys. Res., 89, 4269–4637, 1984. a, b, c

Kohfeld, K. and Chase, Z.: Temporal evolution of mechanisms controlling ocean carbon uptake during the last glacial cycle, Earth Planet. Sci. Lett., 472, 206–215, 2017. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa, ab, ac, ad, ae, af, ag, ah, ai, aj, ak

Kohfeld, K. and Ridgwell, A.: Glacial-Interglacial Variability in Atmospheric CO2, Surface Ocean–Lower Atmosphere Processes, Geophys. Res. Series, 187, 251–286, 2009. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q

Kohfeld, K., Le Quéré, C., Harrison, S., and Anderson, R.: Role of Marine Biology in Glacial-Interglacial CO2 Cycles, Science, 308, 74–78, 2005. a, b

Köhler, P. and Fischer, H.: Simulating changes in the terrestrial biosphere during the last glacial/interglacial transition, Global Planet. Change, 43, 33–55, 2004. a, b, c, d, e

Köhler, P., Fischer, H., and Schmitt, J.: Atmospheric δ13CO2 and its relation to pCO2 and deep ocean δ13C during the late Pleistocene, Paleoceanography, 25, PA1213,, 2010. a, b, c, d, e, f, g, h, i

Kohn, M.: Carbon isotope compositions of terrestrial C3 plants as indicators of (paleo)ecology and (paleo)climate, P. Natl. Acad. Sci., 107, 19691–19695, 2010. a, b, c, d

Kohn, M.: Carbon isotope discrimination in C3 land plants is independent of natural variations in pCO2, Geochem. Perspect. Lett., 2, 35–43, 2016. a, b, c, d

Kurahashi-Nakamura, T., Paul, A., and Losch, M.: Dynamical reconstruction of the global ocean state during the Last Glacial Maximum, Paleoceanography, 32, 326–350, 2017. a, b

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

Lambert, F., Tagliabue, A., Shaffer, G., Lamy, F., Winckler, G., Farias, L., Gallardo, L., and Pol-Holz, D.: Dust fluxes and iron fertilization in Holocene and Last Glacial Maximum climates, Geophys. Res. Lett., 42, 6014–6023, 2015. a, b, c

LeGrand, P. and Wunsch, C.: Constraints from paleotracer data on the North Atlantic circulation during the Last Glacial Maximum, Paleoceanogr. Paleoclim., 10, 1011–1045, 1995. a

Lindgren, A., Hugelius, G., and Kuhry, P.: Extensive loss of past permafrost carbon but a net accumulation into present-day soils, Lett. Nat., 560, 219–222, 2018. a

Lisiecki, L. E. and Raymo, M. E.: A Pliocene-Pleistocene stack of 57 globally distributed benthic δ18O records, Paleoceanography, 20, PA1003,, 2005. a, b, c

Liu, C., Kohl, A., Liu, Z., Wang, F., and Stammer, D.: Deep-reaching thermocline mixing in the equatorial pacific cold tongue, Nat. Commun., 7, 11576,, 2016. a

Lourey, M. J. and Trull, W.: Seasonal nutrient depletion and carbon export in the Subantarctic and Polar Frontal Zones of the Southern Ocean south of Australia, J. Geophys. Res.-Oceans, 106, 31463–31487, 2001. a

Lueker, T. J., Dickson, A. G., and Keeling, C. D.: Ocean pCO2 calculated from dissolved inorganic carbon, alkalinity, and equations for K-1 and K-2: validation based on laboratory measurements of CO2 in gas and seawater at equilibrium, Mar. Chem., 70, 105–119, 2000. a

MacFarling Meure, C., Etheridge, D., Trudinger, C., Steele, P., Langenfelds, R., van Ommen, T., Smith, A., and Elkins, J.: Law Dome CO2, CH4 and N2O ice core records extended to 2000 years BP, Geophys. Res. Lett., 33, L14810,, 2006. a, b

Marchitto, T., Lehman, S., Ortiz, J., Flückiger, J., and van Geen, A.: Marine Radiocarbon Evidence for the Mechanism of Deglacial Atmospheric CO2 Rise, Science, 316, 1456–1459, 2007. a, b, c, d, e

Marcott, S., Bauska, T., Buizert, C., Steig, E., Rosen, J., Cuffey, K., Fudge, T. J., Severinghaus, J. P., Ahn, J., Kalk, M. L., McConnell, J. R., Sowers, T., Taylor, K., White, J. W. C., and Brook, E.: Centennial-scale changes in the global carbon cycle during the last deglaciation, Nature, 514, 616–621, 2014. a, b

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

Martin, J. H., Knauer, G., Karl, D., and Broenkow, W.: VERTEX: carbon cycling in the northeast Pacific, Deep-Sea Res., 34, 267–285, 1987. a, b, c, d, e, f

Martinez-Garcia, A., Sigman, D., H.Ren, Anderson, R., Straub, M., Hodell, D., Jaccard, S., Eglinton, T., and Haug, G.: Iron Fertilization of the Subantarctic Ocean During the Last Ice Age, Science, 343, 1347–1350, 2014. a, b, c, d, e, f, g, h

Marzocchi, A. and Jansen, M. F.: Connecting Antarctic sea ice to deep-ocean circulation in modern and glacial climate simulations, Geophys. Res. Lett., 44, 6286–6295,, 2017. a, b

Matsumoto, K.: Biology-mediated temperature control on atmospheric pCO2 and ocean biogeochemistry, J. Geophys. Res., 34, L20605,, 2007. a, b

Mauritz, M., Celis, G., Ebert, C., Hutchings, J., Ledman, J., Natali, S., Pegoraro, E., Salmon, V., Schädel, C., Taylor, M., and Schuur, E.: Using stable carbon isotopes of seasonal ecosystem respiration to determine permafrost carbon loss, J. Geophys. Res.-Biogeosci., 124, 46–60, 2018. a, b

Menviel, L. and Joos, F.: Toward explaining the Holocene carbon dioxide and carbon isotope records: Results from transient ocean carbon cycle-climate simulations, Paleoceanography, 27, PA1207,, 2012. a, b, c

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, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x

Menviel, L., Timmermann, A., Friedrich, T., and England, M. H.: Hindcasting the continuum of Dansgaard–Oeschger variability: mechanisms, patterns and timing, Clim. Past, 10, 63–77,, 2014. a

Menviel, L., Mouchet, A., Meissner, K. J., Joos, F., and England, M. H.: Impact of oceanic circulation changes on atmospheric δ13CO2, Global Biogeochem. Cy., 29, 1944–1961, 2015. a

Menviel, L., Yu, J., Joos, F., Mouchet, A., Meissner, K. J., and England, M. H.: Poorly ventilated deep ocean at the Last Glacial Maximum inferred from carbon isotopes: A data-model comparison study, Paleoceanography, 31, 2–17, 2016. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s

Menviel, L., Spence, P., Skinner, L. C., Tachikawa, K., Friedrich, T., and Yu, L. M. J.: Enhanced Mid-depth Southward Transport in the Northeast Atlantic at the Last Glacial Maximum Despite a Weaker AMOC, Paleoceanogr. Paleoclim., 35, PA003793,, 2020. a, b, c

Millero, F. J.: Influence of pressure on chemical processes in the sea, in: Chemical Oceanography, 2nd edn, Vol. 8, Academic Press, New York, 1983. a

Monnin, E., Steig, E., Siegenthaler, U., Kawamura, K., Schwander, J., Stauffer, B., Stocker, T., Morse, D., Barnola, J., Bellier, B., Raynaud, D., and Fischer, H.: Evidence for substantial accumulation rate variability in Antarctica during the Holocene, through synchronization of CO2 in the Taylor Dome, Dome C and DML ice cores, Earth Planet. Sci. Lett., 224, 45–54, 2004. a, b

Morrison, A. and Hogg, A.: On the Relationship between Southern Ocean Overturning and ACC Transport, J. Phys. Oceanogr., 43, 140–148, 2013. a, b

Morrison, A., Hogg, A., and Ward, M.: Sensitivity of the Southern Ocean overturning circulation to surface buoyancy forcing, Geophys. Res. Lett., 38, L14602,, 2011. a

Morse, J. W. and Berner, R. A.: Dissolution kinetics of calcium carbonate in sea water. II: A kinetic origin for the lysocline, Am. J. Sci., 272, 840–851,, 1972. a

Muglia, J., Skinner, L., and Schmittner, A.: Weak overturning circulation and high Southern Ocean nutrient utilization maximized glacial ocean carbon, Earth Planet. Sci. Lett., 496, 47–56, 2018. a, b, c, d, e, f, g, h, i, j, k, l

Muscheler, R., Beer, J., Wagner, G., Laj, C., Kissel, C., Raisbeck, G., Yioud, F., and Kubik, P.: Changes in the carbon cycle during the last deglaciation as indicated by the comparison of 10Be and 14C records, Earth Planet. Sci. Lett., 219, 325–340, 2014. a, b, c, d, e

Oliver, K. I. C., Hoogakker, B. A. A., Crowhurst, S., Henderson, G. M., Rickaby, R. E. M., Edwards, N. R., and Elderfield, H.: A synthesis of marine sediment core δ13C data over the last 150 000 years, Clim. Past, 6, 645–673,, 2010. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o

O'Neill, C. M., Hogg, A. McC., Ellwood, M. J., Eggins, S. M., and Opdyke, B. N.: The [simple carbon project] model v1.0, Geosci. Model Dev., 12, 1541–1572,, 2019. a, b, c, d, e, f, g, h, i, j, k

O'Neill, C. M., Hogg, A. McC., Ellwood, M. J., Opdyke, B. N., and Eggins, S. M.: [simple carbon project] model v2.2 (Version V2.2), Zenodo,, 2021. a

Opdyke, B. and Walker, J.: Return of the coral reef hypothesis: Basin to shelf partitioning of CaCO3 and its effect on atmospheric CO2, Geology, 20, 733–736, 1992. a, b, c, d, e, f, g

Otto, D., Rasse, D., Kaplan, J., Warnant, P., and Francois, L.: Biospheric carbon stocks reconstructed at the Last Glacial Maximum: comparison between general circulation models using prescribed and computed sea surface temperatures, Global Planet. Change, 33, 117–138, 2002. 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

Piotrowski, A., Banakar, V., Scrivner, A., Elderfield, H., Galy, A., and Dennis, A.: Indian Ocean circulation and productivity during the last glacial cycle, Earth Planet. Sci. Lett., 285, 179–189, 2009. a, b, c, d

Prentice, I. C., Harrison, S., and Bartlein, P.: Global vegetation and terrestrial carbon cycle changes after the last ice age, New Phytol., 189, 988–998, 2011. a, b, c

Qin, B., Li, T., Xiong, Z., Algeo, T. J., and Chang, F.: Deepwater carbonate ion concentrations in the western tropical Pacific since 250 ka: Evidence for oceanic carbon storage and global climate influence, Paleoceanography, 32, 351–370, 2017. a, b, c

Qin, B., Li, T., Xiong, Z., Algeo, T., and Jia, Q.: Deep-Water Carbonate Ion Concentrations in the Western Tropical Pacific Since the Mid-Pleistocene: A Major Perturbation During the Mid-Brunhes, J. Geophys. Res.-Oceans, 123, 6876–6892,, 2018. a, b, c

Reimer, P., Baillie, M., Bard, E., Bayliss, A., Beck, J., Blackwell, P., Ramsey, C. B., Buck, C., Burr, G., Edwards, R., Friedrich, M., Grootes, P., Guilderson, T., Hajdas, I., Heaton, T., Hogg, A., Hughen, K., Kaiser, K., Kromer, B., McCormac, F., Manning, S., Reimer, R., Richards, D., Southon, J., Talamo, S., Turney, C., van der Plicht, J., and Weyhenmeyer, C.: IntCal09 and Marine09 radiocarbon age calibration curves, 0-50,000 years cal BP., Radiocarbon, 51, 1111–50, 2009. a, b, c

Ridgwell, A.: An end to the ”rain ratio” reign?, Geochem. Geophys. Geosyst., 4, 1051,, 2003. a

Ridgwell, A., Watson, A., Maslin, M., and Kaplan, J.: Implications of coral reef buildup for the controls on atmospheric CO2 since the Last Glacial Maximum, Paleoceanography, 18, 1083,, 2003. a, b, c, d, e

Ridgwell, A. J.: Glacial-interglacial perturbations in the global carbon cycle, Ph.D. thesis, University of East Anglia, Norwich, UK, 2001. a

Rohling, E., Grant, K., Bolshaw, M., Roberts, A., Siddall, M., Hemleben, C., and Kucera, M.: Antarctic temperature and global sea level closely coupled over the past five glacial cycles, Nat. Geosci., 2, 500–504, 2009. a, b, c, d, e

Ronge, T., Tiedemann, R., Lamy, F., Kohler, P., Alloway, B., Pol-Holz, R. D., Pahnke, K., Southon, J., and Wacker, L.: Radiocarbon constraints on the extent and evolution of the South Pacific glacial carbon pool, Nat. Commun., 7, 11487,, 2016. a, b

Rubino, M., Etheridge, D., Trudinger, C., Allison, C., Battle, M., Langenfelds, R., Steele, L., Curran, M., Bender, M., White, J., Jenk, T., Blunier, T., and Francey, R.: A revised 1000 year atmospheric δ13C-CO2 record from Law Dome and South Pole, Antarctica, J. Geophys. Res., 118, 8482–8499, 2013. a, b

Sarmiento, J. L. and Gruber, N.: Ocean biogeochemical dynamics, Princeton University Press, Princeton, New Jersey, 2006. a, b, c, d, e

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

Schmitt, J., Schneider, R., Elsig, J., Leuenberger, D., Lourantou, A., Chappellaz, J., Köhler, P., Joos, F., Stocker, T., Leuenberger, M., and Fischer, H.: Carbon Isotope Constraints on the Deglacial CO2 Rise from Ice Cores, Science, 336, 711–714, 2012. a, b, c, d

Schneider, R., Schmitt, J., Köhler, P., Joos, F., and Fischer, H.: A reconstruction of atmospheric carbon dioxide and its stable carbon isotopic composition from the penultimate glacial maximum to the last glacial inception, Clim. Past, 9, 2507–2523,, 2013. a, b, c, d, e, f, g, h, i

Schubert, B. and Jahren, A.: The effect of atmospheric CO2 concentration on carbon isotope fractionation in C3 land plants, Geochim. Cosmochim. Acta, 96, 29–43, 2012. a, b, c

Shaffer, G. and Lambert, F.: In and out of glacial extremes by way of dust-climate feedbacks, P. Natl. Acad. Sci., 115, 2026–2031, 2018. a, b

Siani, E., Michel, E., Pol-Holz, R. D., DeVries, T., Lamy, F., Carel, M., Isguder, G., Dewilde, F., and Lourantou, A.: Carbon isotope records reveal precise timing of enhanced Southern Ocean upwelling during the last deglaciation, Nat. Commun., 4, 1–9,, 2013. a, b, c, d

Siegel, D. A., Buesseler, K. O., Doney, S. C., Sailley, S. F., Behrenfeld, M. J., and Boyd, P. W.: Global assessment of ocean carbon export by combining satellite observations and foodweb models, Global Biogeochem. Cy., 28, 181–196, 2014. a, b, c

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

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

Sigman, D., Hain, M., and Haug, G.: The polar ocean and glacial cycles in atmospheric CO2 concentration, Nat. Rev., 466, 47–55, 2010. a, b, c

Sikes, E., Samson, C., Guilderson, T., and Howard, W.: Old radiocarbon ages in the southwest Pacific Ocean during the last glacial period and deglaciation, Nature, 405, 555–559, 2000. a, b, c, d

Sikes, E., Cook, M., and Guilderson, T.: Reduced deep ocean ventilation in the Southern Pacific Ocean during the last glaciation persisted into the deglaciation, Earth Planet. Sci. Lett., 438, 130–138, 2016. a, b

Skinner, L. and Shackleton, N. J.: Rapid transient changes in northeast Atlantic deep water ventilation age across Termination I, Paleoceanography, 19, PA2005,, 2004. a

Skinner, L., Waelbroeck, C., Scrivner, A., and Fallon, S.: Radiocarbon evidence for alternating northern and southern sources of ventilation of the deep Atlantic carbon pool during the last deglaciation, P. Natl. Acad. Sci., 111, 5480–5484, 2014. a

Skinner, L., McCave, I., Carter, L., Fallon, S., Scrivner, A., and Primeau, F.: Reduced ventilation and enhanced magnitude of the deep Pacific carbon pool during the last glacial period, Earth Planet. Sci. Lett., 411, 45–52, 2015. a, b, c

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

Skinner, L. C., 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, b, c, d

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

Stott, L., Southon, J., Timmermann, A., and Koutavas, A.: Radiocarbon age anomaly at intermediate water depth in the Pacific Ocean during the last deglaciation, Paleoceanography, 24, PA2223,, 2009. a

Strutz, T.: Data Fitting and Uncertainty. A practical introduction to weighted least squares and beyond., vol. 2 of Wiesbaden, Springer Vieweg, 2016. a

Tagliabue, A., Bopp, L., Roche, D. M., Bouttes, N., Dutay, J.-C., Alkama, R., Kageyama, M., Michel, E., and Paillard, D.: Quantifying the roles of ocean circulation and biogeochemistry in governing ocean carbon-13 and atmospheric carbon dioxide at the last glacial maximum, Clim. Past, 5, 695–706,, 2009. a

Takahashi, T., Sutherland, S., Feely, R., and Cosca, C.: Decadal variation of the surface water PCO2 in the western and central equatorial Pacific, Science, 302, 852–856, 2003. a

Talley, L.: Closure of the global overturning circulation through the Indian, Pacific, and Southern Oceans: Schematics and transports, Oceanography, 78, 257–303, 2013. a, b, c, d, e, f, g, h

Tamburini, F. and Föllmi, K. B.: Phosphorus burial in the ocean over glacial-interglacial time scales, Biogeosciences, 6, 501–513,, 2009. a

Tarnocai, C., Canadell, J., Schuur, E., Kuhry, P., Mazhitova, G., and Zimov, S.: Soil organic carbon pools in the northern circumpolar permafrost region, Global Biogeochem. Cy., 23, GB2023,, 2009. a, b, c, d, e

Toggweiler, J. and Sarmiento, J.: Glacial to interglacial changes in atmospheric carbon dioxide: The critical role of ocean surface water in high latitudes, in: The Carbon Cycle and Atmospheric CO2: Natural Variations Archean to Present, Geophysical Monograph Series, American Geophysical Union, 32, 163–184, 1985. a

Toggweiler, J. R.: Variation of atmospheric CO2 by ventilation of the ocean's deepest water, Paleoceanography, 14, 571–588, 1999. a

Toggweiler, J. R.: Origin of the 100,000-year time scale in Antarctic temperatures and atmospheric CO2, Paleoceanography, 23, PA2211,, 2008. a, b, c

Treat, C., Kleinen, T., Broothaerts, N., Dalton, A., Dommain, R., Douglas, T., Drexler, J., Finkelstein, S., Grosse, G., Hope, G., Hutchings, J., Jones, M., Kuhry, P., Lacourse, T., Lähteenoja, O., Loisel, J., Notebaert, B., Payne, R., Peteet, D., Sannel, A., Stelling, J., Strauss, J., Swindles, G., Talbot, J., Tarnocai, C., Verstraeten, G., C.J. Williams, Z. X., Yu, Z., Väliranta, M., Hättestrand, M., Alexanderson, H., and Brovkin, V.: Widespread global peatland establishment and persistence over the last 130,000 y, P. Natl. Acad. Sci., 116, 4822–4827, 2019. a, b, c, d, e

van der Sleen, P. et al.: No growth stimulation of tropical trees by 150 years of CO2 fertilization but water-use efficiency increased, Nat. Geosci., 8, 24–28, 2015. a, b, c

Vecsei, A. and Berger, W.: Increase of atmospheric CO2 during deglaciation: Constraints on the coral reef hypothesis from patterns of deposition, Global Biogeochem. Cy., 18, GB1035,, 2004. a

Walker, J. and Kasting, J.: Effects of fuel and forest conservation on future levels of atmospheric carbon dioxide, Palaeogeogr. Palaeoclim. Palaeoecol., 97, 151–189, 1992. a, b

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

Weeding, B. and Trull, W.: Hourly oxygen and total gas tension measurements at the Southern Ocean time series site reveal winter ventilation and spring net community production, J. Geophys. Res.-Oceans, 119, 348–358, 2004. a

Winckler, G., Anderson, R., Jaccard, S., and Marcantonio, F.: Ocean dynamics, not dust, have controlled equatorial Pacific productivity over the past 500,000 years, P. Natl. Acad. Sci., 113, 6119–6124, 2016. a

Wolff, E., Barbante, C., Becagli, S., Bigler, M., Boutron, C., Castellano, E., de Angelis, M., Federer, U., Fischer, H., Fundel, F., Hansson, M., Hutterli, M., Jonsell, U., Karlin, T., Kaufmann, P., Lambert, F., Littot, G., Mulvaney, R., Roöthlisberger, R., Ruth, U., Severi, M., Siggaard-Andersen, M., Sime, L., Steffensen, J., Stocker, T., Traversi, R., Twarloh, B., Udisti, R., Wagenbach, D., and Wegner, A.: Changes in environment over the last 800,000 years from chemical analysis of the EPICA Dome C ice core, Quaternary Sci. Rev., 29, 285–95, 2010.  a, b, c, d, e, f, g

Yamamoto, A., Abe-Ouchi, A., Ohgaito, R., Ito, A., and Oka, A.: Glacial CO2 decrease and deep-water deoxygenation by iron fertilization from glaciogenic dust, Clim. Past, 15, 981–996,, 2019. a

Yu, J., Broecker, W., Elderfield, H., Jin, Z., McManus, J., and Zhang, F.: Loss of Carbon from the Deep Sea Since the Last Glacial Maximum, Science, 330, 1084–1087, 2010. a

Yu, J., Anderson, R., Jin, Z., Rae, J., Opdyke, B., and Eggins, S.: Responses of the deep ocean carbonate system to carbon reorganization during the Last Glacial-interglacial cycle, Quaternary Sci. Rev., 76, 39–52, 2013. a, b, c, d, e, f

Yu, J., Anderson, R., Jin, Z., Menviel, L., Zhang, F., Ryerson, F., and Rohling, E.: Deep South Atlantic carbonate chemistry and increased interocean deep water exchange during last deglaciation, Quaternary Sci. Rev., 90, 80–89, 2014a. a, b, c, d, e, f

Yu, J., Anderson, R. F., and Rohling, E. J.: Deep ocean carbonate chemistry and glacial-interglacial atmospheric CO2 changes, Oceanography, 27, 16–25, 2014b. a, b, c, d, e, f, g, h

Yu, J., Menviel, L., Jin, Z. D., Thornalley, D., Barker, S., Marino, G., Rohling, E. J., Cai, Y., Zhang, F., Wang, X., Dai, Y., Chen, P., and Broecker, W. S.: Sequestration of carbon in the deep Atlantic during the last glaciation, Nat. Geosci., 9, 319–325, 2016. a, b, c, d, e, f, g

Yu, J., Menviel, L., Jin, Z., Thornalley, D., Foster, G., Rohling, E., McCave, I., McManus, J., Dai, Y., Ren, H., He, F., Zhang, F., Chen, P., and Roberts, A.: More efficient North Atlantic carbon pump during the Last Glacial Maximum, Nat. Commun., 10, 1–11,, 2019. a

Zeebe, R. E.: LOSCAR: Long-term Ocean-atmosphere-Sediment CArbon cycle Reservoir Model v2.0.4, Geosci. Model Dev., 5, 149–166,, 2012. a, b

Zhao, N., Marchal, O., Keigwin, L., Amrhein, D., and Gebbie, G.: A Synthesis of Deglacial Deep-Sea Radiocarbon Records and Their (In)Consistency With Modern Ocean Ventilation, Paleoceanogr. Paleoclim., 33, 128–151, 2017. a, b, c

Short summary
We undertake a model–data study of the last glacial–interglacial cycle of atmospheric CO2, spanning 0–130 ka. We apply a carbon cycle box model, constrained with glacial–interglacial observations, and solve for optimal model parameter values against atmospheric and ocean proxy data. The results indicate that the last glacial drawdown in atmospheric CO2 was delivered mainly by slowing ocean circulation, lower sea surface temperatures and also increased Southern Ocean biological productivity.