Articles | Volume 15, issue 3
Research article
18 Jun 2019
Research article |  | 18 Jun 2019

Coupled climate–carbon cycle simulation of the Last Glacial Maximum atmospheric CO2 decrease using a large ensemble of modern plausible parameter sets

Krista M. S. Kemppinen, Philip B. Holden, Neil R. Edwards, Andy Ridgwell, and Andrew D. Friend

During the Last Glacial Maximum (LGM), atmospheric CO2 was around 90 ppmv lower than during the pre-industrial period. The reasons for this decrease are most often elucidated through factorial experiments testing the impact of individual mechanisms. Due to uncertainty in our understanding of the real system, however, the different models used to conduct the experiments inevitably take on different parameter values and different structures. In this paper, the objective is therefore to take an uncertainty-based approach to investigating the LGM CO2 drop by simulating it with a large ensemble of parameter sets, designed to allow for a wide range of large-scale feedback response strengths. Our aim is not to definitely explain the causes of the CO2 drop but rather explore the range of possible responses. We find that the LGM CO2 decrease tends to predominantly be associated with decreasing sea surface temperatures (SSTs), increasing sea ice area, a weakening of the Atlantic Meridional Overturning Circulation (AMOC), a strengthening of the Antarctic Bottom Water (AABW) cell in the Atlantic Ocean, a decreasing ocean biological productivity, an increasing CaCO3 weathering flux and an increasing deep-sea CaCO3 burial flux. The majority of our simulations also predict an increase in terrestrial carbon, coupled with a decrease in ocean and increase in lithospheric carbon. We attribute the increase in terrestrial carbon to a slower soil respiration rate, as well as the preservation rather than destruction of carbon by the LGM ice sheets. An initial comparison of these dominant changes with observations and paleoproxies other than carbon isotope and oxygen data (not evaluated directly in this study) suggests broad agreement. However, we advise more detailed comparisons in the future, and also note that, conceptually at least, our results can only be reconciled with carbon isotope and oxygen data if additional processes not included in our model are brought into play.

1 Introduction

Analyses of Antarctic ice core records suggest that the atmospheric CO2 concentration at the Last Glacial Maximum (LGM), about 21 kyr ago, was around 190 ppmv, well below the pre-industrial atmospheric concentration of around 280 ppmv. The most commonly accepted mechanisms to explain the atmospheric CO2 decrease include lower sea surface temperatures, which increase the ocean CO2 solubility (Martin et al., 2005; Menviel et al., 2012), enhanced ocean biological pump due to increased input of aeolian iron to the ocean (Bopp et al., 2003; Oka et al., 2011; Jaccard et al., 2013; Ziegler et al., 2013; Martínez-García et al., 2014; Lambert et al., 2015), capping of air–sea gas exchange by expanding sea ice (Stephens and Keeling, 2000; Sun and Matsumoto, 2010; Chikamoto et al., 2012) and ocean circulation/stratification changes increasing the net CO2 flux into the ocean (Adkins et al., 2002; Lynch-Stieglitz et al., 2007; Skinner et al., 2010, 2014; Lippold et al., 2012; Gebbie, 2014; Tiedemann et al., 2015; de la Fuente et al., 2015; Freeman et al., 2015). These changes may in turn be due to a range of possible mechanisms such as increased brine rejection (Shin et al., 2003; Bouttes et al., 2010, 2011; Zhang et al., 2013; Ballarotta et al., 2014), a shift in/weakening of the westerly wind belt over the Southern Ocean (Toggweiler et al., 2006; Anderson et al., 2009; Völker and Köhler, 2013) and a reduced or reversed buoyancy flux from the atmosphere to the ocean surface in the Southern Ocean (Watson and Naveira Garabato, 2006; Ferrari et al., 2014). A process that is conversely assumed to have contributed to increasing atmospheric CO2 is increasing salinity and ocean total dissolved inorganic carbon (DIC) concentration in response to decreasing sea level (Ciais et al., 2013).

A dominant assumption is also that the terrestrial biosphere carbon inventory was reduced (Crowley, 1995; Adams and Faure, 1998; Ciais et al., 2012; Peterson et al., 2014), in line with independent estimates of an ocean carbon inventory that was enhanced by several hundred petagrams (Goodwin and Lauderdale, 2013; Sarnthein et al., 2013; Allen et al., 2015; Skinner et al., 2015; Schmittner and Somes, 2016). The decrease in terrestrial carbon is generally attributed to unfavourable climatic conditions for photosynthesis and the destruction of organic material by moving ice sheets (e.g. Otto et al., 2002; Prentice et al., 2011; Brovkin et al., 2012; O'ishi and Abe-Ouchi, 2013). The hypothesis that there was an increase in terrestrial carbon has, however, also been put forward (e.g. Zeng, 2003; Zimov et al., 2006), with some studies additionally suggesting little net change (e.g. Brovkin and Ganopolski, 2015). Processes proposed to be responsible for the terrestrial carbon increase include growth in “inert” or permafrost carbon, slower “active” soil respiration rates, continental shelf regrowth and the preservation rather than destruction of terrestrial biosphere carbon in areas to be covered by the expanding Laurentide and Eurasian ice sheets (Weitemeyer and Buffett, 2006; Franzén and Cropp, 2007; Zeng, 2007; Zimov et al., 2009; Zech et al., 2011).

Other mechanisms which may have affected the LGM atmospheric CO2 change include changes in carbonate weathering rate, through its control on the ocean ALK:DIC ratio and consequently the solubility of CO2 (Munhoven, 2002; Jones et al., 2002; Foster and Vance, 2006; Vance et al., 2009; Brovkin et al., 2012; Crocket et al., 2012; Lupker et al., 2013; Simmons et al., 2016). The change in carbonate weathering rate would in turn have been caused by lower sea level, exposing previously submerged rock, the presence of a greater amount of glacial flour, which is more susceptible to weathering (Kohfeld and Ridgwell, 2009), or potentially higher soil carbon content. The lower sea level may also have reduced shallow water carbonate deposition by decreasing the area of shallow ocean (Opdyke and Walker, 1992; Kleypas, 1997; Brovkin et al., 2007), and increased oceanic PO4 inventory, alleviating the PO4 limitation on marine production (Tamburini and Föllmi, 2009; Wallmann, 2014, 2015). Other potential CO2 mechanisms include decreasing dissolved organic carbon inventory due to a more stratified deep ocean (Ma and Tian, 2014) and reduced marine bacterial metabolic rate in response to lower ocean temperatures. The lower metabolic rate acts to decrease the return rate of DIC from the remineralisation of organic material and hence the concentration of CO2 at the ocean surface (Matsumoto et al., 2007; Roth et al., 2014). The net flux of CO2 into the ocean may also have increased due to enhanced diatom production caused by the leakage of silicic acid trapped in the Southern Ocean (Matsumoto et al., 2002, 2014) or increased Si inventory, caused by increased input of Si from wind-born dust or enhanced weathering (Harrison, 2000; Tréguer and Pondaven, 2000).

Mechanisms put forward to explain the LGM atmospheric CO2 decrease arise from paleodata and model studies. The latter most often involve factorial experiments, introducing mechanisms one at a time. There is rarely any investigation of the impact of alternative assumptions regarding parameter values or model structure. An example of a relevant study is Bouttes et al. (2011), which varied model parameters controlling the importance of iron fertilisation, brine rejection and stratification-dependent diffusion in an ensemble setting, assessing the agreement of the model output with data. Here, our aim is conversely to take an uncertainty-based approach to investigating the LGM CO2 drop by simulating it with a large ensemble of parameter sets designed to allow for a wide range of large-scale feedback response strengths (Holden et al., 2013a). The objective is not to definitely explain the causes of the CO2 drop but rather explore the range of possible responses. By “responses” we mean physical and biogeochemical changes in the Earth system (e.g. change in global particulate organic carbon export flux) and how these might be linked to ΔCO2 and to each other, rather than specific mechanisms (e.g. iron fertilisation). Knowledge of these relationships can in turn inform analysis, in the future, of the relationship between the ensemble parameters and model outputs, in order to isolate individual LGM CO2 mechanisms. In this study, we furthermore seek to simulate the LGM atmospheric CO2 drop with the simulated CO2 feeding back to the simulated climate, which is still infrequently done in LGM CO2 experiments, and the first time it is done with GENIE-1. Moreover, rather than assuming that terrestrial carbon is destroyed by the LGM ice sheets, we assume that it is gradually buried. This assumption has not yet been implemented, in GENIE-1 or other models, in an equilibrium set-up.

Despite our ensemble varying many of the parameters thought to contribute to variability in glacial–interglacial atmospheric CO2, not all sources of uncertainty can be captured, and this is reflected in our simulated ΔCO2 distribution. We estimate that up to ∼60 ppmv of ΔCO2 could be due to processes not included in our model and error in our process representations (see Sect. 2.4 for details). We thus treat ΔCO2 between -90 and −30 ppmv as “equally plausible” and focus on describing the physical and biogeochemical changes seen in the subset of simulations with this ΔCO2. We also conduct an initial assessment of how the subset mean and/or dominant (in terms of sign) responses compare against observations and paleoproxies, including temperature, sea ice, precipitation, AMOC and AABW cell strengths, terrestrial carbon, ocean carbon, particulate organic matter export and deep-sea CaCO3 burial.

Finally, to test the robustness of relationships derived from the analysis of the ensemble subset with ΔCO2 between -90 and −30 ppmv, we briefly compare the physical and biogeochemical changes seen therein with the changes seen in the ensemble with no ΔCO2 filter and the ensemble with a more negative ΔCO2 filter (-90 to −60 ppmv) (Sect. 2.4). In general, the same dominant relationships between ΔCO2 and the physical and biogeochemical changes are observed as in the subset with ΔCO2 between -90 and −30 ppmv. In the case of the ensemble subset with ΔCO2 between -90 and −60 ppmv, we additionally look at what proportion of the total terrestrial carbon change comes from within the ice sheet areas and from there draw conclusions for the rest of the ensemble.

The paper is organised as follows. Section 2 describes the model, the ensemble, the simulation set-up and the ensemble subsets to be analysed. Section 3 is the results and discussion section, which includes a brief evaluation of the pre-industrial (control) spin-up simulation to verify reproducibility of Holden et al. (2013a). The majority of the section is devoted to the LGM simulation: namely, diagnosis of the physical and biogeochemical changes (including potential causal relationships) seen in the subset with ΔCO2 between -30 and −90 ppmv, and to a lesser extent, the ensemble with both more and less constrained ΔCO2. Comparison of the first subset against observations and paleoproxies is also included. Section 4 provides the key conclusions.

2 Methods

2.1 The model

The GENIE-1 configuration is as described in Holden et al. (2013a). The physical model consists of a three-dimensional frictional geostrophic ocean model (GOLDSTEIN) coupled to a thermodynamic/dynamic sea ice model (Edwards and Marsh, 2005; Marsh et al., 2011) and a two-dimensional Energy–Moisture Balance Model (EMBM). Atmospheric tracers are a subcomponent of the EMBM, with a simple module (ATCHEM) used to store the concentration of atmospheric gases and their relevant isotopic properties (Lenton et al., 2007). The model land surface physics and terrestrial carbon cycle are represented by an efficient numerical terrestrial scheme (ENTS) (Williamson et al., 2006). The ocean biogeochemistry model (BIOGEM) is as described in Ridgwell et al. (2007) but includes a representation of iron cycling (Annan and Hargreaves, 2010) and the biological uptake scheme of Doney et al. (2006). The model sediments are represented by SEDGEM (Ridgwell and Hargreaves, 2007). GENIE-1 also includes a land surface weathering model, ROKGEM (Colbourn, 2011), which redistributes prescribed weathering fluxes according to a fixed river-routing scheme. The model is on a 36×36 equal-area horizontal grid, with 16 vertical levels in the ocean.

2.2 The simulation ensemble

The GENIE-1 ensemble consists of 471 parameter sets, varying 29 key model parameters over the ranges in Table 1. It derives from the 471-member emulator-filtered plausibility constrained ensemble of Holden et al. (2013a), which varies 24 active parameters and 1 dummy parameter (as a check against overfitting). The parameter values in Holden et al. (2013a) were derived by building emulators of eight pre-industrial climate metrics and applying a rejection sampling method known as approximate Bayesian computation (ABC) to find parameter sets that the emulators predicted were modern plausible. Two parameters were later added to the ensemble, in Holden et al. (2013b), to describe the unmodelled response of clouds to global average temperature change (OL1) (see Appendix A for further information) and the uncertain response of photosynthesis to changing atmospheric CO2 concentration (VPC). The parameters are as described in Holden et al. (2013b). We add two further parameters here that represent uncertain processes specific to the LGM. The first (FFX) scales ice sheet meltwater fluxes to account for uncertainty in unmodelled isostatic depression at the ice–bedrock interface due to ice sheet growth and for assuming a fixed land–sea mask (Holden et al., 2010b). We vary the parameter in the ensemble to capture the uncertainty in the magnitude of the glacial sea level drop and its effects on the carbon cycle. The second (GWS) scales the global average pre-industrial carbonate weathering rates for the LGM, to account for uncertainty in carbonate weathering and unmodelled shallow water carbonate deposition rate changes. For both FFX and GWS, uniform random values were derived using the generation function “runif” in R.

Table 1Ensemble parameters. Ranges are from (a) Holden et al. (2013a), (b) Holden et al. (2013b) and (c) Holden et al. (2010b), with the exception of GWS (see main text). The table also precludes the dummy parameter.

n/a – not applicable

Download Print Version | Download XLSX

2.3 Experimental set-up of the model

The pre-industrial ensemble simulation results were repeated to verify reproducibility of Holden et al. (2013a). The simulations were performed in two stages, each lasting 10 kyr, on the Cambridge high-performance computing (HPC) cluster Darwin. The first stage involved spinning up the model with atmospheric CO2 concentration relaxed to 278 ppmv and a closed biogeochemistry system. This means that there are no sediment–ocean interactions and the model forces the CaCO3 weathering and deep-sea sediment burial rates into balance. An initial CaCO3 weathering flux is prescribed, but this is subsequently rescaled internally to balance the modelled CaCO3 burial rate and conserve alkalinity. In the second stage, atmospheric CO2 was allowed to evolve freely, with interacting oceans and sediments, and the CaCO3 weathering rate is set equal to the CaCO3 burial rate diagnosed from the end of stage 1. To allow the sediments to reach equilibrium as fast as possible, no bioturbation was modelled in either stage 1 or stage 2.

Each parameter set was then applied to LGM simulations. The modelled pre-industrial equilibrium states were used as initial conditions and the ensemble members were integrated for 10 kyr, with freely evolving CO2. These 10 kyr simulations are variously referred to here as the “LGM equilibrium simulation” or “stage 3”, and the LGM equilibrium state refers to the end of stage 3 (see Sect. S1 in the Supplement for more details). After application to stages 2 and 3, the original 471 ensemble members were filtered to 315 ensemble members to exclude those simulations with a stage 2 atmospheric CO2 concentration outside of the range 268 to 288 ppmv (see Prentice et al., 2001), those that entered a snowball Earth state in stage 3 (global annual SAT between -68 and −57C) or those that showed evidence of numerical instability (see Holden et al., 2013b).

Boundary conditions applied in the LGM simulations included orbital parameters (Berger, 1978) and aeolian dust deposition fields (Mahowald et al., 2006). The atmospheric CO2 used in the radiative code is internally generated, rather than prescribed, but the radiative forcing from dust and gases other than CO2 was neglected. The model also requires a detrital flux field to the sediments, containing contributions from opal and material from non-aeolian sources (Ridgwell and Hargreaves, 2007). Weathering fluxes from the pre-industrial simulation were applied, scaled by GWS (the land-to-ocean bicarbonate flux scaling factor).

The representation of the ice sheets is as described in Holden et al. (2010b), using the terrestrial ice sheet fraction and orography from the ICE-4G reconstruction of Peltier (1994). Rather than initialising the ensemble with the ice sheet extent and orography at 21 kyr BP, the ice sheets are configured to grow from their pre-industrial to LGM extent in 1 kyr, at the beginning of the LGM simulation (i.e. 0–1 kyr) in order to account for the impact of sea level change on ocean tracers. Following Holden et al. (2010b), only the Laurentide and Eurasian ice sheets are allowed to change from their pre-industrial form (accounting for ∼80 % of global ice sheet change), and we also route the freshwater to build the ice sheets from the Atlantic, Pacific and Arctic, assuming modern topography, rather than extracting it uniformly. As the ice sheets grow, grid cells on land are gradually covered by ice, and any carbon that remains, or has accumulated, is preserved underneath (“buried”). Once the ice covers a grid cell, there is no more exchange between the land carbon in that grid cell and the atmosphere. Prior to being buried, however, it is subject to the same forcings as carbon in any other grid cell. To determine how sensitive the burial carbon amount (i.e. the amount of carbon that is available for preservation underneath the ice) is to the duration of ice sheet build-up, we test the impact of varying the latter from 1 to 10 kyr for one ensemble member (extending the total simulation length to 11 kyr). Our assumption is that if the difference is negligible, applying the same ensemble member to a transient simulation of the full glacial cycle (and therefore a more realistic ice sheet build-up history) would not have yielded a dramatically different burial carbon inventory. We find that increasing the ice sheet build-up duration indeed changes the burial carbon amount only marginally: an increase of ∼34 PgC. A limitation, however, is that we do not have a way of testing if the response of other ensemble members would be equally subdued.

2.4 Ensemble subsets

Although our ensemble varies many of the parameters thought to contribute to variability in glacial–interglacial atmospheric CO2, not all sources of uncertainty can be captured. We estimate, based on our expert opinion, that up to ∼60 ppmv of ΔCO2 could be due to error in our process representations and processes not included in our model, such as changing marine bacterial metabolic rate, wind speed (via its effect on gas transfer) and Si fertilisation. This is not a comprehensive assessment, however, as our model also does not include processes such as the effect of changing winds on ocean circulation (Toggweiler et al., 2006), Si leakage (Matsumoto et al., 2002, 2013, 2014), the effect of decreasing sea surface temperatures (SSTs) on CaCO3 production (Iglesias-Rodriguez et al., 2002) or changing oceanic PO4 inventory (Menviel et al., 2012). We focus our analyses on the subset of the ensemble with ΔCO2 between -90 and −30 ppmv (Table 2), treating each value in this range as equally plausible. To test the robustness of diagnosed relationships, we also briefly compare the response of this subset (ENS104) with the response of the ensemble with no ΔCO2 filter (ENS315) and the response of the ensemble with a more negative ΔCO2 filter (ENS16). In ENS16, the upper ΔCO2 limit is set to -60 ppmv, roughly equivalent to allowing for an extra atmospheric CO2 decrease due to changing marine bacterial metabolic rate, wind speed (via its effect on gas transfer) and Si fertilisation, between the best and upper estimate of Kohfeld and Ridgwell (2009). The ΔCO2 distribution in each subset or ensemble is shown in Fig. 1.

Table 2Ensemble subsets, including ΔCO2 and number of members in each.

Download Print Version | Download XLSX

Figure 1LGM change in atmospheric CO2 distribution. The ENS315 response is shown in grey, the ENS104 ensemble response in yellow and the ENS16 ensemble response in orange. The same colour legend applies to all figures in the paper.


3 Results and discussion

3.1 Pre-industrial simulations

Comparison of the pre-industrial response of ENS315 (i.e. the original, non-ΔCO2 filtered ensemble) against the pre-industrial ensemble response of Holden et al. (2013a) confirms that the two are very similar. We additionally evaluate ENS315 against a few additional pre-industrial metrics (see Sect. S2) and find responses that can be deemed plausible, following the design principles for the ensemble, outlined in Holden et al. (2013a).

Figure 2LGM change in surface air temperature and sea surface temperature (a–b) distributions.


3.2 LGM simulations

3.2.1 Climate, sea level and ocean circulation


The ENS104 mean LGM surface air temperature (SAT) anomaly (ΔSAT) is -4.6±1.7, and the range is −2.5 to −10.4C. The mean is close to the observed ΔSAT of -4±1C (Annan and Hargreaves, 2013) and the range roughly equivalent to the range of previous model-based estimates (Kim et al., 2003; Masson-Delmotte et al., 2006; Schneider von Deimling et al., 2006; Braconnot et al., 2007; Holden et al., 2010a; Brady et al., 2013). The ENS104 mean LGM SST anomaly (ΔSST) is -1.8±0.8C, and the range is −4.5 and −0.7C. The mean is again close to an observational data-constrained model estimate (Schmittner et al., 2011) and within the range of estimates inferred from proxy data (MARGO Project Members 2009 in Masson-Delmotte et al., 2013). There is a positive correlation between ΔSAT and ΔCO2 (r=0.75, 0.05 significance level henceforth), most likely reflecting the radiative impact of atmospheric CO2 on SAT, as well as the effect of changing SAT on ΔCO2. As suggested above, decreasing SST may contribute to decreasing CO2 via the CO2 solubility temperature dependence. Changing SAT may also affect ΔCO2 via its effects on sea ice, ocean circulation, terrestrial and marine productivity (see below). The positive correlation is reproduced in ENS315 (r=0.74), and as shown in Fig. 2, ΔSAT and ΔSST tend to be less negative in ENS315 than in ENS104. In ENS16, ΔSAT and ΔSST are from the extreme or at least lower end of the ENS104 range.

Figure 3LGM change in surface air temperature (a–b) and sea surface temperature (c–d) (C) ENS104 mean (a, c) and standard deviation (b, d).


The ENS104 mean ΔSAT and ΔSST spatial distributions are shown in Fig. 3. In line with observations (Annan and Hargreaves, 2013), the largest SAT decreases (>10C) are simulated over the Laurentide and Eurasian ice sheets. The Equator-to-pole temperature gradient is also broadly reproduced. The largest SST decreases (≥4C) are found in the North Atlantic and northeast Pacific, with more limited cooling (-2C) in the tropics and polar regions, again consistent with observations. However, the largest SST decreases ought to also be found in the Southern Hemisphere midlatitudes, whereas the simulated cooling is more moderate.


The ENS104 mean percentage increase in LGM salinity (and DIC, ALK, PO4, etc.) due to decreasing sea level is 2.84±0.62 %S, and the range is 2 %S to 4 %S. There is no significant relationship between %S and ΔCO2 in either ENS104 or ENS315, and the distribution of %S is similar in ENS104 and in both ENS315 and ENS16 (Fig. 4).

Figure 4Percentage increase in LGM salinity due to decreasing sea level distribution.


Sea ice

The ENS104 mean LGM global annual sea ice area anomaly (ΔSIA) is 18.6±7.4 million km2, and the range is 9.9 to 44 million km2. There is a negative correlation between ΔSIA and ΔSAT (r=-0.97) and between ΔSIA and ΔCO2 (r=-0.74). The negative correlation between ΔSIA and ΔCO2 likely reflects the impact of changing atmospheric CO2 on ΔSIA but may also include a smaller contribution from changing sea ice area to ΔCO2. Increasing LGM sea ice area could, for instance, have capped the outgassing of CO2 from the ocean, particularly in the Southern Ocean, and also reduced the net ocean–atmosphere CO2 flux by decreasing the AMOC strength (see below). The negative correlation between ΔSIA and ΔSAT, and between ΔSIA and ΔCO2, is reproduced in ENS315 (r=-0.96 and r=-0.74, respectively). ΔSIA in ENS104 also tends to be higher than in ENS315 and smaller than in ENS16 (Fig. 5).

Figure 5LGM change in global sea ice area distribution.


As shown in Fig. 6, fractional sea ice cover increases in all regions where sea ice is present in pre-industrial simulations, although the largest increases take place in the North Atlantic.

Figure 6LGM-PRE (a–b) and PRE (c–d) fractional sea ice cover ENS104 means (a, c) and standard deviations (b, d).



The ENS104 mean spatial distribution of the LGM precipitation rate anomaly (ΔPP) is shown in Fig. 7. The LGM changes are mostly negative but regions of positive ΔPP do exist, notably over Siberia and Australia. The largest LGM precipitation decreases (>1.5 mm d−1, with a maximum of 2.25 mm d−1) are found over northern North America and from around the eastern North Atlantic to northwest Asia, coinciding with the location of the Laurentide and Eurasian ice sheets (and the largest increases in fractional sea ice cover), respectively. Relatively large precipitation decreases (>0.75 mm d−1) are also simulated in eastern Asia, equatorial Africa and other regions of enhanced fractional sea ice cover. Comparison against a pollen-based precipitation reconstruction (Bartlein et al., 2011; Alder and Hostetler, 2015) suggest that the simulated precipitation changes over Europe and equatorial Africa are of the right direction, while precipitation changes over western Siberia at least ought to be negative. The sign of the precipitation changes over North America is mostly consistent with observations, which record negative changes over most of the continent. However, positive changes, which are also observed, are not captured. Although not shown here, comparison of the ENS104 mean against the ENS315 mean suggests that the precipitation patterns in the two are very similar, but the decreases generally tend to be higher in the ENS104 mean. The precipitation decreases in ENS104 conversely tend to be smaller than in ENS16.

Figure 7LGM change in precipitation rate (mm d−1) ENS104 mean (a) and standard deviation (b).


Ocean circulation

The ENS104 mean LGM AMOC strength anomaly (Δψmax) is -2.8±2.8 Sv, and the range is −8 to 4.7 Sv (Fig. 8). These estimates lie at the low end of the Δψmax predicted by nine Paleoclimate Modelling Intercomparison Project phase 2 (PMIP2) (Weber et al., 2007) and eight PMIP3 (Muglia and Schmittner, 2015) coupled model simulations. However, they do not include the more negative Δψmax predicted by Völker and Köhler (2013), for instance. The ENS104 mean LGM-PRE AABW cell strength in the Atlantic Ocean (Δψmin) is 0.1±1.2 Sv. A positive Δψmin represents an LGM decrease in cell strength as we keep the original (negative) sign for anticlockwise flow of Antarctic water. A negative Δψmin conversely represents an LGM increase in cell strength. However, the difference between the LGM and PRE Atlantic AABW here is not statistically significant. The range of Δψmin is −4.3 to 4.3 Sv, roughly comparable to the range of Δψmin predicted in Weber et al. (2007) (see also Muglia and Schmittner, 2015) but excluding the much larger ψmin increase predicted by Kim et al. (2003), for example. As shown in Fig. 9, the northern limit of the ENS104 mean LGM AABW cell is roughly at the same latitude as in the pre-industrial simulations. The maximum depth reached by the ensemble mean AMOC base is also similar to pre-industrial. Observations (Lynch-Stieglitz et al., 2007; Lippold et al., 2012; Gebbie, 2014; Böhm et al., 2015), conversely, suggest that the LGM AMOC shoaled to less than 2 km, raising its base depth by 2500 and 600 m at the north and south ends of the return flow, respectively. The LGM AABW, in turn, is thought to have filled the deep Atlantic below 2 km, reaching as far north as 65 N, which is approximately 25 north of its modern northern limit (Oppo et al., 2015).

Figure 8LGM change in ψmax and ψmin (a–b) distributions.


Figure 9LGM (a–b) and PRE (c–d) Atlantic overturning stream function (Sv) ENS104 means (a, c) and standard deviations (b, d).


Although not shown here, the ENS16 ensemble members tend to exhibit a shoaling of the AMOC and enhanced penetration of AABW. With regard to Δψmax and Δψmin, one can see from Fig. 8 that these tend to be more negative (i.e. weaker AMOC and stronger AABW) than in ENS104. The Δψmax and Δψmin in ENS315 tend to conversely be more positive. In ENS104, we also find a positive relationship between Δψmax and ΔCO2 (r=0.57) and a negative relationship between Δψmin and ΔCO2 (r=-0.42). The relationships are reproduced in ENS315 (r=0.59 and r=-0.36, respectively). We additionally find, in both ENS104 and ENS315, negative correlations between Δψmax and Δψmin (r=-0.62 and −0.63), Δψmin and ΔSAT (r=-0.4 and −0.4), and Δψmax and ΔSIA (r=-0.62 and −0.66), as well as positive correlations between Δψmax and ΔSAT (r=0.68 and 0.66), and Δψmin and ΔSIA (r=0.37 and 0.42). Based on these relationships, we hypothesise that increasing LGM AABW strength led to an expansion of the AABW cell. The latter in turn restricted the AMOC to lower depths and reduced its overturning rate (e.g. Shin et al., 2003). The increase in AABW strength was likely driven by increases in sea ice enhancing brine rejection. Sea ice increases in the North Atlantic may have additionally weakened the AMOC cell by locally reducing deep convection.

Table 3LGM-PRE carbon partitioning scenarios in ENS315, ENS104 and ENS16.

Download Print Version | Download XLSX

The relationships between Δψmax, Δψmin and ΔCO2 are also consistent with increasing ψmin (decreasing ψmax) contributing to decreasing atmospheric CO2. The replacement of North Atlantic Deep Water (NADW) by AABW in the North Atlantic would, for instance, have led to a dissolution of deep-sea sediment CaCO3 due to AABW having a lower bottom water CO32- concentration than NADW (see, e.g. Yu et al., 2014). The increased CaCO3 dissolution flux would in turn have raised the whole ocean alkalinity, lowering the atmospheric CO2. Enhanced AABW production would also have caused the deep ocean to become more stratified, allowing more DIC to accumulate at depth and promoting further CaCO3 dissolution. A decrease in NADW formation could have additionally lowered atmospheric CO2 by reducing CO2 outgassing at the ocean surface and reducing the burial rate of deep-sea CaCO3 due to the concomitant increase in deep-sea DIC accumulation. Further investigation is, however, required to confirm these causal relationships.

3.2.2 Terrestrial biosphere, ocean and lithospheric carbon

As shown in Fig. 10, most of the ensemble members in ENS104 predict an LGM increase in terrestrial biosphere (ΔTerrC) and lithospheric1 (ΔLithC) carbon inventory and a decrease in ocean carbon inventory (ΔOceanC). The remaining ensemble members predict one of four other scenarios of carbon partitioning, with the second most common scenario (11 % of ensemble members) being increasing terrestrial carbon and decreasing ocean and lithospheric carbon (Table 3). Similar patterns can also be observed in ENS315 and ENS16. A likely explanation for scenario 1 (increase in terrestrial biosphere and lithospheric carbon, decrease in ocean carbon) is that reduced soil decomposition (see below) causes a flux of CO2 from the atmosphere to the land, leading to an immediate outgassing of CO2 from the ocean to remove the atmospheric pCO2 difference. The CO2 outgassing also leads to an increase in surface [CO32-] and subsequently deep ocean [CO32-], which reduces CaCO3 dissolution (and increases lithospheric carbon). The increase in CaCO3 burial in turn decreases [CO32-] and increases [CO2], which is communicated back to the surface, with a resultant increase in atmospheric CO2 (Kohfeld and Ridgwell, 2009). The above explanation is of course only part of the explanation for this dominant carbon partitioning scenario, with physical mechanisms also expected to play a role, in addition to any changes in ocean productivity and changes in land carbonate weathering (see below).

Figure 10LGM change in vegetation (a), soil (b), terrestrial (vegetation plus soil) (c), ocean (d) and lithospheric (e) carbon inventory distributions.


The ENS104 mean ΔTerrC, ΔOceanC and ΔLithC, the signs of which are consistent with scenario 1, are reported in Table 4, alongside previous estimates from observational data- and model-based studies. From here, we can see that the mean ΔTerrC is only aligned with a handful of estimates and no studies so far report a negative ΔOceanC. Instead, ΔOceanC is estimated to be positive, primarily based on carbon isotope data. The loss of hundreds of petagrams of carbon from the ocean in response to terrestrial carbon growth has, however, been previously proposed (e.g. Zimov et al., 2006). Moreover, if we assume that 90 % of the atmospheric CO2 perturbation caused by the increase in terrestrial biosphere carbon reported in Table 4 gets removed by the ocean and sediments, the change in ocean carbon would be negative, even after adding the remaining carbon to be lost from the atmosphere to the ocean. We discuss what these results would likely mean for carbon isotope data in Sect. 3.2.6.

Table 4LGM-PI difference in terrestrial (ΔTerrC), ocean (ΔOceanC) and lithospheric (ΔLithC) carbon inventory (PgC) in this study (ENS104 mean, standard deviation and range) and previous studies. CR95 is Crowley (1995), AF98 is Adams and Faure (1998), Z03 is Zeng (2003), ZI09 is Zimov et al. (2009), PR11 is Prentice et al. (2011), ZE11 is Zech et al. (2011), BR12 is Brovkin et al. (2012), CI12 is Ciais et al. (2012), GL13 is Goodwin and Lauderdale (2013), OA13 is O'ishi and Abe-Ouchi (2013), SA is Sarnthein et al. (2013), PE14 is Peterson et al. (2014), AL15 is Allen et al. (2015), BG15 is Brovkin and Ganopolski (2015), SK15 is Skinner et al. (2015), SS16 is Schmittner and Somes (2016), ME17 is Menviel et al. (2017) and JE18 is Jeltsch-Thömmes et al. (2019).

Download Print Version | Download XLSX

The positive ΔTerrC studies in Table 4 attribute the increase in terrestrial carbon to different factors: Zimov et al. (2009) and Zech et al. (2011) predict large increases in permafrost carbon, while Zeng (2003) ignores permafrost. Instead, the author attributes the glacial terrestrial carbon increase to the preservation rather than destruction of carbon in areas to be covered by the LGM ice sheets, lower soil respiration rates (in the active carbon pool) caused by a colder climate, as well as storage of carbon on exposed continental shelves. Here, neither the latter carbon accumulation mechanism, nor that of permafrost growth, are included. However, our model does attempt to capture the very slow rates of soil decomposition characteristic of permafrost (Williamson et al., 2006). We attribute the terrestrial carbon increase in our ensemble to a higher soil carbon inventory caused by a decreasing soil respiration rate. As can be seen from Fig. 10, vegetation carbon tends to conversely decrease at the LGM. As in Zeng (2003), we also do not assume that, as the ice sheets expand over terrestrial carbon, they destroy it. The lack of this potential loss term means that our ΔTerrC estimates may be higher than in many previous studies, irrespective of what the response of the terrestrial biosphere is to the LGM climate and CO2 forcings. O'ishi and Abe-Ouchi (2011), for instance, estimated that the LGM climate is responsible for the loss of 502 PgC of terrestrial carbon, while another 388 PgC is removed by the ice sheets. Zeng (2003) conversely proposed that 431 PgC is preserved under the ice sheets at the LGM. This number includes 315 PgC present during the interglacial and another 116 PgC accumulated in response to the glacial climate forcings, prior to insulation of the terrestrial carbon from the atmosphere by the ice sheet coverage. Here, analysis of ENS16 suggests that during the 1000 years of LGM ice sheet build-up, the terrestrial carbon inventory in the areas to be occupied by the ice sheets increases by between 6 and 444 PgC, yielding LGM “ice sheet or burial” carbon inventories between 318 and 1341 PgC (Table 5). This increase accounts for less than half of the total LGM change in terrestrial carbon (i.e. ΔTerrC) in the majority of simulations. However, if this burial carbon were to have been destroyed rather than preserved, ΔTerrC would be negative in all but three simulations, as opposed to positive in all but one simulation (Table 5).

Table 5Ice sheet and non-ice-sheet terrestrial carbon stocks in ENS16. Columns 2 and 5 show the amount of carbon stored in areas covered by the Eurasian and Laurentide ice sheets (“ice sheet”) during the pre-industrial and LGM periods, respectively. Column 3 is the difference between the two inventories. Column 4 is the LGM change in carbon in ice sheet areas expressed as a percentage of the total LGM terrestrial carbon change. Column 6 is the LGM change in carbon outside of the ice sheets.

Download Print Version | Download XLSX

Most of the terrestrial carbon increase in areas to be covered by the ice sheets in ENS16 is due to soil carbon, with vegetation carbon decreasing in all but one simulation. The range of terrestrial carbon increases and the associated burial carbon amounts include Zeng (2003)'s estimates. Our LGM burial carbon estimates would also accommodate an additional 250–550 PgC (Franzen, 1994) from increased glacial peat accumulation (Zeng, 2003). No observational data-based estimates of the LGM burial carbon inventory are available since only limited evidence exists for organic material being preserved by ice during glaciations (Franzen, 1994, and references in Weitemeyer and Buffett, 2006). Outside of the ice sheets, increases in the terrestrial carbon inventory in ENS16 are mostly due to soil carbon, which increases in all simulations. Vegetation carbon, conversely, decreases in the majority of simulations. Our range of carbon changes outside of the ice sheet areas include the 198 PgC increase predicted by Zeng (2003) as a result of reduced soil respiration.

Although not evaluated directly, it is likely that similar ice sheet/non-ice-sheet terrestrial carbon proportions than in ENS16 are found in ENS104 and ENS315 because of the similar climate change distributions in all three instances (see earlier sections). Although not shown here, the spatial distribution of ΔTerrC in ENS16 is also similar to that of the ENS104 (and ENS315) mean. The spatial distribution of ΔTerrC in ENS104 is shown in Fig. 11.

Figure 11LGM vegetation (a–b), soil (c–d) and total terrestrial carbon changes (e–f) ENS104 mean (a, c) and standard deviation (b, d). Units are kgC m−2.


The largest increases in terrestrial carbon (≥20 kgC m−2) are found in North America and Europe/western Asia, both within and south of the Laurentide and Eurasian ice sheet margins (Fig. 11). Regions with smaller but still relatively large (≥10 kgC m−2) increases include the Andes and Patagonia regions, the southern tip of the African continent, eastern north Siberia and the grid cells just south of the Tibetan Plateau. The largest LGM decreases in terrestrial carbon (≥10 kgC m−2) conversely tend to be found in northwest North America, Beringia and the Tibetan Plateau region. Other regions with relatively large (≥5 kgC m−2) decreases include equatorial Africa and the deserts in central Asia. Everywhere else the LGM terrestrial carbon density increases by between 0 and 10 kgC m−2. Comparison against paleoecological reconstruction studies (Crowley, 1995) suggests that the simulated terrestrial carbon changes within the Laurentide and Eurasian ice sheet areas are of the wrong sign, except in northwest North America, since these studies assume the complete destruction of vegetation and soils in ice sheet areas. Discrepancies between the ENS104 and observations further arise from the rainforest regions, where the ensemble mean predicts terrestrial biosphere carbon density changes between −5 and 10 kgC m−2, well above observed changes of -23 kgC m−2. It is important to note, however, that as suggested in Zeng (2007), the rate of decomposition of soil carbon at the LGM may have been slower than assumed in pollen data-based studies. The largest increases in terrestrial carbon density (∼40 kgC m−2) produced by the ensemble mean are comparable to those found in areas with permafrost growth (Zimov et al., 2006). However, the peaks are potentially misplaced, being located within and south of the Laurentide and Eurasian ice sheet covered areas, rather than in eastern Siberia and Alaska. Alternatively, terrestrial carbon increases in eastern Siberia and Alaska are simply underestimated in the ensemble mean and large increases in terrestrial carbon indeed took place within the ice sheet areas during glacial periods.

The large LGM decreases in terrestrial carbon in northwest North America and adjacent Beringia are likely caused by precipitation decreasing comparatively more than SAT and causing the decrease in photosynthesis to exceed the decrease in soil respiration. However, it is also noteworthy that, although not shown here, the regions with the largest decreases in terrestrial carbon density, namely northwest North America, Beringia and the Tibetan Plateau area, are also the regions with the largest terrestrial carbon densities in the pre-industrial ENS315 mean. We further note that the Tibetan soil carbon peak is overestimated in the latter, and the North American soil carbon peak misplaced, compared to observations. We attribute the first discrepancy to the lack of soil weathering in the model and the inclusion of land use effects in the observational data-based estimate (Holden et al., 2013b; Williamson et al., 2006). The second discrepancy is attributed to the lack of explicit representation of permafrost and the absence of moisture control on soil respiration (Williamson et al., 2006).

3.2.3 Ocean primary productivity

The ENS104 mean LGM total POC export flux anomaly (ΔPOCexp) is -0.19±1 PgC yr−1 and the range is −2.57 to 2.56 PgC yr−1, roughly consistent with previous model-based estimates (e.g. Brovkin et al., 2002, 2007; Bopp et al., 2003; Chikamoto et al., 2012; Palastanga et al., 2013; Schmittner and Somes, 2016; Buchanan et al., 2016). As shown in Fig. 12, compared to ENS104, the POC flux decreases in ENS315 and ENS16 tend to be smaller and larger, respectively. ΔPOCexp is positively correlated with Δψmax (r=0.72 and 0.79) and negatively correlated with Δψmin (r=-0.62 and −0.58) in both ENS104 and ENS315. The correlations potentially suggest that decreasing AMOC strength and increasing AABW production led to decreasing POC export. One possible mechanism is enhanced deep ocean stratification due to increasing AABW formation leading to not only more efficient trapping of DIC at depth (see above) but also nutrients and therefore reduced availability in the euphotic zone. All else held constant, a weaker and shallower AMOC cell would also inhibit the transfer of nutrients from the deep ocean to the surface. A negative correlation can additionally be found between ΔPOCexp and ΔSIA (r=-0.55 and −0.6), probably because no primary production occurs beneath the sea ice surface. Increasing sea ice area at the LGM therefore leads to decreasing POC export flux. This would also explain the largest ENS104 mean decreases in POC export flux, shown in Fig. 13, coinciding with increases in sea ice fraction.

Figure 12LGM change in POC export flux distributions.


Figure 13LGM surface POC export flux change (molC m−2 yr−1) ENS104 mean (a) and standard deviation (b).


The largest ENS104 LGM increases in POC export flux conversely occur at around 50 S, roughly in front of the Antarctic sea ice margins. Increases in POC export are also simulated close to the North Pacific and Atlantic sea ice margins, as well as in the eastern equatorial Pacific and the southwest Atlantic upwelling region. The increases in POC export flux at the sea ice margins are likely caused by the advection of unutilised nutrients from underneath the sea ice. However, they may additionally be due to the enhanced iron availability from the increased supply of aeolian dust, particularly in the Southern Ocean and North Pacific since these are strongly limited by iron (Ridgwell et al., 2007). Iron fertilisation may also explain the increases in POC export flux in the eastern equatorial Pacific and in the southwest Atlantic upwelling region. Comparison against observations suggests that the ensemble mean POC flux changes immediately north, and south of the Antarctic sea ice margins align with observations of increased and reduced marine productivity in the subantarctic (∼45 to 60 N) and Southern Ocean, respectively (Kohfeld, 2005; Kohfeld et al., 2013; Jaccard et al., 2013; Martínez-García et al., 2014). The simulated decreases in export flux in the Arctic and subarctic Atlantic (i.e. above approximately 50 N), and the increases in export flux immediately south of 50 N are also in agreement with previous reconstructions (Kohfeld, 2005; Radi and de Vernal, 2008). The mostly lower LGM export fluxes at the Equator and in the South Atlantic are conversely inconsistent with the observational data of Kohfeld (2005). The decreases may be caused by the increases in productivity in high-nutrient, low-chlorophyll (HNLC) regions reducing the phosphate (the other limiting nutrient in GENIE-1 besides iron) availability for photosynthesis in other regions. They may additionally be due to the model not simulating enhanced nutrient inventories in response to enhanced weathering or reduced shallower water deposition of organic matter. The model also does not vary wind speed, which may have resulted in stronger tropical upwelling in the Atlantic at the LGM. The evidence is more ambiguous (or missing) for the Pacific (Jaccard et al., 2010; Kohfeld and Chase, 2011; Kohfeld, 2005; Costa et al., 2016) and Indian oceans (Kohfeld, 2005; Singh et al., 2011) and is therefore not discussed in more detail here.

Figure 14Land-to-ocean bicarbonate flux scaling factor (GWS) distributions.


3.2.4 Carbonate weathering and shallow water deposition

The ENS104 mean land-to-ocean bicarbonate flux scaling factor (GWS) is 1.16±0.24 (corresponding to a percentage change in the land-to-ocean bicarbonate flux, %LOC, of 38.67), and the range is 0.52 to 1.5 (corresponding to a %LOC between −49.33 and 50). As shown in Fig. 14, the GWS in ENS104 tends to be larger than in ENS315 and smaller than in ENS16. There is also a negative correlation between GWS and ΔCO2 (r=-0.52) in ENS315, suggesting that increasing the input of bicarbonate to the ocean leads to a decrease in CO2 by raising the inventories of ALK and DIC in a 2:1 ratio. In ENS104, however, r is below the 0.05 significance level, suggesting that it is less important.

Figure 15LGM change in deep-sea sediment CaCO3 burial flux distributions.


Figure 16LGM deep-sea CaCO3 burial rate change (mol cm−2 yr−1) ENS104 mean (a) and standard deviation (b).


3.2.5 Deep-sea carbonate burial

The ENS104 mean global deep-sea CaCO3 burial flux anomaly (ΔCaCO3bur) is 0.036±0.045 PgC yr−1 and the range is −0.098 to 0.139 PgC yr−1. The mean value is approximately 3 times larger than the observed value (Catubig et al., 1998), although the latter still falls within the range of simulated values. As shown in Fig. 15, ΔCaCO3bur in ENS104 tends to be higher than in ENS315 and lower than in ENS16. The change in global deep-sea CaCO3 burial flux anomaly is strongly determined by GWS, as suggested by the positive correlation (r=0.88 and 0.9) between the two, in both ENS104 and ENS315. Increasing %LOC (the percentage change in the land-to-ocean bicarbonate flux) should indeed enhance the CaCO3 burial flux as increasing ALK means the deep ocean CO32- will eventually increase. The latter in turn would cause the saturation horizon to fall, allowing CaCO3 to accumulate over greater areas (which are now exposed to undersaturated waters) (Sigman and Boyle, 2000). The input of ALK to the surface ocean would also increase the rate of CaCO3 export production (enhancing the sediment deposition flux of CaCO3), since as discussed in Chikamoto et al. (2008), the latter is proportional to the production rate of POC (which is equal to the POC export flux), together with the sea surface saturation state with respect to CaCO3, in GENIE-1. There is indeed also a positive correlation between %LOC and the global change in CaCO3 export flux (r=0.27 and 0.4), and between the latter and ΔCaCO3bur (r=0.34 and 0.45) in both ENS104 and ENS315.

The ENS104 mean spatial distribution of ΔCaCO3bur is shown in Fig. 16. Relatively large increases in burial flux (0.5×10-5 mol cm−2 yr−1) can be found at around 50 S, in the North Pacific and to a lesser extent the North Atlantic. In other regions, the burial flux is significantly lower or negative, with the largest losses (-0.5×10-5 mol cm−2 yr−1) occurring in the North Atlantic and arctic regions. The only exception is the western North Atlantic, which exhibits a large increase in burial. A comparison of the results against the reconstructions of Catubig et al. (1998) is somewhat difficult, as the coverage is poor but overall CaCO3 burial was higher in the North Atlantic and the Pacific, and lower in the tropical and South Atlantic, and the Indian Ocean and Southern Ocean.

3.2.6 Other paleoproxies

As shown in Table 4, a frequent argument for a lower glacial terrestrial carbon inventory is the reconstructed mean glacial ocean δ13C value of approximately 0.35 ‰ lower than present due to the fact that plants discriminate against 13C during photosynthesis. In our simulations, conversely, it follows that the increase in glacial terrestrial carbon inventory would have resulted in an increase in ocean δ13C. Decreasing SSTs and increasing CaCO3 weathering would have, moreover, likely raised it further (in the first instance, by enhancing fractionation at the air–sea interface, and in the second instance, through the input of isotopically heavy weathering products). However, as noted in Zeng (2007), the interpretation of the δ13C value can be complicated by factors such as the impact of enhanced glacial carbonate ion concentrations on δ13C in foramifera shells (Lea et al., 1999). In addition, there are other processes in our model which may have counteracted at least part of the increase in ocean δ13C. These include reduced marine productivity (e.g. Zimov et al., 2009), as phytoplankton discriminate against 13C during photosynthesis, giving the marine organic carbon reservoir a low δ13C. However, we note that the sign of this impact would additionally depend on the associated changes in organic matter remineralisation and burial. Another relevant process is greater sea ice area, which can lower the ocean δ13C by reducing the air–sea gas exchange and therefore the net transfer of 13C into the ocean (Stephen and Keeling, 2000). Moreover, we propose that adding missing processes could have decreased ocean δ13C even further. These include weaker surface winds (while in our model these are fixed), again through reduced air–sea gas exchange (Menviel et al., 2015), as well as enhanced weathering and reduced deposition of organic carbon at continental margins due to lower sea levels (Wallmann, 2014). Yet, further research is required here as one recent study suggests that taking into account these processes would most likely not (the possibility is not completely ruled out) allow reconciliation of a positive ΔTerrC with the observed mean glacial ocean δ13C value (Jeltsch-Thömmes et al., 2019).

Missing processes would likely also be needed to reconcile our lower POC export fluxes with deep ocean oxygen records. These records tend to indicate there was a decrease in LGM deep ocean oxygen concentration (Jaccard et al., 2016). Lower POC export fluxes would conversely have resulted in an increase in deep ocean oxygen due to reduced oxygen consumption at depth. When observed over sufficiently large areas, lower oxygen concentrations can support the presence of an enhanced ocean carbon inventory as deoxygenation can be explained by reduced ocean ventilation (the sole input of oxygen is from the ocean surface) (Wagner and Hendy, 2015). The reduced ventilation is, in turn, assumed to have led to the accumulation of a significant amount of DIC in the ocean interior. Thus, explaining the lower deep ocean oxygen concentrations without having to reduce ocean ventilation as extensively as suggested by previous studies would as a minimum likely require LGM export production to have increased rather than decreased. However, it may be possible to increase deep ocean oxygen consumption by increasing organic matter at depth but keeping the surface POC export flux constant. This would require adding missing processes such as increasing remineralisation depth with decreasing ocean temperature and increasing ballasting into our model (Kohfeld and Ridgwell, 2009; Menviel et al., 2012).

4 Conclusions

We have used an uncertainty-based approach to investigating the LGM atmospheric CO2 drop by simulating it with a large ensemble of parameter sets and exploring the range of possible responses. Despite our ensemble varying many of the parameters thought to contribute to variability in glacial–interglacial atmospheric CO2, we estimated that up to  60 ppmv of ΔCO2 could be attributed to processes not included in our model and error in our process representations. As a result, we treated ΔCO2 between -90 and −30 ppmv as equally plausible and focused on describing the responses of the subset of simulations with this ΔCO2. We found the range of responses to be large, including the presence of five different ways of achieving a plausible ΔCO2 in terms of the sign of individual carbon reservoir changes. However, several dominant changes could be detected. Namely, the LGM atmospheric CO2 decrease tended to predominantly be associated with decreasing SSTs, increasing sea ice area, a weakening of the AMOC, a strengthening of the AABW cell in the Atlantic Ocean, a decreasing ocean biological productivity, an increasing CaCO3 weathering flux and an increasing deep-sea CaCO3 burial flux. The majority of our simulations also predicted an increase in terrestrial carbon, coupled with a decrease in ocean and an increase in lithospheric carbon. The increase in terrestrial carbon, which is uncommon in LGM simulations, was attributed to reduced soil respiration in response to the climate forcings, as well as our choice to preserve rather than destroy carbon that accumulates in ice sheet areas. The dominant changes were broadly in agreement with observations and paleoproxies other than carbon isotope and oxygen data, which we did not evaluate directly. However, we advise more detailed comparisons in future studies. It is also likely that our results can only be reconciled with carbon isotope and oxygen data if processes currently missing from our model are taken into account.

Code and data availability

GENIE-1 was checked out via (last access: May 2019) using Subversion (SVN). The simulations described here are with release version 2-8-0. In addition to the source code, several packages and applications such as the NetCDF libraries are required by GENIE-1 (University of Bristol Geography Source, 2014).

The way in which GENIE-1 is run manually is as described in Ridgwell (2012) for the GENIE developmental variant cGENIE: the basic flavour and configuration of GENIE-1 is run from /genie/genie-main by issuing the command


genie.job is a shell script which determines the basic (“base”) configuration of the model. A different flavour and configuration of the model is obtained by specifying a different base configuration file:

./genie.job -f example.xml

where example.xml is a specified model configuration (/flavour) .xml file (Ridgwell, 2012).

The ensemble parameter sets required to repeat the experiments in this study are available upon request to the corresponding author (

Appendix A: The OLR feedback parameter

The ensemble parameter OL1 is varied through all stages but stage 1. OL1 describes the unmodelled response of clouds to global average temperature change and corresponds to the KLW1 constant in the equation below (Eq. 1 in Holden et al., 2010a):

(A1) L out = L out ( T , q ) - K LW 0 - K LW 1 Δ T ,

where Lout(T,q) is the unmodified “clear skies” OLR term of Thompson and Warren (1982), KLW0 is the clear-sky outgoing long-wave radiation parameter (OL0), representing the effects of clouds on the unmodified OLR (with KLW0 only ever taking positive values), and ΔT corresponds to the difference between the globally averaged surface air temperature and the equilibrium pre-industrial temperature (Holden et al., 2010a).

In stage 1, OL1 is set to zero and T0 can be set to any value. The temperatures simulated at the end of stage 1 are used to define T0 in stage 2 and in the LGM simulations.


The supplement related to this article is available online at:

Author contributions

KMSK designed the experiments with guidance from PBH, NRE and AR. KMSK carried them out and analysed the data. KMSK prepared the manuscript with contributions from all co-authors.

Competing interests

The authors declare that they have no conflict of interest.


This work made use of the Darwin Supercomputer of the University of Cambridge High Performance Computing Service (HPCS). Krista M. S. Kemppinen thanks staff at HPCS for their technical support, and Antara Banerjee and Alex Archibald for help with R. The manuscript was greatly improved by comments from an anonymous reviewer.

Financial support

This research has been supported by the UK Natural Environment Research Council (NERC) through funding for the project DESIRE (grant no. NE/E007554/1).

Review statement

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


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

Adkins, J. F., Mcintyre, K., and Schrag, D. P.: The Salinity, Temperature, and δ18O of the Glacial Deep Ocean, Science, 298, 1769–1773,, 2002. 

Alder, J. R. and Hostetler, S. W.: Global climate simulations at 3000-year intervals for the last 21 000 years with the GENMOM coupled atmosphere–ocean model, Clim. Past, 11, 449–471,, 2015. 

Allen, K. A., Sikes, E. L., Hönisch, B., Elmore, A. C., Guilderson, T. P., Rosenthal, Y., and Anderson, R. F.: Southwest Pacific deep water carbonate chemistry linked to high southern latitude climate and atmospheric CO2 during the Last Glacial Termination, Quaternary Sci. Rev., 122, 180–191,, 2015. 

Anderson, R. F., Ali, S., Bradtmiller, L. I., Nielsen, S. H. H., Fleisher, M. Q., Anderson, B. E., and Burckle, L. H.: Wind-Driven Upwelling in the Southern Ocean and the Deglacial Rise in Atmospheric CO2, Science, 323, 1443–1448,, 2009. 

Annan, J. D. and Hargreaves, J. C.: Efficient identification of ocean thermodynamics in a physical/biogeochemical ocean model with an iterative Importance Sampling method, Ocean Model., 32, 205–215,, 2010. 

Annan, J. D. and Hargreaves, J. C.: A new global reconstruction of temperature changes at the Last Glacial Maximum, Clim. Past, 9, 367–376,, 2013. 

Ballarotta, M., Falahat, S., Brodeau, L., and Döös, K.: On the glacial and interglacial thermohaline circulation and the associated transports of heat and freshwater, Ocean Sci., 10, 907–921,, 2014. 

Bartlein, P. J., Harrison, S. P., Brewer, S., Connor, S., Davis, B. A. S., Gajewski, K., Guiot, J., Harrison-Prentice, T. I., Henderson, A., Peyron, O., Prentice, I. C., Scholze, M., Seppä, H., Shuman, B., Sugita, S., Thompson, R. S., Viau, A. E., Williams, J., and Wu, H.: Pollen-based continental climate reconstructions at 6 and 21 ka: A global synthesis, Clim. Dynam., 37, 775–802,, 2011. 

Berger, A.: Long-term variations of daily insolation and quaternary climatic changes, J. Atmos. Sci., 35, 2362–2367, 1978. 

Böhm, E., Lippold, J., Gutjahr, M., Frank, M., Blaser, P., Antz, B., Fohlmeister, J., Frank, N., Andersen, M. B., and Deininger, M.: Strong and deep Atlantic meridional overturning circulation during the last glacial cycle, Nature, 517, 73–76,, 2015. 

Bopp, L., Kohfeld, K. E., Le Quéré, C., and Aumont, O.: Dust impact on marine biota and atmospheric CO2 during glacial periods, Paleoceanography, 18, 1046,,, 2003. 

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

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

Brady, E. C., Otto-Bliesner, B. L., Kay, J. E., and Rosenbloom, N.: Sensitivity to glacial forcing in the CCSM4, J. Climate, 26, 1901–1925,, 2013. 

Braconnot, P., Otto-Bliesner, B., Harrison, S., Joussaume, S., Peterchmitt, J.-Y., Abe-Ouchi, A., Crucifix, M., Driesschaert, E., Fichefet, Th., Hewitt, C. D., Kageyama, M., Kitoh, A., Laîné, A., Loutre, M.-F., Marti, O., Merkel, U., Ramstein, G., Valdes, P., Weber, S. L., Yu, Y., and Zhao, Y.: Results of PMIP2 coupled simulations of the Mid-Holocene and Last Glacial Maximum – Part 1: experiments and large-scale features, Clim. Past, 3, 261–277,, 2007. 

Brovkin, V. and Ganopolski, A.: The role of the terrestrial biosphere in CLIMBER-2 simulations of the last glacial CO2 cycles, Nova Act. LC NF, 121, 43–47, 2015. 

Brovkin, V., Hofmann, M., Bendtsen, J., and Ganopolski, A.: Ocean biology could control atmospheric δ13C during glacial-interglacialcycle, Geochem. Geophy. Geosy., 3, 1027,, 2002. 

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

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. 

Buchanan, P. J., Matear, R. J., Lenton, A., Phipps, S. J., Chase, Z., and Etheridge, D. M.: The simulated climate of the Last Glacial Maximum and insights into the global marine carbon cycle, Clim. Past, 12, 2271–2295,, 2016. 

Catubig, N. R., Archer, D. E., Francois, R., DeMenocal, P., Howard, W., and Yu, E. F.: Global deep-sea burial rate of calcium carbonate during the Last Glacial Maximum, Paleoceanography, 13, 298–310,, 1998. 

Chikamoto, M. O., Matsumoto, K., and Ridgwell, A.: Response of deep-sea CaCO3 sedimentation to Atlantic meridional overturning circulation shutdown, J. Geophys. Res.-Biogeo., 113, G03017,, 2008. 

Chikamoto, M. O., Abe-Ouchi, A., Oka, A., Ohgaito, R., and Timmermann, A.: Quantifying the ocean's role in glacial CO2 reductions, Clim. Past, 8, 545–563,, 2012. 

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. 

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

Colbourn, G.: Weathering effects on the carbon cycle in an Earth System Model, PhD thesis, School of Environmental Sciences, University of East Anglia, UK, 2011. 

Costa, K. M., McManus, J. F., Anderson, R. F., Ren, H., Sigman, D. M., Winckler, G., Fleisher, M. Q., Marcantonio, F., and Ravelo, A. C.: No iron fertilization in the equatorial Pacific Ocean during the last ice age, Nature, 529, 519–522,, 2016. 

Crocket, K. C., Vance, D., Foster, G. L., Richards, D. A., and Tranter, M.: Continental weathering fluxes during the last glacial/interglacial cycle: insights from the marine sedimentary Pb isotope record at Orphan Knoll, NW Atlantic, Quaternary Sci. Rev., 38, 89–99,, 2012. 

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

De La Fuente, M., Skinner, L., Calvo, E., Pelejero, C., and Cacho, I.: Increased reservoir ages and poorly ventilated deep waters inferred in the glacial Eastern Equatorial Pacific, Nat. Commun., 6, 7420,, 2015. 

Doney, S. C., Lindsay, K., Fung, I., and John, J.: Natural variability in a stable, 1000-yr global coupled climate-carbon cycle simulation, J. Climate, 19, 3033–3054,, 2006. 

Edwards, N. R. and Marsh, R.: Uncertainties due to transport-parameter sensitivity in an efficient 3-D ocean-climate model, Clim. Dynam., 24, 415–433,, 2005. 

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

Foster, G. L. and Vance, D.: Negligible glacial-interglacial variation in continental chemical weathering rates, Nature, 444, 918–921,, 2006. 

Franzen, L. G.: Are wetlands the key to the ice-age cycle enigma, Ambio, 23, 300–308, 1994. 

Franzén, L. G. and Cropp, R. A.: The Peatland/ice age Hypothesis revised, adding a possible glacial pulse trigger, Geogr. Ann. A, 89, 301–330, 2007. 

Freeman, E., Skinner, L. C., Tisserand, A., Dokken, T., Timmermann, A., Menviel, L., and Friedrich, T.: An Atlantic-Pacific ventilation seesaw across the last deglaciation, Earth Planet. Sc. Lett., 424, 237–244,, 2015. 

Gebbie, G.: How much did Glacial North Atlantic Water shoal?, Paleoceanography, 29, 190–209,, 2014. 

Goodwin, P. and Lauderdale, J. M.: Carbonate ion concentrations, ocean carbon storage, and atmospheric CO2, Global Biogeochem. Cy., 27, 882–893,, 2013. 

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

Holden, P. B., Edwards, N. R., Oliver, K. I. C., Lenton, T. M., and Wilkinson, R. D.: A probabilistic calibration of climate sensitivity and terrestrial carbon change in GENIE-1, Clim. Dynam., 35, 785–806,, 2010a. 

Holden, P. B., Edwards, N. R., Wolff, E. W., Lang, N. J., Singarayer, J. S., Valdes, P. J., and Stocker, T. F.: Interhemispheric coupling, the West Antarctic Ice Sheet and warm Antarctic interglacials, Clim. Past, 6, 431–443,, 2010b. 

Holden, P. B., Edwards, N. R., Müller, S. A., Oliver, K. I. C., Death, R. M., and Ridgwell, A.: Controls on the spatial distribution of oceanic δ13CDIC, Biogeosciences, 10, 1815–1833,, 2013a. 

Holden, P. B., Edwards, N. R., Gerten, D., and Schaphoff, S.: A model-based constraint on CO2 fertilisation, Biogeosciences, 10, 339–355,, 2013b. 

Iglesias-Rodríguez, M. D., Brown, C. W., Doney, S. C., Kleypas, J., Kolber, D., Kolber, Z., Hayes, P. K., and Falkowski, P. G.: Representing key phytoplankton functional groups in ocean carbon cycle models: Coccolithophorids, Global Biogeochem. Cy., 16, 47-1–47-20,, 2002. 

Jaccard, S. L., Galbraith, E. D., Sigman, D. M., and Haug, G. H.: A pervasive link between Antarctic ice core and subarctic Pacific sediment records over the past 800 kyrs, Quaternary Sci. Rev., 29, 206–212,, 2010. 

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

Jaccard, S. L., Galbraith, E. D., Martínez-García, A., and Anderson, R. F.: Covariation of deep Southern Ocean oxygenation and atmospheric CO2 through the last ice age, Nature, 530, 207–210,, 2016. 

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. 

Jones, I. W., Munhoven, G., Tranter, M., Huybrechts, P., and Sharp, M. J.: Modelled glacial and non-glacial HCO3-, Si and Ge fluxes since the LGM: Little potential for impact on atmospheric CO2 concentrations and a potential proxy of continental chemical erosion, the marine Ge∕Si ratio, Global Planet. Change, 33, 139–153,, 2002. 

Joos, F., Gerber, S., Prentice, I. C., Otto-Bliesner, B. L., and Valdes, P. J.: Transient simulations of Holocene atmospheric carbon dioxide and terrestrial carbon since the Last Glacial Maximum, Global Biogeochem. Cy., 18, GB2002,, 2004. 

Kim, S.-J., Flato, G., and Boer, G.: A coupled climate model simulation of the Last Glacial Maximum, Part 2: approach to equilibrium, Clim. Dynam., 20, 635–661,, 2003. 

Kleypas, J. A.: Modeled estimates of global reef habitat and carbonate production since the last glacial maximum, Paleoceanography, 12, 533–545,, 1997. 

Kohfeld, K. E.: Role of Marine Biology in Glacial-Interglacial CO2 Cycles, Science, 308, 74–78,, 2005. 

Kohfeld, K. E. and Chase, Z.: Controls on deglacial changes in biogenic fluxes in the North Pacific Ocean, Quaternary Sci. Rev., 30, 3350–3363,, 2011. 

Kohfeld, K. E. and Ridgwell, A.: Glacial-Interglacial Variability in Atmospheric CO2, Surf. Ocean. Atmos. Process., 187, 251–286,, 2009. 

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

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

Lea, D. W., Bijma, J., Spero, H. J., and Archer, D.: Implications of a carbonate ion effect on shell carbon and oxygen isotopes for glacial ocean conditions, in: Use of Proxies in Paleoceanography: Examples from the South Atlantic, edited by: Fischer, G. and Wefer, G., Springer-Verlag, Berlin Heidelberg, 513–522, 1999. 

Lenton, T. M., Marsh, R., Price, A. R., Lunt, D. J., Aksenov, Y., Annan, J. D., Cooper-Chadwick, T., Cox, S. J., Edwards, N. R., Goswami, S., Hargreaves, J. C., Harris, P. P., Jiao, Z., Livina, V. N., Payne, A. J., Rutt, I. C., Shepherd, J. G., Valdes, P. J., Williams, G., Williamson, M. S., and Yool, A.: Effects of atmospheric dynamics and ocean resolution on bi-stability of the thermohaline circulation examined using the Grid ENabled Integrated Earth system modelling (GENIE) framework, Clim. Dynam., 29, 591–613,, 2007. 

Lippold, J., Luo, Y., Francois, R., Allen, S. E., Gherardi, J., Pichat, S., Hickey, B., and Schulz, H.: Strength and geometry of the glacial Atlantic Meridional Overturning Circulation, Nat. Geosci., 5, 813–816,, 2012. 

Lupker, M., France-Lanord, C., Galy, V., Lavé, J., and Kudrass, H.: Increasing chemical weathering in the Himalayan system since the Last Glacial Maximum, Earth Planet. Sc. Lett., 365, 243–252,, 2013. 

Lynch-Stieglitz, J., Adkins, J. F., Curry, W. B., Dokken, T., Hall, I. R., Herguera, J. C., Hirschi, J. J.-M., Ivanova, E. V, Kissel, C., Marchal, O., Marchitto, T. M., McCave, I. N., McManus, J. F., Mulitza, S., Ninnemann, U., Peeters, F., Yu, E.-F., and Zahn, R.: Atlantic meridional overturning circulation during the Last Glacial Maximum, Science, 316, 66–69,, 2007. 

Ma, W. and Tian, J.: Modeling the contribution of dissolved organic carbon to carbon sequestration during the last glacial maximum, Geo-Mar. Lett., 34, 471–482,, 2014. 

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

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

Martin, P., Archer, D., and Lea, D. W.: Role of deep sea temperature in the carbon cycle during the last glacial, Paleoceanography, 20, 1–10,, 2005. 

Martínez-García, A., Sigman, D. M., Ren, H., Anderson, R. F., Straub, M., Hodell, D. A., Jaccard, S. L., Eglinton, T. I., and Haug, G. H.: Iron Fertilization of the Subantarctic Ocean During the Last Ice Age, Science, 343, 1347–1350,, 2014. 

Masson-Delmotte, V., Kageyama, M., Braconnot, P., Charbit, S., Krinner, G., Ritz, C., Guilyardi, E., Jouzel, J., Abe-Ouchi, A., Crucifix, M., Gladstone, R. M., Hewitt, C. D., Kitoh, A., LeGrande, A. N., Marti, O., Merkel, U., Motoi, T., Ohgaito, R., Otto-Bliesner, B., Peltier, W. R., Ross, I., Valdes, P. J., Vettoretti, G., Weber, S. L., Wolk, F., and Yu, Y.: Past and future polar amplification of climate change: Climate model intercomparisons and ice-core constraints, Clim. Dynam., 26, 513–529,, 2006. 

Masson-Delmotte, V., Schulz, M., Abe-Ouchi, A., Beer, J., Ganopolski, A., González Rouco, J. F., Jansen, E., Lambeck, K., Luterbacher, J., Naish, T., Osborn, T., Otto-Bliesner, B., Quinn, T., Ramesh, R., Rojas, M., Shao, X., and Timmermann, A.: Information from Paleoclimate Archives, in: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M., Cambridge University Press, Cambridge, UK and New York, NY, USA, 2013. 

Matsumoto, K., Sarmiento, J. L., and Brzezinski, M. A.: Silicic acid leakage from the Southern Ocean: A possible explanation for glacial atmospheric pCO2, Global Biogeochem. Cy., 16, 5-1–5-23,, 2002. 

Matsumoto, K., Hashioka, T., and Yamanaka, Y.: Effect of temperature-dependent organic carbon decay on atmospheric pCO2, J. Geophys. Res.-Biogeo., 112, G02007,, 2007. 

Matsumoto, K., Tokos, K., Huston, A., and Joy-Warren, H.: MESMO 2: a mechanistic marine silica cycle and coupling to a simple terrestrial scheme, Geosci. Model Dev., 6, 477–494,, 2013. 

Matsumoto, K., Chase, Z., and Kohfeld, K.: Different mechanisms of silicic acid leakage and their biogeochemical consequences, Paleoceanography, 29, 238–254,, 2014. 

Menviel, L., Joos, F., and Ritz, S. P.: 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. 

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. 

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, 32, 2–17,, 2017. 

Muglia, J. and Schmittner, A.: Glacial Atlantic overturning increased by wind stress in climate models, Geophys. Res. Lett., 42, 9862–9869,, 2015. 

Munhoven, G.: Glacial – Interglacial changes of continental weathering: Estimates of the related CO2 and HCO3- flux variations and their uncertainties, Global Planet. Change, 33, 155–176,, 2002. 

O'ishi, R. and Abe-Ouchi, A.: Influence of dynamic vegetation on climate change and terrestrial carbon storage in the Last Glacial Maximum, Clim. Past, 9, 1571–1587,, 2013. 

Oka, A., Abe-Ouchi, A., Chikamoto, M. O., and Ide, T.: Mechanisms controlling export production at the LGM: Effects of changes in oceanic physical fields and atmospheric dust deposition, Global Biogeochem. Cy., 25, GB2009,, 2011. 

Opdyke, B. N. and Walker, J. C. G.: Return of the coral reef hypothesis: basin to shelf partitioning of CaCO3 and its effect on atmospheric CO2, Geology, 20, 733–736,<0733:ROTCRH>2.3.CO;2, 1992. 

Oppo, D. W., Curry, W. B., and McManus, J. F.: What do benthic δ13C and δ18O data tell us about Atlantic circulation during Heinrich Stadial 1?, Paleoceanography, 30, 353–368,, 2015. 

Otto, D., Rasse, D., Kaplan, J., Warnant, P., and François, 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. 

Palastanga, V., Slomp, C. P., and Heinze, C.: Glacial-interglacial variability in ocean oxygen and phosphorus in a global biogeochemical model, Biogeosciences, 10, 945–958,, 2013. 

Peltier, W. R.: Ice age paleotopography, Science, 265, 195–201,, 1994. 

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. 

Prentice, I., Farquhar, G., and Fasham, M.: The carbon cycle and atmospheric carbon dioxide, Weather, 183–237,, 2001. 

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

Radi, T. and de Vernal, A.: Last glacial maximum (LGM) primary productivity in the northern North Atlantic Ocean, Can. J. Earth Sci., 45, 1299–1316,, 2008. 

Ridgwell, A. and Hargreaves, J. C.: Regulation of atmospheric CO2 by deep-sea sediments in an Earth system model, Global Biogeochem. Cy., 21, GB2008,, 2007. 

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

Ridgwell, A. J.: cGENIE v.0.9 (“muffin”) User Manual [PDF document], available at:, last access: 26 October 2012. 

Roth, R., Ritz, S. P., and Joos, F.: Burial-nutrient feedbacks amplify the sensitivity of atmospheric carbon dioxide to changes in organic matter remineralisation, Earth Syst. Dynam., 5, 321–343,, 2014. 

Sarnthein, M., Schneider, B., and Grootes, P. M.: Peak glacial 14C ventilation ages suggest major draw-down of carbon into the abyssal ocean, Clim. Past, 9, 2595–2614,, 2013. 

Schmittner, A. and Somes, C. J.: Complementary constraints from carbon (13C) and nitrogen (15N) isotopes on the glacial ocean's soft-tissue biological pump, Paleoceanography, 31, 669–693,, 2016. 

Schmittner, A., Urban, N. M., Shakun, J. D., Mahowald, N. M., Clark, P. U., Bartlein, P. J., Mix, A. C., and Rosell-Mele, A.: Climate Sensitivity Estimated from Temperature Reconstructions of the Last Glacial Maximum, Science, 334, 1385–1388,, 2011. 

Schneider von Deimling, T., Held, H., Ganopolski, A., and Rahmstorf, S.: Climate sensitivity estimated from ensemble simulations of glacial climate, Clim. Dynam., 27, 149–163,, 2006. 

Shin, S.-I., Liu, Z., Otto-Bliesner, B., Brady, E., Kutzbach, J., and Harrison, S.: A Simulation of the Last Glacial Maximum climate using the NCAR-CCSM, Clim. Dynam., 20, 127–151,, 2003. 

Sigman, D. M. and Boyle, E. A.: Glacial/Interglacial Variations In Atmospheric Carbon Dioxide, Nature, 407, 859–869,, 2000. 

Simmons, C. T., Mysak, L. A., and Matthews, H. D.: An investigation of carbon cycle dynamics since the Last Glacial Maximum: Complex interactions between the terrestrial biosphere, weathering, ocean alkalinity, and CO2 radiative warming in an Earth system model of intermediate complexity, Clim. Past Discuss.,, 2016. 

Singh, A. D., Jung, S. J. A., Darling, K., Ganeshram, R., Ivanochko, T., and Kroon, D.: Productivity collapses in the Arabian Sea during glacial cold phases, Paleoceanography, 26, PA3210,, 2011. 

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. 

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

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

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

Sun, X. and Matsumoto, K.: Effects of sea ice on atmospheric pCO2: A revised view and implications for glacial and future climates, J. Geophys. Res., 115, 1–8,, 2010. 

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

Thompson S. L. and Warren S. G.: Parameterization of outgoing infrared radiation derived from detailed radiative calculations, J. Atmos. Sci., 39, 2667–2680,<2667:POOIRD>2.0.CO;2, 1982. 

Tiedemann, R., Ronge, T., Lamy, F., Köhler, P., Frische, M., De Pol Holz, R., Pahnke, K., Alloway, B. V., Wacker, L., and Southon, J.: New Constraints on the Glacial Extent of the Pacific Carbon Pool and its Deglacial Outgassing, Nova Act. LC, 121, 229–233, 2015. 

Toggweiler, J. R., Russell, J. L., and Carson, S. R.: Midlatitude westerlies, atmospheric CO2, and climate change during the ice ages, Paleoceanography, 21, PA2005,, 2006. 

Tréguer, P. and Pondaven, P.: Global change: Silica control of carbon dioxide, Nature, 486, 358–359, 2000. 

University of Bristol Geography Source: GENIE, available at: (last access: May 2019), 2014. 

Vance, D., Teagle, D. A. H., and Foster, G. L.: Variable Quaternary chemical weathering fluxes and imbalances in marine geochemical budgets, Nature, 458, 493–496,, 2009. 

Völker, C. and Köhler, P.: Responses of ocean circulation and carbon cycle to changes in the position of the Southern Hemisphere westerlies at Last Glacial Maximum, Paleoceanography, 28, 726–739,, 2013. 

Wagner, M. and Hendy, I. L.: Trace metal evidence for a poorly ventilated glacial Southern Ocean, Clim. Past Discuss., 11, 637–670,, 2015. 

Wallmann, K.: Is late Quaternary climate change governed by self-sustained oscillations in atmospheric CO2?, Geochim. Cosmochim. Ac., 132, 413–439,, 2014. 

Wallmann, A.: Effects of Eustatic Sea-Level Change on Atmospheric CO2 and Glacial Climate, Nova Act. LC NF, 121, 241–245, 2015. 

Watson, A. J. and Naveira Garabato, A. C.: The role of Southern Ocean mixing and upwelling in glacial-interglacial atmospheric CO2 change, Tellus B, 58, 73–87,, 2006. 

Weber, S. L., Drijfhout, S. S., Abe-Ouchi, A., Crucifix, M., Eby, M., Ganopolski, A., Murakami, S., Otto-Bliesner, B., and Peltier, W. R.: The modern and glacial overturning circulation in the Atlantic ocean in PMIP coupled model simulations, Clim. Past, 3, 51–64,, 2007. 

Weitemeyer, K. A. and Buffett, B. A.: Accumulation and release of methane from clathrates below the Laurentide and Cordilleran ice sheets, Global Planet. Change, 53, 176–187, 2006. 

Williamson, M. S., Lenton, T. M., Shepherd, J. G., and Edwards, N. R.: An efficient numerical terrestrial scheme (ENTS) for Earth system modelling, Ecol. Model., 198, 362–374,, 2006. 

Yu, J., Anderson, R., and Rohling, E.: Deep Ocean Carbonate Chemistry and Glacial-Interglacial Atmospheric CO2 Change, Oceanography, 27, 16–25,, 2014. 

Zech, R., Huang, Y., Zech, M., Tarozo, R., and Zech, W.: High carbon sequestration in Siberian permafrost loess-paleosols during glacials, Clim. Past, 7, 501–509,, 2011. 

Zeng, N.: Glacial-interglacial atmospheric CO2 changes – the Glacial Burial Hypothesis, Adv. Atmos. Sci., 20, 677–693, 2003. 

Zeng, N.: Quasi-100 ky glacial-interglacial cycles triggered by subglacial burial carbon release, Clim. Past, 3, 135–153,, 2007. 

Zhang, X., Lohmann, G., Knorr, G., and Xu, X.: Different ocean states and transient characteristics in Last Glacial Maximum simulations and implications for deglaciation, Clim. Past, 9, 2319–2333,, 2013. 

Ziegler, M., Diz, P., Hall, I. R., and Zahn, R.: Millennial-scale changes in atmospheric CO2 levels linked to the Southern Ocean carbon isotope gradient and dust flux, Nat. Geosci., 6, 457–461, 2013. 

Zimov, S. A., Schuur, E. A. G., and Chapin, F. S., III.: Permafrost and the global carbon budget, Science, 312, 1612–1613,, 2006. 

Zimov, N. S., Zimov, S. A., Zimova, A. E., Zimova, G. M., Chuprynin, V. I., and Chapin, F. S.: Carbon storage in permafrost and soils of the mammoth tundra-steppe biome: Role in the global carbon budget, Geophys. Res. Lett., 36, L02502,, 2009. 


The ΔLithC stems from changes in the deep-sea CaCO3 burial flux and/or CaCO3 weathering/shallow water deposition flux and was initially calculated to ensure that carbon was being conserved over the LGM simulation.

Short summary
We simulate the Last Glacial Maximum atmospheric CO2 decrease with a large ensemble of parameter sets to investigate the range of possible physical and biogeochemical Earth system changes accompanying the CO2 decrease. Amongst the dominant ensemble changes is an increase in terrestrial carbon, which we attribute to a slower soil respiration rate, and the preservation of carbon by the LGM ice sheets. Further investigation into the role of terrestrial carbon is warranted.