Response of the carbon cycle in an intermediate complexity model to the different climate configurations of the last nine interglacials

Atmospheric CO₂ levels during interglacials prior to the Mid-Brunhes Event (MBE, ∼ 430 ka BP) were around 40 ppm lower than after the MBE. The reasons for this difference remain unclear. A recent hypothesis proposed that changes in oceanic circulation, in response to different external forcings before and after the MBE, might have increased the ocean carbon storage in pre-MBE interglacials, thus lowering atmospheric CO₂. Nevertheless, no quantitative estimate of this hypothesis has been produced up to now. Here we use an intermediate complexity model including the carbon cycle to evaluate the response of the carbon reservoirs in the atmosphere, ocean and land in response to the changes of orbital forcings, ice sheet configurations and atmospheric CO₂ concentrations over the last nine interglacials. We show that the ocean takes up more carbon during pre-MBE interglacials in agreement with data, but the impact on atmospheric CO₂ is limited to a few parts per million. Terrestrial biosphere is simulated to be less developed in pre-MBE interglacials, which reduces the storage of carbon on land and increases atmospheric CO₂. Accounting for different simulated ice sheet extents modifies the vegetation cover and temperature, and thus the carbon reservoir distribution. Overall, atmospheric CO₂ levels are lower during these pre-MBE simulated interglacials including all these effects, but the magnitude is still far too small. These results suggest a possible misrepresentation of some key processes in the model, such as the magnitude of ocean circulation changes, or the lack of crucial mechanisms or internal feedbacks, such as those related to permafrost, to fully account for the lower atmospheric CO₂ concentrations during pre-MBE interglacials.

1 Introduction

Ice core data have shown that atmospheric CO₂ concentration has been different during interglacials of the last 800 000 years (Lüthi et al., 2008; Bereiter et al., 2015). Older interglacials before the Mid-Brunhes Event (MBE) around 430 ka BP, i.e. Marine Isotope Stages (MIS) 13, 15, 17 and 19, are characterised by relatively lower atmospheric CO₂, around 240 ppm, compared to more recent interglacials, i.e. MIS 1, 5, 7, 9 and 11, which have a higher CO₂ level of around 280 ppm (Fig. 1a).

Proxy data such as the marine benthic foraminifera δ¹⁸O stack record, embedding both deep-sea temperature and ice sheet volume (Lisiecki and Raymo, 2005), indicate that older interglacials (pre-MBE) experienced a colder climate than the more recent ones (post-MBE). This tendency is also supported by individual δ¹⁸O and sea surface temperature (SST, derived from Mg∕Ca paleothermometry, alkenones or foraminifera assemblages) records from marine sediment cores (Lang and Wolff, 2011; Past Interglacials Working Group of PAGES, 2016), although some individual sub-stages such as MIS 7c and 7e were colder than the post-MBE mean interglacial climate.

To explain the different climates of interglacials before and after the MBE, modelling studies have shown that it is necessary to include the change of atmospheric CO₂ (Yin and Berger, 2010, 2012). Indeed, these numerical simulations with an intermediate complexity model have demonstrated that differences in Earth's orbital configuration, and hence seasonal and spatial distribution of insolation, alone cannot explain the colder climate recorded during pre-MBE interglacials, whereby lower atmospheric CO₂ concentration is also necessary to simulate colder climate (Yin and Berger, 2010, 2012). However, the reasons for the lower CO₂ values remain elusive and very few modelling exercises have tackled the issue of different CO₂ levels during interglacials before and after the MBE. Köhler and Fischer (2006) have produced transient simulations of the last 740 000 years using the BICYCLE box model. They used several paleoclimatic records such as ocean temperature, sea ice, sea level, ocean circulation, marine biota, terrestrial biosphere and CaCO₃ chemistry to force their model forward. They run a set of simulations prescribing only one forcing at a time and another with all forcings excluding one at a time, which allows them to analyse which forcings are the most important. In their simulations, they have shown that the lower CO₂ values during pre-MBE are mainly explained by the prescribed lower Southern Ocean sea surface temperature and weaker Atlantic meridional overturning circulation (low North Atlantic Deep Water formation and Southern Ocean vertical mixing) compared to post-MBE interglacials. Using an intermediate complexity model, Yin (2013) conversely simulated vigorous bottom water formation and stronger ventilation in the Southern Ocean during pre-MBE interglacials and suggested this could increase deep oceanic carbon storage and lower atmospheric CO₂. However, this effect on the ocean carbon reservoir and atmospheric CO₂ has not been evaluated yet in a climate model including a carbon cycle representation.

Changes in surface temperature also modify the partition of the carbon cycle: in the ocean, colder SST increases the solubility of CO₂, increasing its potential uptake from the atmosphere during pre-MBE interglacials. In contrast, on land a colder climate might yield a decrease in biomass, reducing CO₂ uptake via lower continental carbon storage. Because the ice sheets in the Northern Hemisphere are different during the interglacials in response to the different values of CO₂ and orbital configurations (Ganopolski and Calov, 2011), they might also have an impact on the carbon cycle, for example by modifying the terrestrial biosphere extent.

Here, we test the impact of the different orbital configurations of the last nine interglacials on the carbon cycle. For this purpose, we use a coupled carbon–climate model to evaluate the changes of carbon storage in the ocean and in the terrestrial biosphere, as well as the impact of different Northern Hemisphere ice sheet volumes.

Table 1 Dates of orbital parameters and CO₂ used for the simulations (Lüthi et al., 2008), and sea level anomalies as compared to present-day conditions (m) corresponding to the prescribed ice sheets (Ganopolski and Calov, 2011).

2 Methods

We use the iLOVECLIM climate model of intermediate complexity, which is a new development branch (code fork) of the LOVECLIM model in its version 1.2 (Goosse et al., 2010). iLOVECLIM has an atmosphere module (ECBILT) with a T21 spectral grid truncation (∼ 5.6^∘ in latitude–longitude in the physical space) and three vertical layers. The ocean component (CLIO) has a horizontal resolution of 3^∘ by 3^∘ and 20 vertical levels. The evolution of the terrestrial biosphere, i.e. the proportion of desert, grasses and tree cover, is computed by the VECODE model (Brovkin et al., 1997). It includes a carbon cycle module on land and in the ocean (Bouttes et al., 2015). iLOVECLIM is an evolution from the LOVECLIM model used in previous model studies of the last nine interglacials focused on climate (Yin and Berger, 2010, 2012; Yin, 2013). It has the same atmospheric and oceanic modules but includes a different carbon cycle representation in the ocean (Bouttes et al., 2015). We have chosen the same dates for the nine orbital configurations as in Yin and Berger (2010, 2012) and Yin (2013), i.e. the maximum of insolation preceding the δ¹⁸O peak values (Table 1, Fig. 1b and c). Contrary to most simulations from these studies, we also use the CO₂ values (as well as CH₄ and N₂O) at the same dates as for the orbital configurations (and not at the CO₂ peak), but as stated in Yin and Berger (2012), this may not affect the main results concerning the simulated climatic changes.

In the model, we separate the atmospheric CO₂ concentration into two distinct variables depending on its physical and chemical impact. The first one is used in the radiative scheme of the atmosphere, for which we prescribe the CO₂ from measured values in all the described simulations (Lüthi et al., 2008; Fig. 1a). Another atmospheric CO₂ is computed interactively in the model as a result of the balance of the carbon fluxes between the different carbon subcomponents (atmosphere, ocean and terrestrial biosphere). We make this choice of keeping the two separated to ensure that the climate simulated by the model is coherent with past measured atmospheric CO₂. In other words, we consider the atmospheric CO₂ concentration as an imposed external forcing, while within the carbon cycle the atmospheric CO₂ concentration is allowed to vary but does not impact the atmospheric radiative forcing. By doing this, we limit the number of degrees of freedom in our climate carbon system, which notably allows us to avoid the complication arising from simulating a different climate when the climate and carbon are fully coupled.

The simulations performed are snapshots, run with constant orbital and atmospheric CO₂ concentration forcing and integrated over 3000 years, allowing the ocean to reach a quasi-equilibrium. All simulations start from the pre-industrial control one, and the average of the last 100 years is used to analyse the results.

Figure 2 Ice sheet elevation (m) in the Northern Hemisphere simulated by the CLIMBER-2 model (Ganopolski and Calov, 2011) and used in the OVIC series for each interglacial simulation, in anomalies with respect to the pre-industrial elevation.

Table 2 Summary of the three series of simulations.

Our strategy is to evaluate the impact of the different climate and carbon compartments to set the atmospheric CO₂ concentration. For this purpose, we consider three series of simulations, which have all been run for the nine interglacials (Table 2). The first series (“Ocean Carbon”, OC) has fixed ice sheets set to the observed pre-industrial ones and a fixed terrestrial biosphere set to the simulated pre-industrial one. This first set of simulations thus provides the response of the ocean alone to the different orbital parameters and CO₂ levels of the nine interglacials. The second series (“Ocean Vegetation Carbon”, OVC) still has fixed ice sheets but includes an interactive terrestrial biosphere, computed by the model. It gives the response of both the ocean and land vegetation reservoirs to the different orbital parameters and CO₂ as well as their interactions for setting the atmospheric CO₂ concentration. Finally, the third series (“Ocean Vegetation Ice sheet Carbon”, OVIC) has different prescribed ice sheets in the Northern Hemisphere for the nine interglacials. The ice sheet distribution change is based on modelling results, given that the uncertainty from data is very large for the interglacials of the last 800 000 years, especially the oldest ones. The prescribed ice sheet distributions are thus taken from an ice sheet simulation of the last 800 000 years (Ganopolski and Calov, 2011) using the intermediate complexity model CLIMBER-2 (Petoukhov et al., 2000; Ganopolski et al., 2001; Brovkin et al., 2002), including a 3-D polythermal ice sheet model (Greve, 1997). This ice sheet model is coupled to the climate component via surface energy and mass balance interface (Calov et al., 2005), which accounts for the effect of aeolian dust deposition on snow albedo. The ice sheet distribution is chosen 2000 years after the chosen interglacial date to account for the long timescale of the ice sheet response during a deglaciation and to ensure that the ice sheet corresponds to an interglacial configuration. The ice sheet elevations for the nine interglacial simulations are shown in Fig. 2 and the corresponding sea level change in Table 1. The terrestrial biosphere is also interactive in this OVIC series of simulations. This last set of simulations thus adds the effect of having different ice sheets in the Northern Hemisphere for the carbon cycle variations.

Figure 3 Annual SST (^∘C) in (a) the pre-industrial control simulation and (b–j) the interglacial simulations of the OC series with fixed vegetation and fixed ice sheets, in anomalies with respect to the pre-industrial control simulation.

Figure 4 (a, b) Annual SST difference (^∘C), (c, d) meridional overturning circulation difference (Sv) and (e, f) dissolved inorganic carbon difference (µmol kg⁻¹) between pre-MBE (MIS 13, 15, 17, 19) and post-MBE (MIS 5, 7, 9, 11) interglacial simulations for (a, c, e) the OC series with fixed vegetation and fixed ice sheets and (b, d, f) the OVC series with interactive vegetation and fixed ice sheets. The vertical black line indicates the limit between the Southern Ocean south of 32^∘ S and the Atlantic Ocean north of 32^∘ S. The dots in panel (a) are SST data differences based on the Past Interglacials Working Group of PAGES (2016) (Table 3).

Table 3 SST data (^∘C) from the Past Interglacials Working Group of PAGES (2016) and shown in Figs. 4a and b and 14.

3 Results and discussion

3.1 Role of the ocean (OC simulations)

Similar to previous numerical studies of the interglacials with the LOVECLIM model (Yin and Berger, 2010, 2012), the changes in orbital configuration and atmospheric CO₂ altered SSTs and oceanic circulation for each interglacial simulation of the OC series. All simulations have warmer SSTs than the control pre-industrial in the Northern Hemisphere high latitudes (Fig. 3). Except for MIS 1, the SSTs in the post-MBE simulations (corresponding to MIS 5, 7, 9 and 11) are also warmer than in the pre-industrial control in large areas in the Northern Hemisphere mid-latitudes. In the MIS 5, 9 and 11 simulations the SSTs are slightly warmer in the Southern Ocean. In the pre-MBE simulations (MIS 13, 15, 17 and 19), the ocean is mainly colder than the pre-industrial ocean, especially in the Southern Hemisphere. To compare the pre-MBE to post-MBE simulations, we built a composite (average) for each period (pre- and post-MBE). We have excluded MIS 1 from the post-MBE composite, for which the date chosen corresponds to a CO₂ much lower than the other post-MBE interglacials. We thus consider MIS 5, 7, 9 and 11 in the post-MBE composite and MIS 13, 15, 17 and 19 in the pre-MBE composite. The difference between the pre- and post-MBE composites shows colder SSTs in the pre-MBE interglacial simulations compared to the post-MBE simulations, especially in the Southern Ocean (Fig. 4a). The difference in global mean simulated SST is −0.30 ^∘C, reaching up to −0.36 ^∘C in the North Atlantic (between 30 and 65^∘ N) and −0.43 ^∘C in the Southern Ocean (south of 54^∘ S). This is in general agreement with SST data, which indicate colder SSTs in the pre-MBE interglacial oceans, especially in the Southern Ocean (Table 3, Fig. 4a), although the comparison is limited by the lack of SST records across the MBE.

Figure 5 Meridional overturning circulation (Sv) in the Southern Ocean and in the Atlantic Ocean north of 32^∘ S in (a) the pre-industrial control simulation and (b–j) the interglacial simulations of the OC series with fixed vegetation and fixed ice sheets, in anomalies with respect to the pre-industrial control simulation. The vertical black line indicates the limit between the Southern Ocean south of 32^∘ S and the Atlantic Ocean north of 32^∘ S.

Compared to the pre-industrial period, the ventilation of the Southern Ocean is increased in all simulations (Fig. 5). The formation of Antarctic Bottom Water (AABW) as well as the wind-driven meridional cell between 40 and 60^∘ S (so-called Deacon cell) are both stronger. On average, the maximum of the Deacon cell is increased by 7 % between pre- and post-MBE simulations, while AABW is increased by 18 % (Fig. 4c). The meridional overturning circulation is also slightly increased by 6 % and deepened in the Atlantic Ocean. All these results concerning oceanic circulation changes are very similar to those of Yin et al. (2013), allowing us to test their hypothesis on the impact of these changes on ocean carbon uptake and atmospheric CO₂ concentrations.

The changes in global ocean circulation and SST modify the carbon storage in the ocean. The colder SST, which increases dilution of CO₂ at the ocean surface, associated with stronger ventilation in the pre-MBE simulations yields a larger carbon uptake by the ocean. Such processes result in higher dissolved inorganic carbon (DIC) concentrations in the Southern Ocean, as well as higher DIC concentration in the upper ocean (first 2 km of the ocean) (Fig. 4e), reflecting the average increase of 4.7 Gt C in pre-MBE simulations. Conversely, the DIC slightly decreases in the deeper ocean, which may be due to the increased ventilation of North Atlantic Deep Water bringing more carbon back from the deep ocean to the surface.

Figure 6 Simulated CO₂ in the interglacial simulations of the OC series as a function of the measured CO₂ from data (Lüthi et al., 2008). The Pearson correlation coefficient and the p value are indicated at the top.

The stronger uptake of carbon by the ocean leads to a decrease in atmospheric CO₂ concentration during the pre-MBE interglacials compared to the post-MBE ones (Fig. 6), in agreement with CO₂ data as shown by the very good correlation between the measured and simulated values (r = 0.91, p < 0.01; Fig. 6). However, the difference in magnitude between pre- and post-MBE values is only a few parts per million in the simulations. Thus, even though simulations qualitatively reproduce the geological CO₂ trend, with lower values during the pre-MBE than the post-MBE interglacials, the magnitude of the difference is much lower in the simulation (∼ 1–5 ppm) than in the data (∼ 30–40 ppm). In fact, the slope of the linear regression between simulated and observed atmospheric CO₂ concentration is 0.07, indicative of a ∼ 14-fold underestimation by the simulations.

Figure 7 Tree cover (%) change with respect to the pre-industrial control simulation for the OVC series with interactive vegetation and fixed ice sheets.

Figure 8 Annual SST difference (^∘C) between simulations with interactive vegetation (OVC) and with fixed vegetation (OC).

Figure 9 Meridional overturning circulation difference (Sv) between simulations with interactive vegetation (OVC) and with fixed vegetation (OC). The vertical black line indicates the limit between the Southern Ocean south of 32^∘ S and the Atlantic Ocean north of 32^∘ S.

Hence the ocean carbon uptake in the simulations is not sufficient to drive a significant lowering of atmospheric CO₂. Either the change in global ocean circulation and SST should be larger or another mechanism and feedbacks need to be taken into account to modify the biological or physical carbon uptake and amplify the initial change. Since the representation of bottom water formation in the Southern Ocean is biased in the model, with an over-representation of open ocean convection, as is also the case for many more complex general circulation models (Heuzé et al., 2013), it is possible that this hinders simulating the full range of carbon storage due to ocean circulation changes, as it is suspected for colder periods such as the Last Glacial Maximum (around 21 000 years ago) (Fischer et al., 2010).

3.2 Role of land vegetation and soils (OVC simulations)

In the first series of simulations, solely the ocean was allowed to respond to the different external forcings while land vegetation and soils were fixed to their pre-industrial distribution. To account for changes in land vegetation and soils on the carbon cycle, a second series of simulations (OVC) was conducted with an interactive terrestrial biosphere module on top of the ocean's module (Table 2).

Compared to the pre-industrial control simulation, more trees develop in North Africa and the southern part of Eurasia in these interglacial simulations, while the tree cover is reduced in central North America and some regions in the northern part of Eurasia (Fig. 7). When compared to OC simulations with fixed vegetation, the OVC simulations demonstrate a surface ocean warming almost everywhere except in the North Atlantic for some interglacials (Fig. 8). In response to the global warmer surface ocean, the stratification in the convection region in the North Atlantic increases (e.g. Swingedouw et al., 2007a) leading to a slowdown of the Atlantic meridional oceanic circulation, especially for MIS 1, 5 7, 9 and 15 (Fig. 9).

Figure 10 Carbon stocks (Gt C) in the three reservoirs (atmosphere, ocean and land) for each simulation. (a) OC series with fixed vegetation and fixed ice sheets, (b) OVC series with interactive vegetation and fixed ice sheets, and (c) OVIC series with interactive vegetation and different prescribed ice sheets. The stocks are given as anomalies with respect to the control pre-industrial simulation.

Figure 11 (a) CO₂ concentration (ppm) at the end of the simulations and (b) composite (average) CO₂ (ppm) in the pre-MBE (MIS 13, 15, 17, 19) and post-MBE (MIS 5, 7, 9, 11) interglacial simulations.

For the carbon cycle, the activation of the terrestrial module results in more carbon stored in land vegetation and soils for all interglacial simulations (Fig. 10b) since the vegetation cover increases compared to the control because of the warmer climate. This tends to lower atmospheric CO₂ concentration, hence the pCO₂ difference at the air–sea interface, leading to an outgassing of carbon from the ocean to the atmosphere, which ultimately decreases the storage of carbon in the ocean. The ocean carbon storage is also diminished compared to the series of simulations with fixed vegetation due to the warmer ocean temperature, which reduces the CO₂ solubility in water. The increase in carbon storage in the terrestrial biosphere is generally larger than the loss of carbon from the ocean so that the carbon content of the atmosphere is also diminished in these simulations compared to the fixed vegetation simulations, and atmospheric CO₂ is slightly lower or not changed (Fig. 11a).

Figure 12 (a, b) Tree cover (%) and (c, d) carbon storage (kg C m⁻²) difference between pre-MBE (MIS 13, 15, 17, 19) and post-MBE (MIS 5, 7, 9, 11) interglacial simulations for (a, c) the OVC series with interactive vegetation and fixed ice sheets and (b, d) the OVIC series with interactive vegetation and different ice sheets. Qualitative indication of tree cover change from data is indicated with dots: blue indicates a reduction of tree cover on average during pre-MBE interglacials compared to post-MBE interglacials, and red indicates an increase.

In terms of difference between pre- and post-MBE interglacial simulations, we find less vegetation cover in most areas (except in North Africa and parts of south Eurasia) and consequently less carbon stored (−48 Gt C) in land vegetation and soils in the pre-MBE simulations compared to post-MBE simulations (Fig. 12). This effect tends to increase atmospheric CO₂ on average in pre-MBE interglacial simulations.

For the ocean, the differences between pre- and post-MBE simulations are similar to the ones for the simulations with fixed vegetation (OC). On average, the SST is lower in the pre-MBE simulations compared to post-MBE simulations (−0.28 ^∘C globally, −0.31 ^∘C in the North Atlantic and −0.44 ^∘C in the Southern Ocean) except in a small area in the North Atlantic (Fig. 4b), and the ventilation is increased in the pre-MBE simulations (Fig. 4d). Hence the ocean can store more carbon in the pre-MBE simulations than the post-MBE simulations with an average increase in carbon storage in the pre-MBE ocean of 43 Gt C compared to the post-MBE ocean. Similarly, the DIC concentration is higher in pre-MBE simulations, especially in the upper ocean and deep Southern Ocean, as in the previous series of simulations with fixed vegetation (Fig. 4f).

As the diminution in land carbon storage is larger than the increase in ocean carbon storage, more carbon is conserved in the atmosphere resulting in higher CO₂ on average for the pre-MBE simulations than the post-MBE simulations (Fig. 11b). There is thus a qualitative disagreement (negative correlation of −0.33 (p = 0.38) between simulated and observed atmospheric CO₂ for the interglacials considered) with the observations.

Nevertheless, it should be noted that permafrost (frozen soil) was not taken into account in these simulations. If there was more permafrost during the colder pre-MBE interglacials, it could store more carbon on land and counteract the loss of carbon due to the reduction of vegetation cover and production (Crichton et al., 2016).

Table 4 Tree cover (%) reconstructed from pollen data at three sites. SW Iberian margin: MIS 5 (MD95-2042, Sanchez Goñi et al., 1999), MIS 7 (MD01-2443; Roucoux et al., 2006), MIS 9, 13, 15, 17 (IODP 1385, unpublished data), MIS 11 (U1385, Oliveira et al., 2016), MIS 19 (U1385, Sanchez Goñi et al., 2016); Tenaghi Philippon, Greece (Past Interglacials Working Group of PAGES, 2016); and Funza, Colombia (Past Interglacials Working Group of PAGES, 2016).

Comparison with pollen data (Table 4) indicates that the model is in qualitative agreement with reconstructed tree cover change in South America where the tree cover was smaller on average in pre-MBE than in post-MBE interglacials (Fig. 12a). In southern Europe, the tree fraction data indicate that on average slightly more tree cover prevailed during pre-MBE than during post-MBE interglacials, also in agreement with simulations, although the variability in the data among interglacials is large.

3.3 Impact of different ice sheets (OVIC simulations)

The last series of simulations (OVIC) has the same design as OVC but also takes into account possible differences in ice sheet distribution in the Northern Hemisphere based on numerical simulations (Ganopolski and Calov, 2011). The simulated ice sheet distributions used for our interglacial simulations mainly differ in North America. On average, the North American ice sheet is more extended in the pre-MBE interglacials compared to the post-MBE interglacials (Fig. 13).

The change of ice sheet extent has a large regional impact on vegetation cover, which is reduced where the ice sheet extends more. On average, it results in a reduction of vegetation in North America in the pre-MBE interglacials, when the ice sheet is more extended, compared to the post-MBE interglacials (Fig. 12b and d).

Figure 13 Ice sheet elevation difference (m) between the average of the pre-MBE (MIS 13, 15, 17, 19) and post-MBE (MIS 5, 7, 9, 11) interglacial simulations.

Figure 14 Annual sea surface temperature difference (^∘C) between the average of the pre-MBE (MIS 13, 15, 17, 19) and post-MBE (MIS 5, 7, 9, 11) interglacials with interactive vegetation and different ice sheets (OVIC). The vertical black line indicates the limit between the Southern Ocean south of 32^∘ S and the Atlantic Ocean north of 32^∘ S. The dots are SST data differences based on the Past Interglacials Working Group of PAGES (2016) (Table 3).

The increase in ice sheet extent and diminution of vegetation cover for pre-MBE simulations has two main impacts for the carbon cycle: (i) it diminishes the terrestrial biosphere carbon storage, increasing atmospheric CO₂ and (ii) it also cools global climate due to the high ice albedo. The SST also decreases by 0.32 ^∘C on the global scale. Sea surface temperature changes are especially pronounced in polar zones, with a drop of 0.49 ^∘C in the North Atlantic and 0.47 ^∘C in the Southern Ocean (Fig. 14 compared to Fig. 4b). Consequently, the ocean carbon storage increases, which lowers atmospheric CO₂. This second effect dominates and the overall result is lower atmospheric CO₂ concentration in pre-MBE simulations compared to post-MBE simulations (Fig. 11). As for the other processes analysed, they only modifies atmospheric CO₂ by a few parts per million, though correcting back (compared to OVC simulations) the difference pre-MBE minus post-MBE towards the observations. Nevertheless, the correlation between simulated and measured CO₂ (accounting for MIS1) is very small (−0.08) and not significant (p = 0.83). The magnitude of the changes of atmospheric CO₂ among the different interglacials is once again largely underestimated as compared to observations.

Accounting for different ice sheets in the OVIC series seems to improve the model–data comparison in southern Europe for tree cover (Fig. 12b) where the data are at the limit between regions of more tree coverage and less tree coverage in the model. This highlights the role of ice sheet extent in setting the vegetation pattern. Nevertheless, the uncertainty in ice sheet distribution is very large and the model-based reconstruction might not be accurate. For example, the lack of ice-rafted debris from North America before MIS16 and the presence of ice-rafted debris from Europe indicate that the ice sheet over Europe (Hoddell et al., 2008) could have been more extended and the Laurentide ice sheets in North America were not. In addition, the model-based reconstruction that we used shows relatively small changes of sea level equivalent between interglacials. Data reconstructions seem to indicate possible larger differences between interglacials (Spratt and Lisiecki, 2016), whose effect on the size of the land surface and the carbon cycle remains to be tested. Sensitivity experiments with prescribed idealised ice sheets designed to be very different would help to evaluate their impact.

4 Conclusions

Using a fully coupled climate model of intermediate complexity including an interactive carbon cycle, we have shown that the difference between pre-MBE and post-MBE cannot be explained by the simulated changes in ocean and vegetation induced by orbital and greenhouse gases forcing. While the oceanic response alone is in qualitative agreement with data (sign of the changes, correlation between each interglacial), it largely underestimates the amplitude of the changes. Past work suggested that either a vigorous AABW (Yin, 2013) or weak Atlantic thermohaline circulation (Köhler and Fischer, 2006) during pre-MBE interglacials could increase the oceanic carbon storage and explain the lower CO₂ than during post-MBE interglacials. Other studies for different background climates have shown opposite results with respect to the effect of ocean circulation on carbon storage. A weaker Atlantic meridional oceanic circulation could either result in more ocean carbon storage with a pre-industrial climate (Obata, 2007; Menviel et al., 2008; Bozbiyik et al., 2011) or glacial climate (Menviel et al., 2008), or it could yield less ocean carbon storage with a pre-industrial climate (Marchal et al., 1998; Swingedouw et al., 2007b; Bouttes et al., 2012) or a glacial climate (Schmittner and Galbraith, 2008; Bouttes et al., 2012; Schmittner and Lund, 2015). Menviel et al. (2014) showed that on top of changes in North Atlantic Deep Water formation, modifications of AABW and North Pacific Deep Water (NPDW) formations could results in different oceanic carbon storage. Data indicate that the modern reduction of carbon uptake in the North Atlantic is due to a reduction in the overturning circulation (Perez et al., 2013). Because the atmospheric CO₂ change that we simulate has a low magnitude of only a few parts per million, it is not yet possible to infer whether stronger or weaker overturning during pre-MBE interglacials could have significantly lowered atmospheric CO₂.

Furthermore, accounting for the vegetation response complicates the simulated response and entirely removes the qualitative agreement. The vegetation response depends on ice sheet extent, and accounting for ice sheet variations limits the disagreement. Comparison of simulated vegetation changes with available pollen data indicates partial agreement, underlying the need to improve vegetation simulations and increase the data coverage to more precisely constrain the change of vegetation cover. The vegetation model in iLOVECLIM only simulates grass and trees; to better evaluate the different vegetation responses to orbital and CO₂ forcings it would be useful to use a more complex terrestrial biosphere model.

We argue that additional processes need to be accounted for or should be better represented in climate models to explain the observations. It is possible that either many different processes, some of them not included in the present model, lead to the observed atmospheric CO₂ concentration or that just a first-order process is misrepresented or not included. In particular, the storage of carbon in frozen soils (permafrost) should be included in future modelling work. Response of the Southern Ocean, the widest oceanic region with large air–sea fluxes of CO₂, is also a good candidate, given the known deficiency in coarse-resolution climate models for the representation of the key element of its dynamics (eddies, katabatic winds, AABW formation, brines, etc.). The impact of ventilation changes could be tested by artificially modifying the buoyancy forcing in the areas of bottom water formation. The use of higher-resolution models in this region could help to better evaluate its response to different interglacial conditions.