Glacial CO 2 cycle as a succession of key physical and biogeochemical processes

Introduction Conclusions References


Introduction
Ice core records of atmospheric CO 2 concentration through the last 800,000 years show the carbon cycle amplifying the climate forcing from variations in Earth's orbit (Augustin et al., 2004).This positive climate-carbon cycle feedback could weaken or even possibly reverse present-day fossil fuel CO 2 uptake by the natural carbon cycle.Despite much effort over the last decades, a process-based explanation of the carbon cycle feedbacks responsible for the glacial-interglacial CO 2 cycles remains elusive.Mechanisms proposed to explain the glacial cycles of atmospheric CO 2 include changes in sea surface temperatures (SSTs), nutrients, circulation, and carbonate chemistry (Archer et al., 2000;Sigman and Boyle, 2000;Toggweiler et al., 2006;Kohfeld and Ridgwell, 2009).The main focus has been to explain the low atmospheric CO 2 concentration of about 190 ppmv at the last glacial maximum (LGM) about 21 000 years ago.Most of these studies were built on the assumption that the carbon cycle was in equilibrium at that time (e.g.Hain et al., 2010;Kurahashi-Nakamura et al., 2010).However, the carbon cycle is most likely never in steady state during a glacial cycle because of the long time scales of many of the biogeochemical processes.Our goal is to explain the CO 2 concentration not only for the LGM snapshot but through the whole glacial cycle, focusing on a possible sequence of mechanisms that drive CO 2 changes (Köhler and Fischer, 2006;Peacock et al., 2006;Watson et al., 2000).

V. Brovkin et al.: Glacial CO 2 cycle as a succession of key physical and biogeochemical processes
The challenge is to explore these mechanisms using the Earth system models that include 2-or 3-dimensional models of oceanic circulation, a key driver of oceanic biogeochemistry.

The model
The CLIMBER-2 model includes six components of the Earth system: atmosphere, ocean, sea ice, land surface, terrestrial vegetation and ice sheets.The first five components are represented by coarse-resolution modules of intermediate complexity (Brovkin et al., 2002;Petoukhov et al., 2000).The ice-sheet component is represented by the three-dimensional polythermal ice-sheet model SICOPOLIS (Greve, 1997).The ice-sheet model is applied only to the Northern Hemisphere (NH).The carbon cycle model includes land carbon, oceanic biogeochemistry, a model for marine biota, and a sediment model (Archer, 1996;Brovkin et al., 2007).Weathering rates scale to runoff from the land surface (Munhoven, 2002), the coral reef growth depends on the sea level rise (Kleypas, 1997), the scale of nutrients utilization in sub-Antarctic ocean is proportional to the dust deposition derived from Antarctic ice core (Augustin et al., 2004), and volcanic outgassing has been assumed to have remained constant through the whole glacial cycle.The carbon cycle components in past studies were run off-line, i.e. simulated changes in atmospheric CO 2 were not fed back into the radiative calculation, which evolved as a function of the history of equivalent CO 2 derived from ice core data.

Experimental setup of the physical model
The experimental set-up followed a recent simulation of the last glacial cycle without carbon fluxes (Ganopolski et al., 2010).We prescribed temporal variations in the orbital parameters (Berger, 1978) and the equivalent CO 2 concentration, which accounts for changes in three major greenhouse gases (GHGs -CO 2 , N 2 O and CH 4 ) derived from the Antarctic ice cores (Barnola et al., 1987;Augustin et al., 2004;Monnin et al., 2004), and aeolian dust deposition (Mahowald et al., 1999) (Fig. 1).As initial conditions for the physical climate system, we used equilibrium climate and ice sheets simulated for the present-day orbital configuration and GHG concentrations.The model was run from 126 kyr BP for 136 000 years, i.e. for 10 000 years into the future.This was done to estimate the equilibration time of the global carbon cycle after the last deglaciation.For the future, GHG concentrations were kept constant at the preindustrial level, i.e. the human influence on climate was not considered.Under the influence of the orbital forcing alone, CLIMBER-2 simulates a slow growth of the ice sheets in the Northern Hemisphere which is explained by favorable orbital configuration but small eccentricity and, therefore, weak orbital forcing.Simulated changes in climate and ice volume over the past 126 000 years agree favourably with paleoclimate reconstructions (Ganopolski et al., 2010, see Fig. 1d).

Ocean carbon cycle
The physical ocean component of the model is essentially the same as the one used in a previous glacial equilibrium simulation (Brovkin et al., 2007).The only important difference is that the background vertical diffusivity of the biogeochemical tracers in the ice-free Southern Ocean south of 50 • S was enhanced by an order of magnitude, i.e. to ca. 10 −3 m 2 s −1 (under the sea ice, the standard values of 10 −4 m 2 s −1 was retained).This was done to obtain CO 2 outgassing in the Southern Ocean (SO) consistent with recent empirical estimates of about 0.4 GtC yr −1 (Gruber et al., 2009).The standard version of CLIMBER-2, similar to many other even much more sophisticated ocean models (Fletcher et al., 2007), considerably underestimates the magnitude of CO 2 outgassing in the SO.We found that enhanced outgassing in the SO considerably increased the sensitivity of the atmospheric pCO 2 to changes in sea ice cover in the SO.For example, in the standard model version (without enhanced vertical diffusivity in the SO) under present-day climate conditions, the complete insulation of the SO (south of 50 • S) from the atmosphere leads to changes in atmospheric pCO 2 by only 4 ppm (Archer et al., 2003), while in the version with enhanced outgassing, the SO insulation leads to a 20 ppm decrease in atmospheric CO 2 .The model version with the enhanced southern outgassing simulates glacial/interglacial CO 2 variations of ca. 10 ppm higher than the standard model version.Enhanced diffusivity was only applied to the passive (biogeochemical) tracers.Although the effect of enhanced vertical diffusivity on active tracers (temperature and salinity) is small, we prefer to keep the physical component of the model unchanged to avoid inconsistencies with our previous results (Brovkin et al., 2007).
To evaluate the effect of modified parameter values on the biogeochemistry, we performed an equilibrium simulation of the model in the pre-industrial setup.In the Annex, the model performance is compared with observations in terms of distribution of essential biogeochemical tracers (PO 4 , dissolved organic carbon, alkalinity, and 14 C).Transient simulations of the oceanic carbon and CFC-11 uptakes throughout the last century are provided for comparison as well.In general, the performance of the oceanic biogeochemistry model for preindustrial state is very similar to the results from the previously published CLIMBER-2 versions (Brovkin et al., 2002(Brovkin et al., , 2007)).

Inventories of 12 C and 13 C
Global inventories of 12 C and 13 C in the model experiments are conservative.The atmospheric CO 2 concentration was calculated from the total carbon inventory equation: where the index i denotes the carbon isotope ( 12 C or 13 C), i C stands for carbon isotope inventories (GtC) and i F is for carbon isotope fluxes (GtC yr −1 ).Indices "atm", "oc", and "land" refer to atmospheric, ocean, and land compartments, respectively, t is the model year starting from 0 at 126 000 yr BP, and the index 0 for carbon compartments is for the initial storages of carbon at the beginning of the simulations.i F volc is for volcanic outgassing of CO 2 fixed to a value 0.066 GtC yr −1 with 13 C content of 6 ‰ (relative to the Pee Dee Belemnite -PDB -standard).i F atm cons is the atmospheric consumption of CO 2 due to terrestrial weathering of both carbonate and silicate rocks, i F weath is the bicarbonate flux from the land to the ocean arising from continental weathering.The consumption of atmospheric CO 2 and the production of riverine HCO − 3 by continental weathering processes are calculated as a function of geographically distributed runoff (interactively simulated by the land surface module) and lithology (Amiotte-Suchet et al., 2003).i F burial is an interactively simulated CaCO 3 burial on shelves (coral reefs) (Kleypas, 1997) and in the deep sea (Archer, 1991).Ocean carbon storage i C oc is for dissolved inorganic and organic carbon components in the ocean, including marine biota biomass.Land carbon compartment i C land is for vegetation biomass and soil carbon storages.In the land carbon cycle, we assumed that the decomposition of soil organic matter is slowed if annual mean temperatures fall below −10 • C.
In Eq. ( 1), the land surface area accounted in the land carbon storage i C land , weathering i F atm cons , and atmospheric CO 2 consumption i F atm cons was linearly scaled between pre-industrial and LGM land surface areas proportionally to the sea level.The ocean volume was also adjusted in accordance with the sea level.The total ocean alkalinity is calculated as a balance of alkalinity fluxes due to terrestrial weathering and burial of CaCO 3 on shelves and in the deep ocean.
Atmospheric CO 2 calculated from Eq. (1) (using a conversion factor of 0.47 ppm GtC −1 ) was not used as a radiative forcing for the atmospheric model component, which was driven by greenhouse gas concentrations reconstructed from Antarctic ice cores (see Sect. 2.2).This methodology assures that the carbon cycle is influenced by changes in the physical climate system.It should, however, be emphasized that there is no feedback of the simulated atmospheric CO 2 concentration to climate.Although this approach is obviously limited, it is a useful first step in the analysis of glacial CO 2 cycles within the model framework avoiding the effect of carbon cycle biases on climate.

Initial conditions of the carbon cycle
Initial conditions of the land and oceanic biogeochemistry were set up following a procedure used recently for initializing interglacial carbon cycle (Kleinen et al., 2010).Firstly, the model was run, with equilibrium conditions of 126 000 BP with CO 2 kept at the level of 280 ppmv.Ocean alkalinity was increased to get a carbonate sedimentation flux of 16 T mol yr −1 in the deep ocean and 2 T mol yr −1 on the shelves in order to simulate the maximum in CaCO 3 preservation in the deep sea before the onset of the interglacial.Atmospheric δ 13 CO 2 was initialized to −6.7 ‰.The model was run with prescribed CO 2 for 5000 years.In the transient simulations through the last glacial cycle, the carbonate accumulation on the shelves and in the deep sea was interactive with a shift towards carbonate accumulation in shallow  (Ahn and Brook, 2008;Barnola et al., 1987;Monnin et al., 2004) (thin red lines) in ppmv; (c) atmospheric CO 2 concentration for simulations P, PC, and PCB (brown, blue, and green lines, respectively) in ppmv; (d) deep-sea δ 13 C for simulation PCBL (thick black line) and geological stack record for the Atlantic Ocean (Lisiecki et al., 2008) (thin pink line) in ‰; (e) the same as (d) but for Pacific Ocean.waters during interglacials.This setup of initial conditions ensures that the ocean biogeochemistry is in equilibrium with the climate at the onset of interglacial while it is in transition from the glacial to interglacial state hereafter.

Transient experiments
The physical part of the CLIMBER-2 model has been used to simulate the last glacial cycle, forced by changes in orbital forcing and greenhouse gases concentrations (Ganopolski et al., 2010).The model reproduces the temporal and spatial dynamics of the major Northern Hemisphere (NH) ice sheets, including the rapid glacial inception started at 120 to 118 kyr BP (thousand years before present) (Fig. 1d).The global mean annual temperature decreases by about 6 K at the LGM relative to pre-industrial (Fig. 1e) while the annual mean sea ice area in the Southern Hemisphere increases three-fold during the same period (Fig. 1f).
Here, we focus on the responses of various parts of the carbon cycle to the changes in the physical system.Simulation P includes the effects of physical mechanisms only (changes in SST, circulation, and sea level), PC adds changes in the carbonate chemistry (carbonate compensation, shift of carbonate sedimentation from shallow waters to the deep sea and vice versa, and weathering), PCB complements with changes in nutrient utilization by marine ecosystems (iron fertilization due to aeolian dust deposition), and PCBL finally adds land carbon dynamics.All simulations start from the same non-equilibrium conditions in oceanic carbonate chemistry at 126 kyr BP and continue to 10 000 years after the present day without considering anthropogenic forcings.The atmospheric CO 2 response of the full PCBL simulation is shown in black in Fig. 2b, and the responses are broken down into component mechanisms in Fig. 2c.
During the last interglacial, CO 2 dynamics in the full experiment, PCBL, increases from an initial value of 280 ppmv at 125 kyr BP to 285 ppmv at 118 kyr BP (Fig. 2b).This is driven by the response of the terrestrial biosphere to climate trends during the late interglacial, as the drying of the northern subtropical regions and the southward retreat of boreal forests decrease land carbon storage by about 300 GtC during this period (Fig. 3).In the P simulation, the increased CO 2 solubility in seawater due to lowered SSTs leads to an atmospheric CO 2 decrease during the late interglacial, while carbonate chemistry processes, including the accumulation of shallow-water carbonates, force CO 2 in the opposite direction leading to CO 2 stagnation till 117 kyr BP in the PC and PCB simulations (Fig. 2c).
The glacial inception results in an abrupt CO 2 drop by ca.50 ppmv to the level of 230-235 ppmv (Fig. 2b) by 110 kyr BP, which is in a good agreement with the Vostok ice core data (Barnola et al., 1987).The main mechanism here is physical (simulation P).It is the reduction in the SSTs and the substitution of North Atlantic Deep Water by colder and saltier waters of Antarctic origin that lead to a drop of ca. 30 ppmv (Brovkin et al., 2007).The essence of the later mechanism, called the "standing volume effect" (Skinner, 2009), is that the dissolved inorganic carbon (DIC) content of southern-source water is higher than that of North Atlantic deep water, leading to increased DIC storage in the deep ocean.The model during glacial inception simulates a decrease in δ 13 C in deep tropical Atlantic by 0.2 ‰ qualitatively in agreement with the benthic δ 13 C data (Lisiecki et al., 2008;Oliver et al., 2010) (Fig. 2d).
During the rest of the glacial period, physical mechanisms only have a substantial impact on atmospheric CO 2 during Heinrich and Dansgaard-Oeschger type events, when the Atlantic meridional overturning circulation (AMOC) is greatly perturbed (Fig. 2a).During these events, the model produces spikes of CO 2 of about 10 ppmv (Fig. 2c) due to outgassing of oceanic CO 2 , mostly via the Southern Ocean (Anderson et al., 2009).At the same time, the deep ocean is enriched in negative δ 13 C due to increased water mass ages (Fig. 2d  and e).Millennial-scale variability in atmospheric CO 2 during abrupt changes in AMOC is visible in the high-resolution Byrd ice core data (Ahn and Brook, 2008).The specific timings of the abrupt climate events are different from the observed ones as they depend on freshwater discharge into the North Atlantic, which are random output resulting from the instability of the Laurentide ice sheet within the model.
By 110 kyr BP, carbonate compensation amplifies the physically driven changes in atmospheric CO 2 by an additional 10 ppmv (Fig. 2c, the difference between P and PC).Between 110 and 100 kyr BP, atmospheric CO 2 rises by 20 ppmv as a consequence of partial deglaciation due to increasing summer insolation in boreal region (Fig. 1d).This CO 2 response to the precessional cycle is not evident in the ice core data (Barnola et al., 1987), although the coarse temporal resolution of the only available Vostok record that might be insufficient for resolving the cycle.Between 90 and 70 kyr BP, the model yields atmospheric CO 2 concentrations at the level of 230-240 ppmv, in line with the Byrd ice-core record (Ahn and Brook, 2008).After 70 kyr BP, the effect of fertilization of marine biological production in sub-Antarctic Atlantic ocean by aeolian dust leads to a CO 2 drop to a level of 200-210 ppmv (Fig. 2c).The timing of the dust effect on the carbon cycle is in line with box model simulations (Hain et al., 2010;Martinez-Garcia et al., 2009;Ridgwell and Watson, 2002;Watson et al., 2000), although the magnitude of the impact of iron fertilization differs among models (e.g, Parekh et al., 2008).Enhanced marine productivity leads to an additional CO 2 drop with minima of pCO 2 at 63 and 20 kyr BP in simulations PCB and PCBL (Fig. 2c).These minima are also pronounced in the deep ocean δ 13 C, both in the model and in the data (Fig. 2d and e).During the remaining period of 60-20 kyr BP, modelled atmospheric CO 2 stays in the range of 200-220 ppmv, with one exceptional spike at 54 kyr BP when CO 2 quickly rises to 230 ppmv during a Heinrich-type event and relaxes slowly to the level of 220 ppmv thereafter.The model atmospheric CO 2 concentration during the Last Glacial Maximum was about 200-210 ppmv, ca. 10 ppmv higher than in the EPICA Dome C (EDC) record (Monnin et al., 2004).In the equilibrium model simulation (Brovkin et al., 2007), the CO 2 level was higher (about 220 ppmv).In equilibrium experiments, the total ocean alkalinity is implicitly prescribed by initial conditions; whereas in transient simulation, it is calculated interactively as a cumulative balance between weathering and sedimentation fluxes, which vary in time.
Deglaciation in the model starts by ca.20 kyr BP in line with the sea level data (Fig. 1d).The increased freshwater flux to the North Atlantic shuts down the AMOC (Fig. 2a).In response, atmospheric CO 2 rises by about 10 ppmv (Fig. 2b).Between 18 kyr BP and 11 kyr BP, CO 2 rises from about 200 to 235 ppmv, which is less than the steep increase in the EDC data by about 50 ppmv.A substantial part of the CO 2 increase (ca.20 ppmv) is due to reduced marine productivity in the sub-Antarctic region (Fig. 2c, PCB).After the AMOC recovers, the ocean circulation returns to an interglacial condition with enhanced North Atlantic Deep Water formation, interrupted by a short shutdown at 11 kyr BP.
In the PCBL simulation, accumulation of CaCO 3 in the deep tropical Pacific decreases during both the last interglacial (126-115 kyr BP) and the Holocene in response to enhanced shallow-water CaCO 3 burial on shelves (Fig. 4).During the glacial period, CaCO 3 burial increased, with a peak of preservation during deglaciation at about 15 kyr BP, in line with data from the deep South Atlantic (Hodell et al., 2001;Anderson et al., 2008).A key diagnostic of ocean carbonate chemistry, the concentration of carbonate ion, [CO in the data for depths below 3 km (Fig. 6, left panel).In the tropical Atlantic, data for the LGM reveals increased dissolution at almost all depths.The model simulates a small decrease in % CaCO 3 below 4 km depth and a small increase above this depth (Fig. 6, right panel).Although the model does not capture well the relative changes in CaCO 3 dissolution between the LGM and the pre-industrial in tropical Atlantic, it does reproduce a different response of carbonate preservation in the deep Atlantic and Pacific to glacial conditions.
The main forcing of atmospheric CO 2 after 10 kyr BP is a reduction of ocean alkalinity resulting from coral reef regrowth (Kleypas, 1997) and its impact on the CaCO 3 cycle (Ridgwell et al., 2003).The terrestrial biosphere ameliorates the CO 2 rise by taking about 500 GtC in between 18 and 5 kyr BP (Fig. 3) in response to changes in climate, atmospheric CO 2 concentration, and sea level rise.By 0 kyr BP, the carbonate chemistry is still not fully equilibrated to interglacial conditions, as atmospheric CO 2 , now at 274 ppmv, continues to rise.Extending the model runs by 10 kyr into the future (without anthropogenic perturbations) allows the model to find an interglacial equilibrium of 279 ppm at 4 kyr in the future.The main difference between this state and the previous interglacial state at 126 kyr is that the terrestrial biosphere carbon storage is about 250 GtC smaller, due to more favourable conditions of the past interglacial relative to  the Holocene (Fig. 3).Recent evidence (Zimov et al., 2009) suggests that several hundred GtC might have been stored in the high latitude regions in permafrost soils during the glacial period.Therefore, although model experiments suggest that the land biosphere had 500-800 Gt less carbon at the LGM than during the Holocene (Kaplan et al., 2002), the simulated reduction in glacial land carbon storage could possibly be overestimated since these models neglect changes in permafrost carbon storage.
A substantial part of the glacial-interglacial changes in the carbonate chemistry in the PC, PCB, and PCBL simulations is due to the changes in the continental weathering (Fig. 7).The global riverine bicarbonate (HCO − 3 ) flux increases from 36-38 T mol yr −1 in interglacial states to almost 50 T mo yr −1 in the LGM.Changes in the HCO − 3 flux result from changes in the runoff distribution (globally reduced during glacial periods) and in the area of exposed continental shelves (increased during glacials).The relevant parts of the latter are mainly located in the tropical regions where shelves are mainly covered with carbonates.The main increase in the global weathering during glacials comes from carbonate weathering, as documented by the atmospheric CO 2 consumption rates depicted in Fig. 7.The increase in carbonate weathering is reflected by the difference between the global CO 2 consumption and silicate weathering, or, equivalently, between the global river HCO − 3 flux and the global CO 2 consumption.The change in the CO 2 consumption rate due to silicate weathering is less pronounced.It varies between 12-13 T mol yr −1 for interglacial periods and 10 T mol yr −1 for the last glacial maximum (Fig. 7).Silicate weathering varies out of phase with total CO 2 consumption and HCO − 3 production rates (both dominated by carbonate weathering variations).This is because carbonate weathering variations are mainly driven by the shelf exposure mainly located in the tropics, while silicate weathering fluxes are related to the expansion and retreat of the ice-sheets as well as the runoff change over the continental areas that remain exposed throughout the entire simulation.

"Standing volume" effect on the atmospheric CO 2 concentration
For a better understanding of the physical mechanisms responsible for the glacial CO 2 drawdown, it is useful to compare the standard experiment P, which only includes physical processes, with an even more schematic experiment P NOI.
In that experiment the ice sheet component was switched off and therefore neither climate changes due to growth and decay of the ice sheets nor variations in the sea level and ocean volume affect the carbon cycle.Therefore, all climate changes that affect the carbon cycle in this experiment are only driven by orbital forcing and radiative forcing of the changing CO 2 concentration.Note that the direct effect of orbital forcing on climate and CO 2 is rather small.In all other respects, the experiment P NOI is equivalent to the experiment P. Simulated CO 2 variations for the P and P NOI simulations are shown in Fig. 8 by black and pink lines, respectively.The simulated atmospheric pCO 2 in the P NOI experiment closely resembles the prescribed equivalent CO 2 concentration but with about a three times smaller amplitude of variation (Fig. 8).The maximum glacial/interglacial CO 2 difference in this experiment is about 40 ppm.Since variations in the global mean SST are about 2 • C (similar changes are simulated in the global surface air temperature and the deep ocean temperatures), this implies a CO 2 sensitivity to temperature changes of about 20 ppm • C −1 , which is too high for the solubility effect alone (Archer et al., 2004).In fact, the solubility effect explains only ca. 10 ppm, while the rest results from changes in the redistribution of major water masses and their ventilation.The first effect is now commonly known as the "standing volume effect" (Skinner, 2006(Skinner, , 2009)).This effect was already discerned in previous equilibrium glacial CLIMBER-2 simulations (Brovkin et al., 2007).The essence of this effect is that under glacial climate conditions, the volume of the water masses in the World Ocean originating from the North Atlantic decreases and the volume of the water masses forming in the SO increases.Since the DIC concentration is higher in the water masses originating in the SO than in NADW, this leads to an increase of carbon storage in the deep ocean.To illustrate this effect we have added a dye tracer to the ocean which allows us to compute the relative contributions of the water masses of each origin at any location.Figure 9 shows that in the interglacial climate state, most of the deep Atlantic is filled with water masses originating in the North Atlantic.Under glacial climate conditions, the upper cell of the AMOC becomes shallower and more water masses of southern origin penetrate into the deep Atlantic.As a result, the fraction of the northern water masses in the deep Atlantic decreases from almost 100 % to 50 %.Smaller but non negligible changes were also found in the Pacific and Indian oceans (not shown).Note that the increase of the fraction of the water masses of southern origin is not related to the changes in the rate of Antarctic bottom water (AABW) formation but to a relative increase in the density of AABW compared to NADW.This, in turn, is primarily related to the enhanced sea ice formation in the Southern Ocean, stronger brine rejection and an increase of the AABW salinity (Brovkin et al., 2007).Increased salinity of AABW qualitatively agrees with empirical reconstructions (Adkins et al., 2002).Since glacial/interglacial climate change in the Southern Ocean is primarily explained by changes of GHG concentrations, the increase of the sea ice area in the P NOI experiment is similar to that in the P experiment with the interactive ice sheets.As a result, the fraction of northern water masses in the World Ocean (below 1.5 km) in the P NOI experiment decreases at the LGM by ca. 10 % compared to the interglacial state.During the first 100 ka, the changes in the fraction of the North Atlantic water masses are similar in the P and in the P NOI experiments, but around the LGM, the fraction of water masses originated in the North Atlantic falls drastically in the P experiment, which is related to further weakening and frequent shutdowns of the AMOC at this period.This further reduction in the fraction of North Atlantic water masses, however, does not result in any further CO 2 drop.In fact, the CO 2 concentration at this period of time in the P NOI experiment is even lower than in the P simulation.This can be explained by the fact that the weakening of the AMOC in the P experiment results in a reduction of the CO 2 uptake in the North Atlantic, which in turn leads to a decrease of DIC in the intermediate water masses and compensates for the increase in the fraction of the water masses of southern origin.Recently, Bouttes et al. (2011a) used another version of CLIMBER-2 model to analyse equilibrium sensitivity of the glacial carbon cycle to brine-induced stratification, stratification-dependant diffusion, and iron fertilization.They demonstrated that, with certain parameterizations of these processes, the model is able to reproduce the atmospheric CO 2 level and the ocean carbon isotopic composition observed for the LGM.The results of Bouttes et al. (2011a) can possibly be interpreted in terms of the "standing volume effect".However, such an analysis would require additional experiments which go beyond the purpose of this study.
Changes in the "standing volume" also have an effect on the oxygen content of the intermediate and deep waters in the tropical Pacific (Fig. 10).In the interglacial states, these waters are well ventilated (Fig. 10, top panel).During the glacial state, the oxygen concentration declines as a result of the increased fraction of waters originating in the SO, the decrease of Southern water ventilation and the slowdown of the circulation in the Pacific basin.This effect is slightly more pronounced in the case where we take the enhancement by the marine biology into account (PCBL simulation, see Fig. 10, bottom panel).The large-scale anoxia in the tropical Pacific during glacial states is not supported by paleoproxies, although there are data that suggest that the bulk of the deep ocean was less oxygenated during the LGM (Jaccard and Galbraith, 2012).It is important to note that oxygen concentration simulated in the model deep Pacific for the preindustrial is about two times less than observed.Therefore, we think that simulated trends are more reliable than the absolute values for the oxygen content and that the simulated lowering of oxygen content during glacial times is not necessarily in contradiction with the available data.Relative to Pacific, the oxygen content in the Atlantic basin is much higher throughout the simulations.A considerable reduction in the oxygen content (but still well above anoxia level) in the tropical Atlantic is simulated only during periods of a shutdown of the AMOC.The oxygen content restores rapidly after the circulation recovers (not shown).

Millennial-scale variability in atmospheric CO 2
Apart from climate changes on the orbital time scales, the CLIMBER-2 model simulates pronounced millennial scale variability associated with changes in the AMOC strength caused by the changing freshwater flux into the northern North Atlantic (Fig. 2a).The largest fluctuations of the freshwater flux are associated with the internal instability of the Laurentide ice sheet and the development of fast ice streams over the Hudson Strait (Calov et al., 2002).These events of enhanced ice flux resemble real Heinrich events and result in additional freshwater fluxes of about 0.1 Sv during several hundreds to one thousand years.Since Heinrich-like events simulated in the model result from internal instability of the ice sheet component, their timing is random and does not coincide with the timing of real Heinrich events.During these intervals of enhanced freshwater flux into the North Atlantic, the AMOC weakens considerably.The NADW formation area retreats southward and the AMOC cell shoals (Ganopolski and Rahmstorf, 2001).These changes of the AMOC are closely associated with millennial scale variations in the atmospheric CO 2 concentration (Fig. 2b and c).During the phases of weak AMOC, the CO 2 concentration rises rapidly by 10-20 ppm and equally rapidly drops after the resumption of strong AMOC.The temporal dynamics and the magnitude of the simulated millennial CO 2 variations resemble CO 2 variations in the Byrd ice core record during Heinrich stadials (Fig. 2b).
The strong positive excursions in the CO 2 concentration during Heinrich-like events can already be seen in the P experiment.The magnitude of these variations is essentially the same in the experiment with and without terrestrial biosphere.Therefore, the rise of CO 2 concentration during the weak AMOC in our model can be primarily attributed to the ocean component of the carbon cycle.The small contribution of the terrestrial biosphere to the CO 2 change can be explained by the absence of terrestrial vegetation in the areas primarily affected by the reorganization of the AMOC under the glacial climate conditions.The detailed analysis of the mechanisms behind the millennial scale variability of the atmospheric CO 2 concentration reveals considerable complexity and a number of counteracting mechanisms involved.

V. Brovkin et al.: Glacial CO 2 cycle as a succession of key physical and biogeochemical processes
During periods of weakened AMOC, the carbon uptake in the northern North Atlantic is significantly reduced and DIC in the core of NADW decreases.However, a reduction of the carbon content in the upper part of the Atlantic is compensated by an increase of DIC in the deep Atlantic due to the substitution of NADW by water masses of southern origin which have relatively high DIC concentration.As a result, the Atlantic Ocean, although characterized by the largest changes in the atmosphere-ocean carbon fluxes, does not contribute much to the CO 2 rise during stadials.In the Southern Ocean, neither the surface fluxes nor the total carbon inventory are appreciably affected by changes in the AMOC in spite of significant variations in the surface temperature and sea ice cover caused by the reduction of the northward ocean heat flux (i.e. the "bipolar seesaw" mechanism).This is because the increase of the oceanic outgassing in the area where the southern sea ice cover is retreated polarwards is counteracted by an increase in the strength of the biological pump.
There are little changes in the surface ocean-atmosphere carbon fluxes over the Indo-Pacific but the carbon inventory indicates that the rise of the atmospheric CO 2 during periods of weakened AMOC is explained by the loss of carbon from the Indo-Pacific.The decrease in DIC occurs here primarily between 1 and 3 km depth and can be explained by a weakening of the reverse cell of the Indo-Pacific overturning circulation during periods of reduced AMOC.This leads to a reduction of the vertical nutrient transport and, as a result, to a weakening of the biological carbon pump.This mechanism is consistent with the results of Schmittner (2005) who found similar changes in the overturning circulation, nutrient distribution and primary production in the Indo-Pacific.Therefore, the mechanisms of the millennial and of the orbital scale variability of the atmospheric CO 2 in our model are fundamentally different.
Previous studies of the carbon cycle response to a weakening or a complete shutdown of the AMOC produced controversial results.While Marchal et al. (1999), Schmittner andGalbraith (2008), Chikamoto et al. (2008) and Bouttes et al. (2011b) found that the ocean was releasing CO 2 during weakenings of the AMOC, Obata (2007), Menviel et al. (2008) and Bozbiyik et al. (2011) reported opposite results.We conclude that the complexity of the response of atmospheric pCO 2 to such changes in the AMOC and a number of counteracting processes taken into account in our model may explain why this response is so strongly model-dependent.

Conclusions
According to our results, physical mechanisms -reduction in SSTs and changes in standing volume in response to the expansion of ice-sheets, especially in the Northern Hemisphere, and the resulting lower sea level -lead to the initial drop of CO 2 during a glacial inception.This drop, amplified by the CaCO 3 cycle, is followed by an increased buildup of ocean alkalinity during the rest of the glacial cycle, as exposed tropical shelves serve as a source of CaCO 3 to the ocean, and CaCO 3 burial is shifted from shallow waters to the deep sea.Increased nutrient utilization in the sub-Antarctic plays the dominant role in further CO 2 decline.The role of the land carbon is less certain, as recently published data suggest globally smaller changes in the land carbon storage between glacial and interglacial times than previously thought.
In the second part of the glacial cycle, our model simulates millennial scale variations in the atmospheric CO 2 concentration that appear to be closely associated with AMOC changes.In response to AMOC weakenings, the CO 2 concentration rapidly rises, similarly to CO 2 variations observed during Heinrich stadials (Ahn and Brook, 2008).During these variations, the Indo-Pacific Ocean serves as a source of CO 2 , while the land carbon storage increases in response to the elevated CO 2 concentration.Distributions of carbon sinks and sources in response to glacial AMOC changes differ among models due to counteracting effects of physical and biogeochemical mechanisms.A detailed analysis of the carbon cycle response to AMOC changes is certainly an interesting topic for future model intercomparison exercises.
During the deglaciation and the Holocene, the model simulates an increase in the atmospheric CO 2 concentration that is of the right amplitude but less abrupt than the rapid CO 2 rise recorded in the ice cores (Monnin et al., 2004).This mismatch could be due to limitations in the representation of ocean biogeochemistry, especially carbonate chemistry, in our zonally averaged model or to abrupt processes not considered in our land carbon model, such as rapid thawing of frozen organic matter stored in permafrost soils.Implementation of permafrost carbon could possibly also improve the CO 2 dynamics during glacial inception.Constraining the amplitude of such a permafrost carbon feedback from paleo simulations is essential for projecting the response of the carbon cycle at high latitudes to anthropogenic climate change over the next centuries (Schaefer et al., 2011).
In this study, we did not consider the feedback from simulated atmospheric CO 2 to the physical climate model.This partial decoupling between climate and carbon cycle is necessary to evaluate a role of different mechanisms in the glacial CO 2 changes within the model framework.These numerical experiments serve as a basis for further simulations of glacial cycles driven solely by changes in orbital forcing.presents present-day observations (GEOSECS) zonally averaged and interpolated onto the model grid.In general, the performance of the oceanic biogeochemistry model is very similar to the previously published versions (Brovkin et al., 2002(Brovkin et al., , 2007)).The current version tends to simulate older waters in the deep Northern Pacific and to accumulate more PO 4 in the North Pacific.
V. Brovkin et al.: Glacial CO 2 cycle as a succession of key physical and biogeochemical processes
The ocean carbon uptake from the year 1800 to 1994 with a prescribed history for the concentration of CO 2 in the atmosphere is 106 GtC.This is lower than the estimate of 118 GtC from Sabine et al. (2004), but within the range of uncertainty of ±19 GtC.The average carbon uptake in the 1990s is 2.2 GtC yr −1 , which is fully in agreement with the IPCC AR4 estimate of 2.2 ± 0.4 GtC yr −1 (Denman et al., 2007).In the decade from 2000 to 2009, the average uptake rises to 2.5 GtC yr −1 (Fig. A3).

Fig. 1 .
Fig. 1.Prescribed forcings and results of physical model simulation of the last glacial cycle.(a) Orbital forcing (illustrated by the June insolation at 65 • N); (b) radiative forcing of GHGs, (c) concentration of dust in the EPICA ice core; (d) changes in the simulated ice sheet volume (thick black line) and reconstructed sea level change (Waelbroeck et al., 2002) (thin red line), in meters of sea level change relative to 126 kyr BP; (e) changes in the simulated global surface air temperature; (f) simulated annual mean sea ice area in the Southern Hemisphere.

Fig. 2 .
Fig. 2. Comparison of model results with geological records: (a) simulated Atlantic meridional overturning circulation in Sv; (b) atmospheric CO 2 concentration, simulation PCBL (thin black line), and Antarctic ice core reconstructions(Ahn and Brook, 2008;Barnola et al., 1987;Monnin et al., 2004) (thin red lines) in ppmv; (c) atmospheric CO 2 concentration for simulations P, PC, and PCB (brown, blue, and green lines, respectively) in ppmv; (d) deep-sea δ 13 C for simulation PCBL (thick black line) and geological stack record for the Atlantic Ocean(Lisiecki et al., 2008) (thin pink line) in ‰; (e) the same as (d) but for Pacific Ocean.

Fig. 3 .
Fig. 3. Dynamics of the total terrestrial carbon storage (soil plus biomass in GtC) in the PCBL simulation.

Fig. 4 .
Fig. 4. Dynamics of the vertical distribution of the simulated CaCO 3 fraction (%) in top layer of seafloor sediments averaged over 30 • S-30 • N in the Pacific Ocean.

Fig. 7 .
Fig. 7. Dynamics of the riverine bicarbonate flux (blue), the global atmospheric CO 2 consumption rate (brown) and silicate weathering CO 2 consumption rate (green) in the PCBL simulation (in T mol yr −1 ).

Fig. 8 .
Fig. 8. Relative volume of the water masses of the northern origin in the Atlantic ocean (a) and simulated atmospheric CO 2 concentration (b).Black line -experiment P , pink line -experiment P NOI.

Fig. 9 .
Fig. 9. "Dye" tracer concentration (shading) and meridional overturning circulation (in Sv, contours) in the Atlantic for interglacial (a) and LGM (b) states.Northern water masses have a "dye" concentration of 1 and the southern water masses of 0.

Fig. 10 .
Fig. 10.Dynamics of the vertical distribution of the simulated O 2 concentration (µmol kg −1 ) averaged over 30 • S-30 • N in the Pacific Ocean.