Simulation of climate , ice sheets and CO 2 evolution during the last four glacial cycles with an Earth system model of intermediate complexity

In spite of significant progress in paleoclimate reconstructions and modeling of different aspects of the past glacial cycles, 10 the mechanisms which transform regional and seasonal variations in solar insolation into long-term and global-scale glacialinterglacial cycles are still not fully understood, in particular, for CO2 variability. Here using the Earth system model of intermediate complexity CLIMBER-2 we performed simulations of co-evolution of climate, ice sheets and carbon cycle over the last 400,000 years using the orbital forcing as the only external forcing. The model simulates temporal dynamics of CO2, global ice volume and other climate system characteristics in good agreement with paleoclimate reconstructions. These 15 results provide strong support to the idea that long and strongly asymmetric glacial cycles of the late Quaternary represent a direct but strongly nonlinear response of the Northern Hemisphere ice sheets to the orbital forcing. This direct response is strongly amplified and globalized by the carbon cycle feedbacks. Using simulations performed with the model in different configurations, we also analyze the role of individual processes and sensitivity to the choice of model parameters. While many features of simulated glacial cycles are rather robust, some details of CO2 evolution, especially during glacial 20 terminations, are sensitive to the choice of model parameters. Specifically, we found two major regimes of CO2 changes during terminations: in the first one, when the recovery of the Atlantic meridional overturning circulation (AMOC) occurs only at the end of the termination, a pronounced overshoot in CO2 concentration occurs at the beginning of the interglacial and CO2 remains almost constant during interglacial or even decline towards the end, resembling Eemian CO2 dynamics. However, if the recovery of the AMOC occurs in the middle of the glacial termination, CO2 concentration continues to rise 25 during interglacial, similar to Holocene. We also discuss potential contribution of the brine rejection mechanism for the CO2 and carbon isotopes in the atmosphere and the ocean during the past glacial termination.


Introduction
Antarctic ice cores reveal that during the past 800 kyr, the atmospheric CO 2 concentration (Petit et al., 1999;Jouzel et al., 2007) varied synchronously with the global ice volume (Waelbroeck et al., 2002;Spratt and Lisiecki, 2016).The most straightforward explanation for this fact is that CO 2 drives glacial cycles together with orbital variations, and the longest, 100 kyr component of the late Quaternary glacial cycles, which is absent in the orbital forcing, is the direct response to CO 2 forcing where the 100 kyr component is the dominant one.However, simulations with climate-icesheet models of different complexity (e.g.Berger et al., 1999;Crowley and Hyde, 2008;Ganopolski and Calov, 2011;Abe-Ouchi et al., 2013) show that long glacial cycles (i.e.cycles with typical periodicity of ca. 100 kyr) can be simulated even with a constant CO 2 concentration if the latter is sufficiently low.Moreover, these model simulations show that not only the dominant periodicity but also the timing of glacial cycles, can be correctly simulated without variable CO 2 forc-ing.This fact strongly suggests an opposite interpretation of close correlation between global ice volume and CO 2 during Quaternary glacial cycles -namely that glacial cycles represent a strongly nonlinear response of the Earth system to the orbital forcing (Paillard, 1998), while variations in CO 2 concentration are directly driven by ice sheet fluctuations.In turn, CO 2 variations additionally strongly amplify and globalise the direct response of ice sheets to the orbital forcing.
In spite of the significant number of studies aimed at explaining low glacial CO 2 concentrations (e.g.Archer et al., 2000;Sigman and Boyle, 2000;Watson et al., 2000), the influence of ice sheets on the carbon cycle remains poorly understood.It is also unclear the extent to which CO 2 variations represent a direct response to ice sheet forcing and how much is the result of additional amplification of CO 2 variations through the climate-carbon cycle feedback.Indeed, although radiative forcing of ice sheets contributes about 50 % to glacial-interglacial variations in global temperature (Brady et al., 2013), most of the cooling associated with ice sheets is restricted to the area covered by ice sheets and their close proximity.Thus, the direct contribution of ice sheets to glacial ocean cooling is rather limited, and therefore the effect of ice sheets on CO 2 drawdown through the solubility effect can explain only a fraction of the reduction in glacial CO 2 .At the same time, the direct effect of ice sheets on atmospheric CO 2 concentration through ca. 3 % changes in the ocean volume and global salinity is rather well understood but works in the opposite direction and leads to a glacial CO 2 rise of about 10-20 ppm (Sigman and Boyle, 2000;Brovkin et al., 2007).Another direct effect of ice sheet growth on the carbon cycle through reducing the area covered by forest (e.g.Prentice et al., 2011) also operates in the opposite direction.However, several other processes could potentially contribute to glacial CO 2 drawdown through ice sheet growth and related lowering of sea level.One such mechanism is enhanced biological productivity in the Southern Ocean due to the iron fertilisation effect (Martin, 1990;Watson et al., 2000).The latter is attributed to enhanced dust deposition over the Southern Ocean seen in the paleoclimate records (Martinez-Garcia et al., 2014;Wolff et al., 2006).At least part of this enhanced deposition is associated with the dust mobilisation from exposed Patagonian shelf and glaciogenic dust production related to the Patagonian ice cap (Mahowald et al., 1999;Sugden et al., 2009).A number of studies on the effect of iron fertilisation suggested a contribution of 10 to 30 ppm to the glacial CO 2 decrease (e.g. Watson et al., 2000;Brovkin et al., 2007).Another effect is related to the brine rejection mechanism, more specifically, to a much deeper penetration of brines produced during sea ice formation in the Southern Ocean during glacial time.The latter is explained by shallowing and significant reduction of the Antarctic Shelf area.According to Bouttes et al. (2010) this mechanism, in combination with enhanced stratification of the deep ocean, can contribute up to 40 ppm to the glacial CO 2 lowering.
Apart from the mechanisms mentioned above, many other processes have been proposed to explain low glacial CO 2 concentration.Among them are changes in the ocean circulation (Watson et al., 2015) and an increase in Southern Ocean stratification (e.g.Kobayashi et al., 2015), increase in sea ice area in the Southern Ocean (Stephens and Keeling, 2000) and a shift in the westerlies (Toggweiler et al., 2006), increase in nutrient inventory or change in the marine biota stoichiometry (Sigman and Boyle, 2000;Wallmann et al., 2016), changes in coral reefs accumulation and dissolution (Opdyke and Walker, 1992), accumulation of carbon in the permafrost regions (Ciais et al., 2012;Brovkin et al., 2016), variable volcanic outgassing (Huybers and Langmuir, 2009), and several other mechanisms.Most of these processes are not directly related to ice sheet area or volume, and thus should be considered as amplifiers or modifiers of the direct response of CO 2 to ice sheets operating through the climate-carbon cycle feedbacks.Although paleoclimate records provide some useful constraints, the relative role of particular mechanisms at different stages of glacial cycles remains poorly understood.
Most studies of glacial-interglacial CO 2 variations performed to date have been aimed at explaining the low CO 2 concentration at the Last Glacial Maximum (LGM, ca.21 ka).In these studies, both continental ice sheets and the radiative forcing of low glacial CO 2 concentration were prescribed from paleoclimate reconstructions.Only a few have attempted to explain CO 2 dynamics during a part of (usually the glacial termination) or the entire last glacial cycle, with models of varying complexity -from simple box-type models (e.g.Köhler et al., 2010;Wallmann et al., 2016), to models of intermediate complexity (Brovkin et al., 2012;Menviel et al., 2012), or a stand-alone complex ocean carbon cycle model (Heinze et al., 2016).In all these studies, radiative forcing of CO 2 (or total GHGs) was prescribed based on paleoclimate reconstructions.Similarly, ice sheets' distribution and elevation were derived from paleoclimate reconstructions or model simulations where radiative forcing of GHGs has been prescribed.Thus, in all these studies, CO 2 was treated as an external forcing rather than an internal feedback.Here we for the first time perform simulations of the Earth system dynamics during the past four glacial cycles using fully interactive ice sheet and carbon cycle modelling components, and therefore the only prescribed forcing in this experiment is the orbital forcing.

CLIMBER-2 model description
In this study we used the Earth system model of intermediate complexity CLIMBER-2 (Petoukhov et al., 2000;Ganopolski et al., 2001).CLIMBER-2 includes a 2.5-dimensional statistical-dynamical atmosphere model, a three-basin zonally averaged ocean model coupled to a thermodynamic sea ice model, the three-dimensional thermomechanical ice sheet model SICOPOLIS (Greve, 1997), the dynamic model of the terrestrial vegetation VECODE (Brovkin et al., 1997), and the global carbon cycle model (Brovkin et al., 2002(Brovkin et al., , 2007)).Atmosphere and ice sheets are coupled bidirectionally using a physically based energy balance approach (Calov et al., 2005).The ice sheet model is only applied to the Northern Hemisphere.The contribution of the Antarctic ice sheet to global ice volume change is assumed to be constant during glacial cycles and equal to 10 %.The model also includes parameterisation of the impact of aeolian dust deposition on snow albedo (Calov et al., 2005;Ganopolski et al., 2010).The CLIMBER-2 model in different configurations has been used for numerous studies of past and future climates -in particular, simulations of glacial cycles (Ganopolski et al., 2010(Ganopolski et al., , 2016;;Ganopolski and Calov, 2011;Willeit et al., 2015) and carbon cycle operation during the last glacial cycle (Brovkin et al., 2012).
As has been shown by Ganopolski and Roche (2009), the temporal dynamics of the Atlantic meridional overturning circulation (AMOC) during glacial terminations in CLIMBER-2 are very sensitive to the magnitude of freshwater flux to the North Atlantic.To explore different possible deglaciation evolutions, together with the standard model version, we performed an additional suite of simulations, in which the component of freshwater flux into the ocean originated from melting of ice sheets was uniformly scaled up or down by up to 10 %.This rather small change in the freshwater forcing (typically smaller than 0.02 Sv) does not affect AMOC dynamics appreciably most of the time but does induce a strong impact during deglaciations (see below).Other modifications of the climate-ice-sheet component of the model are described in the Appendix.
The ocean carbon cycle model includes modules for marine biota, oceanic biogeochemistry, and deep ocean sediments.Biological processes in the euphotic zone (the upper 100 m in the model) are explicitly resolved using the model for plankton dynamics by Six and Maier-Reimer (1996).The sediment diagenesis model (Archer, 1996;Brovkin et al., 2007) calculates burial of CaCO 3 in the deep sea, while shallow-water CaCO 3 sedimentation is simulated based on the coral reef model (Kleypas, 1997) driven by sea level change.Silicate and carbonate weathering rates are scaled to the runoff from the land surface; they are also affected by sea level change (Munhoven, 2002).Compared to Brovkin et al. (2012), the carbon cycle model has been modified in several aspects.Similar to Brovkin et al. (2012), the efficiency of nutrient utilisation in the Southern Ocean is set to be proportional to the dust deposition rate (see Appendix), which in the case of one-way coupling is prescribed to be proportional to the dust deposition in the EPICA ice core.However, in the fully interactive experiment, the dust deposition rate over the Southern Ocean has been computed from simulated sea level (see Appendix).This means that in the fully interactive experiments (see below) we did not explicitly use any paleoclimate data to drive the model, and the orbital forcing was the only prescribed forcing.In the marine carbon cycle component, we also account for a dependence of the remineralisation depth on ocean temperature following Segschneider and Bendtsen (2013) (see Appendix).In our previous studies the remineralisation depth was kept constant.
The CLIMBER-2 model used in earlier studies of glacial carbon cycle did not include long-term terrestrial carbon pools such as permafrost carbon, peat, and carbon buried beneath the ice sheets.In the present version of the model, these pools are included.The model also accounts for peat accumulation.Modification of the terrestrial carbon cycle components is described in detail in the Appendix.For simulation of atmospheric radiocarbon during the last glacial termination we used the rate of 14 C production following the scenario of Hain et al. (2014), which is based on the production model of Kovaltsov et al. (2012).

One-way coupled and fully interactive experiments
In our previous experiments performed with the CLIMBER-2 model (Brovkin et al., 2012;Ganopolski and Brovkin, 2015) we prescribed not only variations in the Earth's orbital parameters (eccentricity, precession, and obliquity) but also the radiative effect of GHGs (CO 2 , CH 4 , and N 2 O) computed using their concentrations from the ice cores records (Lüthi et al., 2008;Petit et al., 1999).In these experiments, which we denote hereafter as "one-way coupled" (Fig. 1a, Table 1), atmospheric CO 2 concentration was computed by the carbon cycle model but not used as the radiative forcing for the climate component.Similarly, in these experiments, the CO 2 fertilisation effect on vegetation was computed using reconstructed CO 2 concentration.Therefore, in one-way coupled experiments there were no feedbacks of the simulated atmospheric CO 2 concentration on climate.In the present study, we performed a suite of one-way coupled experiments for the last four glacial cycles but we also performed simulations in which the orbital forcing was the only prescribed external forcing.Since CLIMBER-2 does not include methane and N 2 O cycles and does not account for these GHGs in its radiative scheme, we made use of the fact that CO 2 is the dominant GHG and that temporal variations of the other two closely follow CO 2 .To account for the effect of methane and N 2 O forcings, we computed the effective CO 2 concentration used in the radiative scheme of the model in such a way that radiative forcing of effective CO 2 exceeds radiative forcing of simulated CO 2 by 30 % at any time.This type of experiment we refer to as "fully interactive" (Fig. 1b).In the fully interactive experiment we also use computed CO 2 concentration in the terrestrial component to account for the CO 2 fertilisation effect.As stated above, the dust deposition rate over the Southern Ocean which is used in the parameterisation of the iron fertilisation effect was computed from the global sea level.The radiative forcing of aeolian dust and dust deposition on ice sheets (apart from the glaciogenic dust sources) in both types of experiments were obtained identically to Calov et al. (2005) and Ganopolski and Calov (2011) by scaling the fields computed with global climate models, where the scaling parameter was proportional to the global ice volume.

Model spin-up
The model spin-up and the proper choice of model parameters for simulation of multiple glacial cycles represent a challenge when using the model with long-term components of the carbon cycle because inconsistent initial conditions or even a small imbalance in carbon fluxes could lead to a large drift in simulated atmospheric CO 2 concentration (in the case of one-way coupling) or the state of the entire Earth system (climate, ice sheets, CO 2 ) in the case of the fully interactive experiments.Note that, in the latter case, the negative climate-weathering feedback will eventually stabilise the system but this occurs at a timescale of several glacial cycles and over this time climate could drift far away from its realistic state.To avoid such drift, volcanic outgassing should be carefully calibrated.Based on a set of sensitivity experiments, we found that the value of 5.3 Tmol C yr −1 allows us to simulate quasiperiodic cycles without a long-term trend in atmospheric CO 2 .Note that even a ±10 % change in volcanic outgassing leads to a significant (order of 100 ppm) drift in CO 2 concentration simulated over the last four glacial cycles.
When the carbon cycle model incorporates such long-term processes as terrestrial weathering, marine sediment accumulation, and permafrost carbon burial, the assumption that the system is close to equilibrium at the pre-industrial period or at any other moment of time is not valid even if the CO 2 concentration was relatively stable during a certain time interval.To produce proper initial conditions at 410 ka we performed a sequence of 410 kyr long one-way coupled runs with the identical forcings.In the first run we used as the initial conditions the final state obtained in simulation of the last glacial cycles (Brovkin et al., 2012).Then we launched each 410 kyr experiment from the final state obtained in the previous model run.The results of such sequence of experiments reveal a clear tendency to converge to the solution with similar initial and final states of the Earth system.We then used the state of climate and carbon cycle obtained at the end of the last run as the initial conditions for all experiments presented in this paper.In the analysis of all experiments described below, we exclude the first 10 000 years.

Simulations of the last four glacial cycles
Realistic simulation of climate and carbon cycle evolution during the last four glacial cycles is more challenging in the case of the fully interactive configuration, because in this case a number of additional positive feedbacks tend to amplify initial model biases.Therefore we begin our analysis with the one-way coupled simulations similar to that performed in Brovkin et al. (2012).This configuration was also used for the calibration of new parameterisations (see Sect. 4) and a set of sensitivity experiments for the last glacial termination (Sect.5).

Experiments with one-way coupled climate-carbon cycle model
Simulated climate and ice sheet variations in the one-way coupled experiments are similar to the ones presented in Ganopolski and Calov (2011), which is not surprising since the only difference between model versions used in these studies is related to the coupling between ice sheet and climate components (see Appendix).Simulated glacial cycles are characterised by global surface air temperature variations of about 5 • C (not shown) and maximum sea level drops by more than 100 m during several glacial maxima.Simulated global ice sheet volume most of the time is close to the reconstructed one (Spratt and Lisiecki, 2016) (Fig. 2d).In general, differences between simulated and reconstructed global sea level are comparable to the uncertainties in sea level reconstructions obtained using different methods.Simulated CO 2 concentration (Fig. 2e) is also in a good agreement with reconstructions based on several Antarctica ice cores (Barnola et al., 1987;Monnin et al., 2004;Petit et al., 1999;Lüthi et al., 2008).The model correctly reproduces the magnitude of glacial-interglacial CO 2 variability of about 80 ppm.Results of simulations with the standard  Spratt and Lisiecki (2016), and in (e) is the compiled Antarctic CO 2 record from Lüthi et al. (2008).Radiative forcing of GHGs in (b) is from Ganopolski and Calov (2011).Antarctic dust is from Augustin et al. (2004).Blue lines in (d, e) correspond to the baseline experiment ONE_1.0 and purple lines to experiment ONE_1.1,where meltwater flux into the Atlantic was scaled up by a factor of 1.1.model version (ONE_1.0)and model with 10 % enhanced meltwater flux (ONE_1.1)are essentially identical most of the time, except for glacial terminations.During glacial terminations even small differences in the freshwater forcing cause pronounced differences in the temporal evolution of the AMOC, and as a result, of CO 2 concentration.As seen in Fig. 2d, in the experiment ONE_1.0,CO 2 concentration grows monotonically during the last glacial termination (TI, midpoint at ca. 15 ka) and Termination IV (ca.330 ka) while it rises faster and overshoots the interglacial level during Terminations II (ca.135 ka) and III (ca.240 ka).In contrast, in experiment ONE_1.1, similar overshoots occur during Terminations I and III but not Termination IV.In all cases, sim- ulated CO 2 lags behind the reconstructed one but this lag is smaller in the cases when overshoot is simulated.Experiments with CO 2 overshoots are clearly in better agreement with empirical data for MIS7 and MIS9.Analysis of model results shows that pronounced CO 2 overshoot occurs in the case when the AMOC is suppressed during the entire glacial termination and recovers only after the cessation of meltwater flux (Fig. 3).In contrast, if the AMOC recovers well before the end of deglaciation, simulated CO 2 experiences only local overshoot and continues to rise during most of the interglacial.The latter behaviour is similar to that was observed in reality during MIS11 and the Holocene, while the former is typical for MIS5, 7, and 9. Thus our model is able to reproduce both types of CO 2 dynamics during recent interglacials.
The rise of CO 2 by 10-20 ppm on the millennial timescale during AMOC shutdowns is a persistent feature of CLIMBER-2 and the mechanism of this rise has been explained in Brovkin et al. (2012) by a weakening of the reverse cell of the Indo-Pacific overturning circulation during periods of reduced AMOC.A similar rise in atmospheric CO 2 concentration during periods of AMOC shutdown has been simulated in some other (but not all) similar modelling experiments.Incorporation of the temperature-dependent remineralisation depth additionally contributes to the CO 2 overshoots at the beginning of several interglacials (see below) but the mechanism described in Brovkin et al. (2012) remains the dominant one.
Comparison of simulated deep ocean δ 13 C with paleoclimate reconstructions (Fig. 4) show that the model correctly simulates larger δ 13 C variability in the deep Atlantic in comparison to the deep Pacific but underestimates the amplitude of glacial-interglacial δ 13 C variability.Simulated atmospheric δ 13 CO 2 shows a rather complex behaviour with an amplitude of variability of 0.6 ‰.The agreement between www.clim-past.net/13/1695/2017/Clim.Past, 13, 1695-1716, 2017 simulated and reconstructed (Eggleston et al., 2016) atmospheric δ 13 CO 2 is rather poor.Both model and data show a drop in atmospheric δ 13 CO 2 during the last and penultimate deglaciations but the data also suggest the strong drop at the end of Eemian interglacial while the model simulated a continuous rise of δ 13 CO 2 during this interval.In addition, temporal variability of the reconstructed δ 13 CO 2 is significantly larger than the simulated one.A more detailed comparison with empirical data during the last deglaciation is presented in Sect. 5. Changes in the ocean oxygenation is considered to be an important indicator of respired carbon storage in the deep ocean, and therefore the proxy for the strength of ocean biological pump.Jaccard et al. (2016) inferred a significant decline in the deep Southern Ocean oxygenation and interpreted it as the result of a combined effect of iron fertilisation by dust and decreased deep ocean ventilation.Our results (Fig. 5) are fully consistent with such an interpretation.The model simulates significant reduction in the dissolved oxygen in the deep Southern Ocean during glacial period.Roughly two-thirds of this reduction is already simulated in the experiment without iron fertilisation and can be solely attributed to reduced deep ocean ventilation.It is noteworthy that changes in the oxygen concentration in this experiment are strongly anticorrelated with the area of sea ice in the Southern Hemisphere (Fig. 5c).This explained by the fact that sea ice directly and indirectly (through stratification of the upper ocean layer) affects gas exchanges between the ocean and the atmosphere.Oxygen concentration is additionally reduced during periods with high dust deposition rate in the experiment, which accounts for the iron fertilisation effect (Fig. 5d).

Experiments with the fully interactive model
In the fully interactive experiments, orbital forcing is the only prescribed forcing and the model does not use any time-dependent paleoclimatological information (such as the Antarctic dust deposition rate used in the one-way coupled experiment).Results of the fully interactive experiment IN-TER_1.0are shown in Fig. 6.For the first experiment of such type ever, the agreement between model simulations and empirical reconstructions is reasonably good.The model simulates correct magnitude and timing of the last four glacial cycles in respect to both sea level and CO 2 concentration.It also reproduces the strong asymmetry of glacial cycles.Naturally, the mismatch between the simulated and reconstructed characteristic in the fully interactive experiments is larger than in the one-way coupled experiments.In particular, in the fully interactive experiment, simulated ice volume is underestimated by 10-20 m in sea level equivalent compared to the reconstructed one.Although the magnitude of glacial-interglacial CO 2 variability in the fully interactive experiment INTER_1.0 is similar to that in the one-way coupled experiment ONE_1.0 and in the reconstructions, the lag between simulated and reconstructed CO 2 during glacial terminations increases additionally in comparison to the oneway coupled experiments.Interestingly, the last glacial cycle and the first 150 kyr of the INTER_1.0 and ONE_1.0 experiments are in very good agreement, while during the time interval between 300 and 150 ka discrepancies are larger.This period corresponds to higher eccentricity and therefore to the larger magnitude of orbital forcing.Similarly to the results of one-way coupled experiments, the fully interactive experiments also show strong sensitivity to magnitude of freshwater flux during glacial terminations.
Comparison of simulated ice sheet spatial distribution and elevation (Fig. 7) shows that the results of one-way coupled (ONE_1.0,Fig. 6a) and the fully interactive experiments (IN-TER_1.0,Fig. 7b) are almost identical during the LGM (the same is true for the previous glacial maxima, not shown) and in reasonable agreement with the reconstructions.Dur-ing glacial terminations, the difference between two experiments increases since in the fully interactive experiment the radiative forcing of GHGs lags considerably behind the reconstructed one used in the one-way coupled experiment.As a result, at 7 ka continental ice sheets melted completely in the one-way coupled experiment (Fig. 7c) while in the fully interactive experiment, a relatively large ice sheet is still present in northeastern Canada (Fig. 7d).
It is instructive to compare frequency spectra of simulated and reconstructed global ice volume in the one-way and fully interactive experiments (Fig. 8).In addition, Fig. 8 shows results from the experiment ONE_240 performed with constant radiative forcing of GHGs corresponding to equivalent CO 2 concentration of 240 ppm.As already shown by Ganopolski and Calov (2011), even with constant CO 2 , the model computes pronounced glacial cycle with 100 kyr periodicity, although it has much weaker amplitude than the reconstructed ice volume.Both model experiments with varying CO 2 radiative forcing (ONE_1.0 and INTER_1.0)reveal much stronger 100 kyr periodicity, which has only a slightly weaker amplitude than the spectrum of the reconstructed ice volume.Interestingly, frequency spectra of ice volume simulated in the one-way and fully interactive experiments have similar powers in 100 kyr and obliquity (40 kyr) bands, but in the precessional band (ca.20 kyr) the one-way coupled experiment reveals a much higher spectral power.This cannot be explained by the prescribed radiative forcing of GHGs because the latter contain very little precessional variability.The stronger precessional component in the ONE_1.0experiment is explained by the fact that in the one-way coupled experiment, the model simulates faster ice sheet growth during the initial part of each glacial cycle and the modelled ice volume variability at the precessional frequency is very sensitive to the global ice volume.

The composition of "the carbon stew" and factor analysis
In this section, we discuss the contribution of different factors to simulated variations in CO 2 concentration.Because none of the mechanisms could explain the CO 2 dynamics in isolation from the other factors (e.g.Sigman and Boyle, 2000; Archer et al., 2000), we call the composition and timing of the mechanisms leading to the glacial CO 2 cycle the "carbon stew".As has been shown in Brovkin et al. (2012), the role of different mechanisms controlling CO 2 concentration at different phases of glacial cycles is different.However, even if we consider only the LGM (as most previous works did), the composition of the "carbon stew" remains highly uncertain, even though there is a growing awareness that both physical and biological processes must have played a comparably important role in glacial CO 2 drawdown (e.g.Schmittner and Somes, 2016;Galbraith and Jaccard, 2015).Obviously, the choice of the "carbon stew" is crucially important for suc- Blue lines represent Ice-5g reconstruction at the LGM (Peltier, 2004).
cessful simulations of glacial cycles.The aim of our paper is not to present the ultimate solution for the "carbon stew" problem since at present this is impossible.Rather we want to demonstrate that with a reasonable representation of physical, geochemical, and biological processes in the model, it is possible to reproduce the main features of Earth system dynamics over the past 400 kyr, including the magnitude and timing of climate, ice volume, and CO 2 variations.Similar to the study by Brovkin et al. (2012), we performed a set of experiments using the one-way coupling technique (see Table 1 for details).We use this approach instead of fully interactive coupling to exclude complex and strongly nonlinear interactions associated with the ice sheet dynamics, which significantly complicates the factor analysis.In the case of the one-way coupled experiments, climate, ice sheets, and other external factors are identical, and these experiments only differ from each other by the parameters of the carbon cycle model.Since CO 2 simulated in the oneway coupled experiment with 10 % enhanced meltwater flux (ONE_1.1) is in a slightly better agreement with observational data than the standard one (ONE_1.0),for the factor analysis we used the experiment ONE_1.1 as the reference one and performed all sensitivity experiments with 10 % enhanced meltwater flux.

The standard carbon cycle model setup
We begin our analysis with the experiment that incorporates only the standard ocean biogeochemistry as described in Brovkin et al. (2007) (Fig. 9).This experiment does not include the effect of the terrestrial carbon cycle.In this configuration, the model is able to explain only about 45 ppm of CO 2 reduction during glacial cycles.Note that this experiment accounts for glacial-interglacial changes in ocean volume of ca. 3 % and corresponding changes in the total biogeochemical inventories including salinity.These volume changes are often neglected in simulations with 3-dimensional ocean models (e.g.Heinze et al., 2016), although these changes have substantial impact on the global carbon cycle, and in our simulations they counteract glacial CO 2 drawdown by ca. 12 ppm.Without the effect of ocean volume reduction, the combination of physical processes and carbonate chemistry can explain up to 57 ppm at the LGM and on average 38 ppm during the entire 400 kyr time interval (see Table 2).This is consistent with the recent results by Buchanan et al. (2016) and Kobayashi et al. (2015).Note that simulated changes in silicate weathering and its impact on atmospheric CO 2 are small, as has already been shown in Brovkin et al. (2012).
Accounting for the land carbon changes does not help to explain the CO 2 concentration changes, since terrestrial carbon contains ca.350 Gt less carbon at the LGM compared to the pre-industrial state.This reduces the glacial-interglacial CO 2 difference by 10-15 ppm compared to the ocean-only experiment (Fig. 9b).Enabling the parameterisation for the iron fertilisation effect in the Southern Ocean results in additional glacial CO 2 drawdown of up to 30 ppm (22 ppm at the LGM), mostly towards the end of each glacial cycle (Fig. 9c).This value is close to that reported by Lambert et al. (2015).With all these processes considered in our previous study (Brovkin et al., 2012), we are still short of ca. 25 ppm to explain the full magnitude of glacial-interglacial variability.

Additional processes included in the carbon cycle model
There are a number of other proposed mechanisms which potentially can explain several tens of ppm of glacial CO 2 decline.Our choice of two processes to obtain the observed magnitude of glacial-interglacial CO 2 variations is somewhat subjective.The chosen mechanisms are explained below, while an alternative one (brine rejection) is discussed in Sect.4.3.
The first additional mechanism added to the mechanisms already described in Brovkin et al. (2012) is temperature-dependent remineralisation depth.In the standard CLIMBER-2 version, remineralisation depth is spatially and temporally constant.Since in the colder ocean remineralisation depth increases, this enhances the efficiency of the carbon pump and contributes to a decrease in atmospheric CO 2 concentration (e.g.Heinze et al., 2016;Menviel et al., 2012;Matsumoto, 2007).Details of the mechanism implementation are described in the Appendix.As seen in Fig. 9e, making remineralisation depth temperature dependent introduces additional glacial-interglacial variability of CO 2 with a magnitude of about 20 ppm.Roughly half of this value is clearly attributed to the CO 2 overshoots seen at the beginning of some interglacials.The reason is that the AMOC shutdowns due to meltwater flux that happened during glacial terminations lead not only to surface cooling in the North Atlantic but also to significant thermocline warming that occurs over the entire Atlantic Ocean (e.g.Mignot et al., 2007).This subsurface warming causes pronounced shoaling of the remineralisation depth and the release of carbon from the ocean into the atmosphere.This process reverses after recovery of the AMOC at the beginning of interglacials.Burley and Katz (2015) and Huybers and Langmuir (2009) proposed that the rate of volcanic outgassing varies during glacial cycle due to variable load of the ice sheet and ocean on the Earth crust.Therefore we assume that volcanic outgassing has a variable component (about 30 % of its averaged value of 5.3 Tmol yr −1 ), which represents the delayed response to the change in global ice volume.This simple parameterisation explained in the Appendix does not affect cumulative volcanic outgassing over the glacial cycle, but contributes to glacial-interglacial CO 2 variability by an additional 10 ppm (Fig. 9d).With varying volcanic out- gassing and temperature-dependent remineralisation depth, CLIMBER-2 reproduces glacial-interglacial CO 2 cycles in a good agreement with paleoclimate records (Fig. 9a).

Brine rejection mechanism
Using While we believe that the brine rejection mechanism belongs to a class of plausible mechanisms contributing to glacial CO 2 drawdown, we did not use brine parameterisation in our simulations for several reasons.Firstly, the parameterisation for brine rejection cannot be tested against observational data.For present-day climate conditions, brine rejection efficiency should be below 0.1, otherwise modern Antarctic bottom water becomes saltier than the North Atlantic deep water, which is at odds with reality.This means that to be an efficient mechanism for glacial CO 2 drawdown, the brine efficiency should increase under glacial conditions, at least by an order of magnitude.Whether this is physically plausible is not clear.The only paleoclimate constraint on the brine efficiency is the reconstruction of paleosalinity based on the pore water (Adkins et al., 2002), which suggests an increase in the deep water salinity in the Southern Ocean by more than 2 psu during the LGM.Such an increase in salinity is indeed difficult to reproduce without contribution of brines.However, the accuracy of salinity reconstructions based on such a method remains uncertain (Wunsch, 2016).Second, there is a problem with the temporal dynamics of brine rejection ef-ficiency.Mariotti et al. (2016) assumed an abrupt decrease in brine rejection efficiency from 0.7 to 0 in a very short interval between 18 and 16 ka.However, both sea level and the size of the Antarctic ice sheets remained essentially constant during this period and therefore there is no obvious reason for such large variations in the brine rejection efficiency.According to the interpretation of Roberts et al. (2016), brine rejection remained efficient during most of the glacial termination and ceased only after 11 ka, when most of the glacial-interglacial CO 2 rise had already been accomplished.In the view of these uncertainties, we decided not to include parameterisations of the brine rejection mechanism in simulations of glacial cycles.However, for simulations of the last glacial termination discussed below, we analysed the potential effect of brine rejection on radiocarbon and other paleoclimate proxies.

Simulation of climate, CO 2 , and carbon isotopes during the last termination
The last glacial termination provides a wealth of paleoclimate records with a potential to better constrain the mechanisms of glacial CO 2 variability.In this section, we discuss the last glacial termination in more detail.Similarly to the previous section, to exclude nonlinear interaction with ice sheets, we discuss here only one-way coupled experiments.
To reduce computational time, we performed experiments only for the last 130 000 years starting from the Eemian interglacial and using the same initial conditions as in the experiments discussed above.
In the standard ONE_1.0_130Kexperiment, the model simulates climate variability across Termination I rather realistically.In particular, it reproduces the temporal resumption of the AMOC in the middle of the termination resembling the Bølling-Allerød warm event (Fig. 10a).However, the timing of this event in our model is shifted by ca.1000 years compared to the paleoclimate records.Results of our experiments reveal high sensitivity of the timing of the AMOC resumption to the magnitude of freshwater flux.A change of the flux by 10 % in the ONE_1.1_130Kexperiment significantly alters millennial-scale variability during the last glacial termination (Fig. 10).This result suggests that simulated millennial-scale variability during the Termination I is not robust -i.e. it is unlikely that a single model run through the glacial termination would reproduce the right timing or even the right sequence of millennial-scale events.
Although simulated CO 2 concentration at the LGM and pre-industrial state are close to observations, simulated CO 2 appreciably lags behind reconstructed CO 2 during the termination (Fig. 10b).This is primarily related to the fact that simulated CO 2 does not start to grow at ca. 18 ka BP as reconstructed, but only after the end of the simulated analogue of the Bølling-Allerød event.At the same time, in agreement with paleoclimate reconstructions, CO 2 concentration  (Lüthi et al., 2008;Schmitt et al., 2012); (d) IntCal13 radiocarbon calibration curve (Reimer et al., 2013).reaches a local maximum at the end of the North Atlantic cold event, which resembles the Younger Dryas event.Simulated CO 2 concentration also reveals a continuous CO 2 rise during the Holocene towards its pre-industrial value of 280 ppm.This result confirms that such CO 2 dynamics could be explained by natural mechanisms alone and do not require early anthropogenic CO 2 emissions until ca. 2 ka (Kleinen et al., 2016).This result also demonstrates that the temporal dynamics of CO 2 during interglacials critically depend on the timing of final AMOC recovery.Late recovery during glacial termination causes a strong overshoot of CO 2 at the beginning of the interglacial followed by some decrease or a stable CO 2 concentration.However, if the complete AMOC recovery occurs well before the end of termination, only temporal CO 2 overshoot occurs and CO 2 continues to rise during the entire interglacial.
It is generally believed that atmospheric δ 13 C provides useful constraint on the mechanisms of deglacial CO 2 rise (Schmitt et al., 2012;Joos et al., 2004;Fischer et al., 2010).
A. Ganopolski and V. Brovkin: Simulation of the last four glacial cycles Simulated atmospheric δ 13 C drops from the LGM level of about −6.4 ‰ to the minimum value of −6.7 ‰ between 16 and 14 ka (Fig. 10c).This is primarily related to the reduction of marine biological productivity which, in turn, is explained by the decrease in iron fertilisation effect over the Southern Ocean during the first part of Termination I.The magnitude of the δ 13 C drop is in a good agreement with empirical data (Fig. 10c).The model is also able to simulate W-shaped δ 13 C evolution associated with reorganisation of the AMOC.However, this W shape is shifted in time compared to the reconstructed one by ca.1000 years because the model analogue of the Bølling-Allerød event occurs earlier than the real one by the same amount of time.Note that this local maximum in δ 13 C is completely absent in the experiment ONE_1.1_130K,where temporal resumption of the AMOC during glacial termination does not occur.δ 13 C rise after 12 ka is primary attributed to the accumulation of carbon in terrestrial carbon pools (forest regrowth and peat accumulation).At the same time, simulated present-day atmospheric δ 13 C is underestimated compared to ice core data by ca.0.15 ‰.The model simulates an almost monotonic decrease in atmospheric 14 C from the LGM to present.Most of this decrease (ca.200 ‰) is caused by a prescribed production rate which was about 20 % higher during LGM.Only about 80 ‰ of 14 C is attributed to a difference in climate state between LGM and the present, primarily due to a less ventilated deep ocean.As shown in Fig. 10d, simulated atmospheric 14 C is significantly underestimated before 12 ka compared to the reconstruction by Reimer et al. (2013) and at the LGM this difference reaches more than 100 ‰.It is possible but unlikely that such big differences can be attributed to uncertainties in reconstructed production rate.An alternative hypothesis for explaining this mismatch is discussed below.
Figure 11 shows LGM time slice anomalies and temporal evolution of δ 13 C and radiocarbon ventilation age during Termination I in the Atlantic Ocean simulated in experiment ONE_1.0_130K.The spatial distribution of glacial anomalies and temporal dynamics of δ 13 C and radiocarbon ventilation age during termination are qualitatively very similar.Glacial δ 13 C in the deep Atlantic at the LGM is 0.6-1 ‰ lower than at present; this is primarily related to a shoaling of the AMOC and reduced ventilation in the Southern Ocean.The vertical distribution of δ 13 C anomalies at the LGM is consistent with the paleoclimate reconstructions (e.g.Hesse et al., 2011).
Simulated ventilation age at the LGM can be directly compared with Skinner et al. (2017) (their Fig. 4a and c).Both models and data show significant increase in radiocarbon ventilation age in the deep Atlantic.However, the spatial patterns of ventilation age changes are rather different.In the model, the largest increase in the ventilation age occurs in the deep North Atlantic, which is explained by shoaling of the AMOC cell and increased presence of the poorly ventilated Southern Ocean water masses.At the same time, the reconstructed radiocarbon ages in the deep North Atlantic are characterised by very large scattering (from 1000 to 3000 14 C years) and it is unclear whether their average values can be directly compared to the results of the zonally averaged ocean model.During glacial termination, both δ 13 C and radiocarbon ventilation age show pronounced response at all depths to the millennial-scale reorganisations of the AMOC (Fig. 12b  and d).The ventilation age in the deep Atlantic, which is about 2000 years prior to the model analogue of the warm Bølling-Allerød event, rapidly reaches nearly the modern level after the AMOC resumption and drops again to the glacial level during the model analogue of the cold Younger Dryas event.Such evolution of ventilation age in the North Atlantic is in agreement with paleoclimate reconstructions (Robinson et al., 2005;Skinner et al., 2014).

Brine rejection mechanisms and radiocarbon in the ocean and atmosphere
As discussed above, our version of the CLIMBER-2 model is not able to accurately reproduce atmospheric 14 C decline during the first part of the last glacial termination.At the same time, Mariotti et al. (2016) demonstrated that their version of CLIMBER-2, which incorporates the mechanism of brine rejection, is able to simulate larger atmospheric 14 C decrease from the LGM to the present, which is in better agreement with the observational data (Reimer et al., 2013).By introducing a similar parameterisation for brine rejection and stratification-dependent vertical diffusivity in our model, we are able to reproduce the results of Mariotti et al. (2016) (Fig. 12).It is noteworthy that we use different temporal dynamics of the efficiency of brine rejections during termination.Instead of abrupt and non-monotonic changes in the brine efficiency prescribed in Mariotti et al. (2016), in the ONE_BRINE_130K experiment we assume that this efficiency is 0.75 at the LGM, 0 at present, and in between it follows the global ice volume.We do not claim that this scenario is more realistic, but at least it is more consistent with the findings of Roberts et al. (2016).Figure 11 shows that the model with the brine rejection parameterisation and stratification-dependent vertical diffusivity simulates atmospheric 14 C in better agreement with empirical data then the standard version.This is explained by the fact that the brine rejection mechanism in combination with stratificationdependent vertical mixing produces very salty and dense deep water masses which are almost completely insulated from the surface.The comparison of the vertical profiles of ventilation age (Fig. 13) with the basin-averaged data of Skinner et al. (2017)   LGM, which is in good agreement with the 689 ±53 14 C yr reported in Skinner et al. (2017).At the same time, the model with the brine rejection parameterisation simulates an increase in the ventilation age of more than 1300 14 C yr.
Interestingly, the two model versions do not differ much with respect to simulated deep ocean δ 13 C. Finally, the two model versions differ significantly with respect to the deep Southern Ocean salinity.Change in salinity in the standard model version is only about 1 psu, which is close to the global mean salinity change due to ice sheet growth.The model version with the brine rejection parameterisation simulates glacial deep South Ocean salinity of 37 psu, which is in good agreement with the reconstruction by Adkins et al. (2002).Thus we found that including additional effects (brines and stratification-dependent diffusion) helps to bring atmospheric 14 C and the deep Southern Ocean salinity into better agreement with available reconstructions, but at the expense of very old (likely to be at odds with paleoclimate data) water masses in the deep ocean.Of course, these results are obtained with a very simple ocean component and it is possible that more realistic ocean models would be able to resolve this apparent contradiction.

Changes in the terrestrial carbon cycle
The evolution of the carbon cycle in the one-way coupled simulation is presented in Fig. 13.The conventional components of the land carbon cycle (vegetation biomass, soil carbon stored in non-frozen and non-flooded environment) change between 1400 Gt C during the LGM and 2000 Gt C during interglacial peaks.Such an amplitude of 600 Gt C of glacial-interglacial changes is typical for the models of the land carbon cycle without long-term components (Kaplan et al., 2002;Joos et al., 2004;Brovkin et al., 2002).However, when we account for permafrost, peat, and buried carbon, the magnitude decreases to 300-400 Gt C.This is due to the counteracting effect of the permafrost and buried carbon pools relative to the conventional components.Both these pools vary between 0 and 350 Gt C and reach their maxima during glacials.The peat storage also reaches about 350 Gt C, but it grows only during interglacials or warm stadials.Let us note that during glacial inceptions, while biomass and mineral soil carbon decrease, terrestrial carbon storage increases due to an increase in buried and permafrost carbon.As a result, total land carbon did not change much during the period of large ice sheet initiation.
During the last deglaciation (Fig. 14, right panels), the peat storages increase monotonically reaching ca.350 Gt C at the pre-industrial period.The conventional carbon pools increase from 1400 to 1800 Gt C at the peak of interglacial (ca. 9 ka), and then start to decline due to the orbital forcing effect on climate in the Northern Hemisphere.The permafrost and buried carbon pools show the opposite behaviour, experiencing minimum at 10 and 5 ka, respectively, and grow afterwards.The combined effect on the total land carbon is a monotonic increase during interglacials, mostly because of peat accumulation.

Conclusions
We present here the first simulations of the last four glacial cycles with the one-way and two-way coupled carbon cycle model.The model is able to reproduce the major aspects of glacial-interglacial variability of climate, ice sheets, and of atmospheric CO 2 concentration even when driven by the orbital forcing alone.These results provide strong support to the idea that long and strongly asymmetric glacial cycles of the late Quaternary represent a direct but strongly nonlinear response of the Northern Hemisphere ice sheets to the orbital forcing which, in turn, is amplified and globalised by the carbon cycle feedback.
The model simulates correct timing of the past glacial terminations in terms of ice volume while the simulated CO 2 concentration lags behind the reconstructed one by several thousand years.The model is also able to simulate temporal evolution of the stable carbon isotope in the ocean.At the same time, the agreement between simulated and reconstructed atmospheric δ 13 C is rather poor.Similarly, the magnitude of simulated atmospheric 14 C decline during the last glacial termination is underestimated.Introducing the brine rejection parameterisation and stratification-dependent vertical diffusivity allows us to improve the agreement for atmospheric 14 C but leads to unrealistically "old" glacial deep ocean water masses.
The temporal dynamics of CO 2 during interglacial depend strongly on the timing of the AMOC recovery during glacial termination.If the AMOC recovers only at the end of glacial termination, CO 2 concentration experiences an overshoot at the beginning of interglacial and then CO 2 declines.In contrast, early recovery of the AMOC leads to a monotonic rise of CO 2 during interglacials.In our simulations, millennialscale variability during the last glacial termination is very sensitive to the magnitude of meltwater flux, and the sequence and timing of simulated millennial-scale events are not robust even when the model is forced by prescribed radiative forcing of GHGs.
Adding new long-term carbon pools (peat, buried and permafrost carbon) decreases the amplitude of glacialinterglacial changes in the total land carbon storage.It helps to reduce an effect of terrestrial biosphere on the CO 2 change during glacial inception and to a lesser extent during glacial terminations.
This work demonstrates that simulation of glacial cycles with Earth system models driven by orbital forcing alone is possible.This does not mean that we have presented here the ultimate recipes for all processes and feedbacks and, in particular, for "the carbon stew".The understanding of how the global carbon cycle operates on orbital and suborbital timescales still remains incomplete, and large uncertainties remain in the choice of individual processes and their parameterisations.Paleoclimate data provide some useful constraints but the proxy data syntheses are far from being in a perfect state, with some proxies telling contradictory stories and others being difficult to interpret.
The CLIMBER-2 model is a rather simple and coarse resolution Earth system model.This allows us to perform a large ensemble of model simulations on orbital and even longer timescales.Obviously, such a fast model has significant limitations.In particular, it employs the zonally averaged ocean model.Many essential processes, such as iron fertilisation effect, are parameterised.The development of a high-resolution state-of-the-art Earth system model suitable for simulation of the interaction between climate, ice sheets, and carbon cycle at orbital timescales is crucial to make the next step forward in understanding of Earth system dynamics during the Quaternary.

A1 Modifications of terrestrial carbon cycle model
The old version of the CLIMBER-2 carbon cycle model described in Brovkin et al. (2002) considers two vegetation types -trees and grass.Each of the two vegetation types has four carbon pools -leaves, stems, and fast and slow soil carbon.Each of these four pools occupies the same fraction of grid cell as the respective vegetation type.Crichton et al. (2014) in their version of permafrost carbon implementation in CLIMBER-2 have not changed the pool structure but modified turnover time, assuming that it is increasing under permafrost conditions.In the new version of the carbon cycle model which we use in the present work, we have introduced three new carbon pools: boreal peat, permafrost, and carbon buried under ice sheets.The fractions of land covered by grass and trees are computed in the vegetation model following Brovkin et al. (1997), the fraction of land covered by ice sheets is computed by the ice sheet model and the fraction of permafrost f pm for the temperature range −5 • C < T ts < 5 • C is computed in the land surface module as where T ts is annual mean top soil layer temperature.It is assumed that grass (in boreal latitudes this mean tundra) is located north of forest and therefore freezes first.Only if permafrost exceeds grass fraction can the permafrost expand over the area covered with trees.During ice sheet growth, all carbon under ice sheets apart from the living biomass is reallocated into the buried carbon pool.Buried carbon remains intact until it is covered by ice sheets.During deglaciation, this buried carbon is transformed into the permafrost pool.The fraction of land covered by peat is defined as where f * pt is the potential fraction of peat in each grid cell prescribed from modern observational data and f gc is the fraction of land covered by ice sheets.Note that we do not consider peatlands in low latitudes.Although peat and permafrost have certain areal fractions, they are considered to be parts of grid cells covered by vegetation.Net primary production and fluxes between the fast carbon pools (leaves, stems, and fast soil pool) are computed the same way as in Brovkin et al. (2012).The downward flux of carbon from the fast soil is partitioned between slow soil pool and permafrost proportionally to their relative fractions.The rate of peat accumulation is equal to a fixed fraction of net primary production in the respective vegetation type.Evolution of carbon content p i in slow carbon pools is described by the equation where p i = b i f i , the b i is the concentration of carbon in the ith carbon pool (in kgC m −2 ), and f i is the fraction of the ith pool.The value q i represents the difference between local accumulation and decay of carbon in the pool, b i is carbon concentration in the pool by which the ith pool is expanded, and b i = b i if the ith pool is shrinking.For peat, b i = 0.For the permafrost, the situation is more complex, because it can gain/lose carbon from/to slow soil, peat, and buried carbon pools.The source term for the permafrost pool q pm consists of the sum of fluxes from the fast grassland and tree soil pools into the respective slow pools (see Brovkin et al., 2002, for details) minus the decay term, where the decay timescale is set to 20 000 years.Apart from carbon, the terrestrial carbon model also computes carbon isotope ( 13 C and 14 C) contents in all carbon pools.Since carbon isotopes are also computed in the oceanic carbon cycle model, we can compute δ 13 C and δ 14 C in the atmosphere and compare modelling results with available paleoclimate data.

A2 Modifications of the ocean carbon cycle model A2.1 Dust deposition in the Southern Ocean
In the one-way coupled experiments, similar to Brovkin et al. (2012), we used the concentration of aeolian dust in the Antarctic ice cores as the proxy for iron deposition over the Southern Ocean.Such choice is supported by the recent measurements of iron content in Southern Ocean sediment cores (Lamy et al., 2014).In the fully interactive run, the iron flux over the Southern Ocean (D) in arbitrary units is parameterised through the global sea level change as D = 100 dS dt + 10 max(S − 50; 0) + 1.5S, where S is the ice volume expressed in metres of sea level equivalent and time t is in years.This formula gives significant increase in iron flux for the case when sea level drops below 50 m, which is likely related to the fact that Patagonian dust source is very sensitive to the area of exposed shelf and glacial erosion processes.Numerical parameters in this formula were obtained by fitting simulated D to dust concentration in the Antarctic record.This allows us to use the same parameterisation for the iron fertilisation effect in the one-way and fully interactive experiments.To prevent large fluctuations in the iron flux related to fluctuations of the time derivative of S, the dust deposition D computed by this equation has been smoothed by applying a relaxation procedure.Namely, at each time step n, the dust deposition D n is computed as where ε = 0.001, which is approximately equivalent to introducing a 1000-year filter.

A2.2 Dependence of remineralisation depth on temperature
In CLIMBER-2 the vertical profile of carbon below the euphotic zone is given by , where remineralisation depth z r is held constant equal to 100 m.To take into account the dependence of remineralisation rate on ambient temperature, following Segschneider and Bendtsen (2013), we now use the dependence of z r on the thermocline temperature (300 m) T : where T 0 = 9 • C and z ro = 100 m.The value of T 0 was selected in such a way that introducing the temperaturedependent remineralisation depth does not affect atmospheric CO 2 concentration under pre-industrial climate conditions.During glacial times, temperature in the thermocline decreases by 2-3 • C, which causes an increase in z r by 20-30 %.This results in additional CO 2 drawdown by ca. 15 ppm.

A2.3 Parameterisation of the iron fertilisation effect
The rate of dust deposition which is prescribed from the ice cores in the one-way coupled experiment or computed from global ice volume in the fully interactive experiments is considered to be a proxy for iron flux and is used in the parameterisation of iron the fertilisation mechanism.This parameterisation is only applied to the Southern Ocean (south of 40 • S).As described in Brovkin et al. (2002), net primary production of phytoplankton in the model is computed as = c(D)r(T , R) where C p is the phytoplankton concentration, P is the phosphorus concentration in the euphotic zone, r is a function of temperature T and photosynthetic active insolation R, f is the fraction of grid cell covered by sea ice, P 0 is a constant, and c is a function of normalised dust deposition rate D. Note that in the case of prescribed dust deposition rate, D was obtained by multiplying observed dust concentration in mg g −1 units by factor 10 −3 .North of 40 • S, parameter c is set to 1 and south of 40 where c d = 2.During glacial maxima the value of c reaches a value of about 1, which implies that during these periods there is no iron limitation in the Southern Ocean.During interglacials, when D is much smaller than 100, c is close to 0.1.The parameters of this parameterisation have been selected (i) to reproduce present-day nutrient concentration in the Southern Ocean, and (ii) to obtain 20-30 ppm additional CO 2 drawdown during glacial maxima due to the iron fertilisation effect.

A3 Variable volcanic outgassing
Following the idea of Huybers and Langmuir (2009), which has been tested already in Roth and Joos (2012), we introduced a dependence of volcanic CO 2 outgassing O on the rate of sea level change.Namely, we assume that volcanic outgassing linearly depends on sea level derivative with the time delay of about 5000 years: Here O 1 = 5.3 Tmol C yr −1 and O 2 = 50 Tmol C m −1 .With these parameters, volcanic outgassing does not change by more than 30 % during all glacial cycles.Note that over one glacial cycle the average value of O is very close to O 1 .

A4 Calculation of radiocarbon ventilation age
For the direct comparison of model results with the empirical reconstruction of marine radiocarbon age offset (Skinner et al., 2017), we calculate the model analogue of this characteristic using simulated relative concentration of the radiocarbon in the atmosphere c atm = 14 C atm / 12 C atm and in the ocean c ocn = 14 C ocn / 12 C ocn using the formula t = −8033 ln c ocn c atm .

A5 Modifications of the energy and surface mass balance interface
In our previous simulations with CLIMBER-2 we found that if maximum ice sheet volume in the Northern Hemisphere exceeds 100 m, the AMOC remains in the off mode during the entire deglaciation.Although this may be realistic for some recent deglaciations, such long AMOC shutdown prevents simulation of complete deglaciation of North America.This is related to the fact that, due to a very coarse spatial resolution of CLIMBER-2, linear interpolation of surface temperature between neighbouring sectors (American and Atlantic) causes a strong cooling over the eastern part of the Laurentide ice sheet, due to the AMOC shutdown (see for example Arz et al., 2007).In the high-resolution climate models, the effect of AMOC shutdown on North America is rather limited compared to Europe (e.g.Zhang et al., 2014;Swingedouw et al., 2009), which is explained by the predominantly eastward direction of air mass transport.To compensate for this resolution-related problem, we made the magnitude of temperature anomaly correction over eastern North America (see Fig. 2 in Ganopolski et al., 2010) dependent on AMOC strength.Namely, for the AMOC strength below 10 Sv, the amplitude of temperature correction is scaled down by a factor of max /10, where max is the maximum of the meridional overturning stream function (in Sv, 1 Sv = 10 6 m 3 s −1 ) in the Atlantic Ocean.With this parameterisation, during complete shutdown of the AMOC, cooling over eastern North America is compensated for by reducing the temperature correction.Introducing this procedure minimises the impact of the AMOC on Laurentide ice sheet mass balance.As the result, now even prolonged AMOC shutdown does not prevent complete melting of the Laurentide ice sheet during glacial terminations.

Figure 2 .
Figure 2. Transient simulations of the last four glacial cycles forced by orbital variations, observed concentration of well-mixed GHGs, and dust deposition rate (one-way coupled experiments).(a) Maximum summer insolation at 65 • N, W m −2 ; (b) radiative forcing (relative to pre-industrial) of well-mixed GHGs, W m −2 ; (c) Antarctic dust deposition rate in relative units; (d) global ice volume expressed in sea level equivalent (m); (e) atmospheric CO 2 concentration (ppm).Dark red colour in (a)-(c) represents prescribed forcings.Black dashed line in (d) is the sea level stack fromSpratt and Lisiecki (2016), and in (e) is the compiled Antarctic CO 2 record fromLüthi et al. (2008).Radiative forcing of GHGs in (b) is fromGanopolski and Calov (2011).Antarctic dust is fromAugustin et al. (2004).Blue lines in (d, e) correspond to the baseline experiment ONE_1.0 and purple lines to experiment ONE_1.1,where meltwater flux into the Atlantic was scaled up by a factor of 1.1.

Figure 3 .
Figure 3. Temporal evolution of the AMOC, Sv (a), and atmospheric CO 2 concentration, ppm (b) during the last four glacial terminations.Blue lines correspond to experiment ONE_1.0 and purple lines to experiment ONE_1.1,where meltwater flux into the Atlantic was scaled up by a factor of 1.1.

Figure 5 .
Figure 5. (a) Simulated CO 2 concentration (ppm); (b) prescribed Antarctic dust deposition rate in relative units; (c) simulated annual mean sea ice area in the Southern Hemisphere (10 6 km 2 ); (d) simulated oxygen concentration in the deep Southern Ocean in (µmol kg −1 ) in the ONE_1.1 experiment (solid line) and the identical experiment but without iron fertilisation effect (dashed line).

Figure 6 .
Figure 6.Transient simulations of the last four glacial cycles forced by orbital variations only (fully interactive experiments).(a) Maximum summer insolation at 65 • N (W m −2 ); (b) radiative forcing (relative to pre-industrial) of well-mixed GHGs (W m −2 ); (c) Antarctic dust deposition rate in relative units; (d) global ice volume expressed in sea level equivalent (m); (e) atmospheric CO 2 concentration, ppm.Black line in (b) is radiative forcing of GHGs from Ganopolski and Calov (2011).Black dashed line in (c) is Antarctic dust from Augustin et al. (2004), in (d) is sea level stack from Spratt and Lisiecki (2016), and in (e) is compiled Antarctic CO 2 record from Lüthi et al. (2008).Blue lines in (d, e) correspond to the fully interactive experiment INTER_1.0 and purple lines to the experiment INTER_1.1,where meltwater flux into the Atlantic was scaled up by a factor of 1.1.

Figure 8 .
Figure 8. Transient simulations of the last four glacial cycles forced by orbital variations, with prescribed, interactive, and fixed concentrations of well-mixed GHGs.(a) Maximum summer insolation at 65 • N (W m −2 ); (b) temporal evolution of reconstructed and simulated sea level (m); (c) frequency spectra of the global ice volume; (d) frequency spectra of boreal summer insolation.Black line is for the data(Spratt and Lisiecki, 2016), blue line corresponds to the one-way coupled experiment ONE_1.0,red line to the fully interactive experiment INTER_1.0, and green line to the ONE_240 experiment with constant (240 ppm) CO 2 concentration.

Figure 11 .
Figure 11.δ 13 C and radiocarbon ventilation age distribution in the Atlantic Ocean in the ONE_1.0_130Ksimulation of Termination I. (a, b) δ 13 C (‰), (c, d) radiocarbon ventilation age in yr 14 C. (a, c) Differences between LGM (21 ka) and pre-industrial in the Atlantic Ocean.(b, d) Temporal evolution of anomalies during the past 20 kyr at 20 • N in the Atlantic.

Figure 13 .
Figure 13.Vertical profile of ventilation age in 14 C years for Atlantic (a), Pacific (b), and Southern Ocean (c).Red line represents modern conditions; solid blue -LGM in ONE_1.0_130Kexperiment using the standard version of the model; dashed blue -LGM in ONE_BRINE_130K experiment with the model version, which includes brine parameterisation and stratification-dependent vertical mixing.Red (blue) squares represent basin-averaged radiocarbon age for the modern (LGM) state based on the data of Skinner et al. (2017).

Figure 14 .
Figure 14.Dynamics of terrestrial carbon pools (Gt C) in the one-way coupled ONE_1.0 simulation.Left panels: the whole 400 kyr period; right panels: the Termination I period.(a) Black line -total carbon storage; magenta line -conventional carbon pools (biomass and mineral soils).(b)-(d) Peat, permafrost, and buried carbon storages, respectively.

Table 1 .
Model experiments performed in this study.P denotes prescribed characteristic; I -interactive; STD -standard model configuration; RD -variable remineralisation depth; VO -variable volcanic outgassing; IF -iron fertilisation in the Southern Ocean; TC -terrestrial carbon cycle; BR -brine rejection mechanism.A minus sign means that the process is excluded and a plus sign means that process is included.Ice sheets are interactive in all simulations.

Table 2
. "The carbon stew" at the LGM and the entire 400 kyr period in the ONE_1.1 experiment.Positive (negative) sign indicates that the process contributes to glacial CO 2 reduction (increase) as compared to the pre-industrial concentration simulated at the end of the experiment.* Deep ocean and shallow water carbonate sediments, carbonate and silicate weathering.
shows that in the Atlantic and Pacific oceans, even the standard model version overestimates the radiocarbon ventilation age of glacial water masses.In turn, the model version with the brine rejection parameterisation simulates water masses which are 500 to 1000 years older than in the standard version.Only in the Southern Ocean is the reconstructed ventilation age consistent with both model versions.As a result, the standard model version simulates ca.800 14 C yr increase in glacial ocean ventilation age at the