Articles | Volume 17, issue 1
Clim. Past, 17, 361–377, 2021

Special issue: Oldest Ice: finding and interpreting climate proxies in ice...

Clim. Past, 17, 361–377, 2021

Research article 01 Feb 2021

Research article | 01 Feb 2021

Reconstructing the evolution of ice sheets, sea level, and atmospheric CO2 during the past 3.6 million years

Reconstructing the evolution of ice sheets, sea level, and atmospheric CO2 during the past 3.6 million years
Constantijn J. Berends1, Bas de Boer2, and Roderik S. W. van de Wal1,3 Constantijn J. Berends et al.
  • 1Institute for Marine and Atmospheric research Utrecht, Utrecht University, Utrecht, the Netherlands
  • 2Earth and Climate Cluster, Faculty of Science, Vrije Universiteit Amsterdam, Amsterdam, the Netherlands
  • 3Faculty of Geosciences, Department of Physical Geography, Utrecht University, Utrecht, the Netherlands

Correspondence: Constantijn J. Berends (


Understanding the evolution of, and the interactions between, ice sheets and the global climate over geological timescales is important for being able to project their future evolution. However, direct observational evidence of past CO2 concentrations, and the implied radiative forcing, only exists for the past 800 000 years. Records of benthic δ18O date back millions of years but contain signals from both land ice volume and ocean temperature. In recent years, inverse forward modelling has been developed as a method to disentangle these two signals, resulting in mutually consistent reconstructions of ice volume, temperature, and CO2. We use this approach to force a hybrid ice-sheet–climate model with a benthic δ18O stack, reconstructing the evolution of the ice sheets, global mean sea level, and atmospheric CO2 during the late Pliocene and the Pleistocene, from 3.6 million years (Myr) ago to the present day. During the warmer-than-present climates of the late Pliocene, reconstructed CO2 varies widely, from 320–440 ppmv for warm periods to 235–250 ppmv for the early glacial excursion ∼3.3 million years ago. Sea level is relatively stable during this period, with maxima of 6–14 m and minima of 12–26 m during glacial episodes. Both CO2 and sea level are within the wide ranges of values covered by available proxy data for this period. Our results for the Pleistocene agree well with the ice-core CO2 record, as well as with different available sea-level proxy data. For the Early Pleistocene, 2.6–1.2 Myr ago, we simulate 40 kyr glacial cycles, with interglacial CO2 decreasing from 280–300 ppmv at the beginning of the Pleistocene to 250–280 ppmv just before the Mid-Pleistocene Transition (MPT). Peak glacial CO2 decreases from 220–250 to 205–225 ppmv during this period. After the MPT, when the glacial cycles change from 40 to 80 120 kyr cyclicity, the glacial–interglacial contrast increases, with interglacial CO2 varying between 250–320 ppmv and peak glacial values decreasing to 170–210 ppmv.

1 Introduction

Understanding the response of ice sheets and the global climate as a whole to changes in the concentrations of atmospheric CO2 is important for understanding the future evolution of the climate system. Since large-scale changes in ice-sheet geometry typically occur over thousands to tens of thousands of years, sources of information other than direct observational evidence are required. In order to gain more insight in the relation between these components of the Earth system, studying their evolution during the geological past is useful.

One particularly rich source of information is presented by ice cores. The chemical and isotopical content of the ice itself can provide valuable information on the state of the global climate, and in particular on the temperature, at the time the ice was formed (Dansgaard, 1964; Jouzel et al., 1997; Alley, 2000; Kindler et al., 2014). Air bubbles, trapped when the snow compresses first into firn and ultimately into ice, contain tiny samples of the atmosphere from the time the air got trapped. The oldest ice core presently available, the EPICA Dome C core, contains ice, air bubbles, and hydrates dating back 800 kyr (Bereiter et al., 2015). The information obtained from these ice cores has greatly improved our understanding of the dynamics of the glacial cycles of the Late Pleistocene.

Different methods of relating chemical and isotopic properties of ocean sediments to the atmospheric CO2 concentration have been used to create proxy data extending back further in time than the ice-core record. Many studies have measured δ11B of different fossil foraminifera, using the observed relation to seawater pH to calculate atmospheric CO2 concentrations (Hönisch et al., 2009; Seki et al., 2010; Bartoli et al., 2011; Martínez-Botí et al., 2015; Foster and Rae, 2016; Stap et al., 2016; Chalk et al., 2017; Dyez et al., 2018; Sosdian et al., 2018). Another line of evidence has focused on the δ13C of alkenones in fossil foraminifera (Seki et al., 2010; Badger et al., 2013; Zhang et al., 2013), although the reliability of this proxy for CO2 concentrations lower than ∼400ppmv has recently been called into question (Badger et al., 2019). A different line of work relates the density of stomata on fossil plant leaves to atmospheric CO2 concentration, providing data throughout the Cenozoic (Kürschner et al., 1996; Wagner et al., 2002; Finsinger and Wagner-Cremer, 2009; Beerling and Royer, 2011; Stults et al., 2011; Bai et al., 2015; Hu et al., 2015). This proxy has been the subject of discussion regarding its reliability (Indermühle et al., 1999; Jordan, 2011; Porter et al., 2019), based on its discrepancies with the ice-core record and presently unresolved problems due to the influence of other effects such as evolution, extinction, changes in local environment, and immigration of species, although recent work has gone some way towards resolving these issues (Franks et al., 2014).

Information about the global glacial state is furthermore obtained from the δ18O of fossil benthic foraminifera, which is influenced by both total ice volume and deep-sea temperature. Ocean sediment cores containing foraminiferal shells have been used to create stacks of benthic δ18O records dating back 65 Myr (million years; Lisiecki and Raymo, 2005; Zachos et al., 2001, 2008). Figure 1 shows the LR04 stack of benthic δ18O records (Lisiecki and Raymo, 2005), together with the EPICA Dome C CO2 record (Bereiter et al., 2015) and the Earth's orbital parameters and Northern Hemisphere summer insolation (Laskar et al., 2004), during the last 800 kyr.

Figure 1Earth's orbital parameters (obliquity, eccentricity, and precession) and the resulting Northern Hemisphere summer insolation (Laskar et al., 2004), benthic δ18O (Lisiecki and Raymo, 2005), and ice-core CO2 (Bereiter et al., 2015) during the past 800 kyr.


Since benthic δ18O contains both an ice-based and a climate-based signal, using it to understand the evolution of either of these two Earth system components is only possible after disentangling the two contributions. Several studies have aimed at performing this separation by deriving either of the two signals from other independent proxies. One approach has been to derive deep-water temperatures from foraminiferal MgCa ratios (Sosdian and Rosenthal, 2009; Elderfield et al., 2012; Shakun et al., 2015). Global mean sea level is then solved as a closure term. Alternatively, temperature proxies have been obtained from records of oxygen and hydrogen isotope abundances from ice cores from Greenland (Alley, 2000) and Antarctica (Jouzel et al., 2007). Focusing on the signal from ice volume rather than temperature, Rohling et al. (2014) used planktic δ18O in the Mediterranean Sea as a proxy for sea level at the Strait of Gibraltar, translating the result into a global mean sea-level record spanning the last 5.5 Myr. Grant et al. (2014) applied the same method to the Red Sea, producing a record of sea level at the Bab el-Mandeb Strait with a higher accuracy, but going back only 500 kyr.

These studies are generally viewed as “data” studies; observed variables (chemical concentrations, isotope ratios, etc.) are used to derive climate parameters (ice volume, global mean temperature, pCO2, etc.) using relatively simple time-independent concepts. Discussions about the validity of the results generally focus less on the physical processes described by these simple models and more on the properties of the data, such as sample contamination and diagenesis, statistical limitations, and measurement uncertainty.

A different family are the “model” studies, where physics-based models are used to describe the evolution of (components of) the Earth system (atmosphere, ocean, cryosphere, carbon cycle, etc.) through time. Typically, the aim of a model study is to determine if our understanding of a physical process is good enough to explain the observations, or to use that understanding to make predictions of that process in the future. Several recent studies have aimed to reproduce the evolution of global climate and the cryosphere throughout the Pleistocene, using primarily insolation as forcing (Brovkin et al., 2012; Willeit et al., 2015, 2019). By studying the differences between model results and observations from proxy data, such studies can help with interpreting these data, providing insight into the physical processes that govern the relation between observed and derived variables. In turn, the differences between data and model results can help in assessing the importance of physical processes that are not included in the model.

Inverse modelling is a hybrid method, using an approach that combines elements from both these families of research. This method aims to derive the evolution of the entire global climate–cryosphere system through time, based on observations of benthic δ18O (Bintanja et al., 2005; Bintanja and van de Wal., 2008; de Boer et al., 2010, 2012, 2013, 2014, 2017; van de Wal et al., 2011; Stap et al., 2016, 2017; Berends et al., 2018a, 2019). This is done by using ice-sheet models, with climate forcing components varying in complexity from simple parameterizations to fully coupled general circulation models (GCMs), to reconstruct ice-sheet evolution. Such models can be used to calculate the contributions to the observed benthic δ18O signal from both ice volume and deep-water temperature over time. By using a tool called an “inverse routine”, the model can be forced to (almost) exactly reproduce the benthic δ18O signal, providing a reconstruction of global climate as it should have evolved in order to explain the observations. The advantage of this approach over other proxy-based reconstructions is that the simulated changes in global climate and ice-sheet evolution are mutually consistent, building on the physical equations within the model framework. Earlier studies adopting this approach used simple one-dimensional ice-sheet models to represent all ice on Earth or on a single hemisphere, and different parameterizations of the relation between ice volume and climate (e.g. de Boer et al., 2010; Stap et al., 2017). More recent studies have used more elaborate 3-D ice-sheet models covering different regions of the Earth, using more comprehensive mass balance parameterizations and representations of the global climate (e.g., Bintanja and van de Wal, 2008; de Boer et al., 2013, 2017; Berends et al., 2018a, 2019). The most recent of these studies, by Berends et al. (2019), used a 3-D ice-sheet–shelf model, forced with output from several different GCM simulations by a so-called “matrix method” of model coupling, to simulate the last four glacial cycles. They showed that their results agreed with geomorphological and proxy-based evidence of ice-sheet volume and extent, benthic δ18O, deep-water temperature, ice-sheet temperature, and atmospheric CO2.

The work presented here builds on the work by Berends et al. (2019), extending their results to 3.6 Myr ago to produce a time-continuous, self-consistent reconstruction of atmospheric CO2, global climate, ice-sheet geometry, and global mean sea level over the period where ice-core records of CO2 are not available. The inverse modelling approach that was adopted to force this model with a benthic δ18O record is described in Sect. 2.1. The ice-sheet and climate models are described in Sects. 2.2 and 2.3, respectively, while the matrix method used to couple these two models is described in Sect. 2.4. The resulting reconstructions of CO2, ice volume, temperature, and sea level over the past 3.6 Myr are presented in Sect. 3 and discussed in Sect. 4.

2 Methodology

2.1 Inverse modelling

In order to disentangle the contributions to the benthic δ18O signal from global land ice volume and deep-water temperature changes, we use the inverse modelling method proposed by Bintanja and van de Wal (2008) and refined by de Boer et al. (2013, 2014, 2017) and Berends et al. (2019). The two contributions to the benthic δ18O signal are not independent; both result from changes in the Earth's climate, driven by changes in insolation and CO2. A cooling of the global climate will, over time, affect the temperature in the deep ocean, especially when the cooling occurs at high latitudes, where deep water formation occurs. Simultaneously, such a global cooling leads to an increase in ice volume, which affects the δ18O of the seawater itself. The inverse modelling method is a tool that determines how modelled CO2 should have evolved over time in order to affect the global climate in such a way that the observed benthic δ18O record is reproduced. The modelled value for CO2 is used, together with the ice sheets simulated by the ice-sheet model, to determine the global climate by using a matrix method of model forcing. This method was presented by Berends et al. (2018a) and refined by Berends et al. (2019), and is described in more detail in Sect. 2.4. The resulting climate is used to determine the surface mass balance over the ice sheets, which forces the ice-sheet model. The resulting global ice volume and deep-water temperature anomaly (which is derived from the high-latitude annual mean surface temperature anomaly, using a moving average of 3000 years) are used to calculate a modelled value of benthic δ18O. This approach is visualized conceptually in Fig. 2.

Figure 2A conceptual visualization of the inverse forward modelling approach. The model is forced externally by an insolation reconstruction and a benthic δ18O record (black boxes). pCO2 is changed over time based on the difference between observed and modelled δ18O. The climate matrix interpolates between the pre-calculated GCM snapshots, based on the prescribed pCO2 value and the modelled state of the cryosphere (ice thickness and albedo), to determine the climate that is used to calculate the surface mass balance over the ice sheets. The modelled ice sheets and climate are then used to calculate a modelled benthic δ18O value, which is used to update modelled CO2 in the next time step. Figure adapted from Berends et al. (2019).

The modelled value δ18Omod is compared to the observed value δ18Oobs at the model time t. A modelled value that is too low indicates that the modelled climate is too warm or that the ice sheets are too small. pCO2 is then decreased in the next time step (with an increment proportional to the δ18O discrepancy). This relationship is quantified by the following equation:

(1) p CO 2 = p CO 2 + α δ 18 O mod - δ 18 O obs .

Here, pCO2 is the mean modelled pCO2 over the preceding 8.5 kyr, and α is a scaling parameter which controls how fast modelled CO2 is allowed to change (serving as a low-pass filter, suppressing overshoot that could result from large changes in the benthic δ18O record). In order to calculate δ18Omod, the contributions from modelled ice volume and deep-sea temperature are calculated separately. The spatially variable isotope content of the individual ice sheets is tracked through time, with the surface isotope balance based on the observed present-day relation between precipitation rates and isotope content according to Zwally and Giovinetto (1997). Benthic δ18O is assumed to be linearly dependent on the global mean deep-water temperature anomaly, which is calculated by temporally smoothing the high-latitude annual mean surface temperature anomaly over the preceding 3 kyr and multiplying with a scaling factor of 0.25. Berends et al. (2019) showed that this yields a deep-sea temperature anomaly of 2–2.5 K at the Last Glacial Maximum (LGM), in general agreement with proxy-based reconstructions (Shakun et al., 2015). The values of 8.5 kyr for the length of the CO2 averaging window, 3 kyr for the deep-water temperature averaging window, and 120 ppmv ‰−1 for the δ18O–CO2 scaling parameter were determined experimentally to accurately reproduce the observed benthic δ18O record.

The combination of this inverse modelling method to reconstruct pCO2 and the matrix method to determine the global climate has been shown to accurately reproduce changes in benthic δ18O and the individual contributions from both global ice volume and deep-water temperature, as well as ice-sheet volume and extent, ice-sheet surface temperature, and pCO2 during the last four glacial cycles (Berends et al., 2019).

2.2 Ice-sheet model

The evolution of the ice sheets is simulated using ANICE, a coupled 3-D ice-sheet–shelf model (de Boer et al., 2013; Berends et al., 2018a). ANICE uses a combination of the shallow-ice approximation (SIA; Morland and Johnson, 1980) for grounded ice and the shallow-shelf approximation (SSA; Morland, 1987) for floating ice to solve the ice mechanical equations. Basal sliding is described by a Coulomb sliding law and solved using the SSA, using the hybrid approach by Bueler and Brown (2009), where basal friction is determined by bedrock elevation. Internal ice temperatures, used to calculate ice viscosity, are calculated using a coupled thermodynamical module. The surface mass balance is parameterized based on monthly mean surface temperature and precipitation, where ablation is calculated using the surface-temperature–albedo–insolation parameterization, as explained in more detail by Berends et al. (2018a). The solution by Laskar et al. (2004) is used to prescribe time- and latitude-dependent insolation at the top of the atmosphere. A combination of the temperature-based formulation by Martin et al. (2011) and the glacial–interglacial parameterization by Pollard and DeConto (2009) is used to calculate sub-shelf melt underneath the Antarctic ice shelves, calibrated by de Boer et al. (2013) to produce realistic present-day Antarctic shelves and grounding lines. A more detailed explanation is provided by de Boer et al. (2013) and references therein. A simple threshold thickness of 200 m is used to describe ice calving, where any shelf ice below this threshold thickness is removed. The model is run on four separate grids simultaneously, covering North America, Eurasia, Greenland, and Antarctica, as shown in Fig. 3. The horizontal resolution is 40 km for Antarctica, North America, and Eurasia, and 20 km for Greenland.

Figure 3The areas of the world covered by the four model domains of ANICE. In the North America and Eurasia domains, Greenland is omitted. Figure adapted from Berends et al. (2018a).

2.3 Climate model

HadCM3 is a coupled atmosphere–ocean GCM (Gordon et al., 2000; Valdes et al., 2017), which has been shown to accurately reproduce the present-day climate heat budget (Gordon et al., 2000). It has been used for future climate projections in the IPCC AR4 (IPCC, 2007), and palaeoclimate reconstructions such as PlioMIP (Haywood and Valdes, 2003; Dolan et al., 2011, 2015; Haywood et al., 2013) and PMIP2 (Braconnot et al., 2007). Atmospheric circulation is calculated at a resolution of 2.5 latitude by 3.75 longitude. The ocean is modelled at a horizontal resolution of 1.25 by 1.25, with 20 vertical layers. We use results from several different steady-state time slice simulations with HadCM3 of different climate states to force our ice-sheet model, using the matrix method explained in Sect. 2.4.

2.4 Matrix method

According to Pollard (2010), a climate matrix is a collection of pre-calculated output data from several steady-state GCM simulations, called “snapshots”. These snapshots differ from each other in one or more key parameters, such as orbital configuration, prescribed atmospheric pCO2, or ice-sheet configuration. Each of these constitutes a separate dimension of the matrix. When using an ice-sheet model to simulate the evolution of an ice sheet over time, the prescribed climate is determined in every model time step by combining the GCM snapshots according to the position of the ice-sheet model state in the climate matrix, which is determined by the modelled values of the parameters describing the snapshots. This approach occupies the middle ground between offline forcing and fully coupled ice-sheet–climate models. The different GCM snapshots contain the key feedback effects of the altitude and albedo on the temperature. In addition, the matrix method captures the effect of ice-sheet geometry on large-scale atmospheric circulation and precipitation. The matrix method creates a spatially variable linear interpolation of these snapshots, providing a first-order approach to the strength of the feedback and the effects of ice-sheet geometry on atmospheric circulation and precipitation.

In this study, we use the matrix method developed by Berends et al. (2018a), where temperature fields from the different climate states are combined based on a prescribed value for pCO2 and on the internally modelled ice sheets. The feedback of the ice sheets on the climate is calculated via the effect on absorbed insolation through changes in surface albedo. This interpolation is carried out separately for all four ice sheets. The altitude–temperature feedback is parameterized by a constant lapse rate derived from the GCM snapshots. Precipitation fields are combined based on changes in surface elevation, reflecting the orographic forcing of precipitation and the plateau desert effect caused by the presence of a large ice sheet. The equations describing these calculations are presented in Appendix A and explained in more detail by Berends et al. (2018a), who demonstrated the viability of this method by simulating the evolution of the North American, Eurasian, Greenland, and Antarctic ice sheets throughout the entire last glacial cycle at the same time. They showed that model results agree well with available data in terms of ice-sheet extent, sea-level contribution, ice-sheet surface temperature, and contribution to benthic δ18O. Here, we apply this matrix method to the climate matrix created by Berends et al. (2019), consisting of 11 pre-calculated GCM snapshots, created with HadCM3. Two of these, produced by Singarayer and Valdes (2010), respectively represent the pre-industrial period (PI) and the LGM. The other nine, produced by Dolan et al. (2015), represent the global climate during Marine Isotope Stage (MIS) M2 (3.3 Myr ago), for four different possible ice-sheet geometries and two different pCO2 concentrations, plus one Pliocene control run. The total set of 11 snapshots allows the climate matrix to disentangle the effects on climate of changes in pCO2 and changes in ice-sheet extent, and provides information on climates that are both colder and warmer than present-day.

3 Experimental set-up and results

Here, we describe a set of simulations from 3.6 Myr ago to the present day. The choice for this starting point aims to include the end of the late Pliocene and the inception of the Pleistocene glacial cycles. The model was initialized with the same PRISM3 ice sheets that were used as boundary conditions in several of the HadCM3 simulations by Dolan et al. (2015). The model was forced with the LR04 stack of benthic δ18O records (Lisiecki and Raymo, 2005). We chose here to perform three simulations: one “default” run, one with the δ18O forcing adjusted upwards by 0.1 , and one with the forcing adjusting downwards by 0.1 . Berends et al. (2019) showed that the contribution to the uncertainty in the reconstructed CO2 and sea level arising from the uncertainty in the forcing is the largest of all model parameters. Since simulating such long periods of time is very computationally demanding, we decided that the added value of repeating the sensitivity analysis for other model parameters by Berends et al. (2019) did not outweigh the computational cost. Margins of error reported here therefore describe only the uncertainty arising from the uncertainty in the benthic δ18O record. The true uncertainty in our reconstruction, which includes contributions from many other parameterizations made in our model, is significantly larger but remains difficult to quantify.

Figure 4Observed benthic δ18O (black; Lisiecki and Raymo, 2005) and reconstructed CO2 (brown) and global mean sea level (blue) for the entire 3.6 Myr simulation period. The present-day values for the three variables (pre-industrial value of 280 ppmv for CO2) are shown by horizontal dashed lines. Shaded areas indicate the uncertainty in the LR04 benthic δ18O stack, and the resulting uncertainty in the reconstructed CO2 and sea level.


The resulting reconstructions of CO2 and global mean sea level, together with the δ18O forcing, are shown in Fig. 4. The 40 kyr glacial cycles of the Early Pleistocene, between 2.6 and 1.2 Myr ago, are clearly visible, changing to 80–120 kyr cycles after the Mid-Pleistocene Transition (MPT). The visibly larger uncertainty for higher CO2 concentrations, occurring during the warm periods of the late Pliocene, is likely due to the relative sparsity of the climate matrix for warmer-than-present states (three snapshots only) compared to the colder-than-present part (seven snapshots), which might lead to a bias towards larger ice sheets in warm worlds.

Figure 5Reconstructed global mean sea level for the entire simulation period (black), compared to reconstructions based on Red Sea δ18O (Grant et al., 2014; red); Mediterranean δ18O (Rohling et al., 2014; yellow); an ice-sheet-model-based inversion of the global benthic δ18O similar to this study (de Boer et al., 2014; blue); an insolation-forced, fully coupled ice-sheet–climate–carbon-cycle model (Willeit et al., 2019; green); a separation of the ice- and temperature-induced δ18O signals based on MgCa ratios (Elderfield et al., 2012; purple); and a direct scaling of benthic δ18O scaled to match results from geological backstripping (Miller et al., 2012; Naish et al., 2009; dark red). The present-day value of 0 is shown by a dashed line.


The modelled sea level is compared to several other reconstructions in Fig. 5. Shown are the reconstructions based on planktic δ18O in the Red Sea (Grant et al., 2014) and in the Mediterranean δ18O (Rohling et al., 2014), which cover the last 500 kyr and 5.5 Myr, respectively. The reconstruction by de Boer et al. (2014) is based on an ice-sheet-model-based decomposition of the benthic δ18O signal, very similar to the work presented here. The main difference is the climate forcing applied to the ice-sheet model, which by de Boer et al. (2014) was described by a simple globally uniform temperature offset, and which was replaced by the climate matrix approach plus inversely modelled CO2 in our work. Despite this improvement in the climate forcing, which resulted in a simulated ice sheet at LGM which agrees better with geomorphological evidence (Berends et al., 2018a), the sea-level reconstructions are virtually indistinguishable. The reconstruction by Willeit et al. (2019) shown in Fig. 5 is based on a fully coupled ice-sheet–climate–carbon-cycle model, forced only with insolation. However, their ice-sheet model only simulated the Northern Hemisphere; the contribution from Antarctica is assumed to be 10 % of that of the northern ice sheets. While this assumption seems to produce good results during the Pleistocene, its validity during the warmer-than-present late Pliocene is questionable. The reconstruction by Elderfield et al. (2012) was made using a MgCa-based reconstruction of sea surface temperature, which was then used to disentangle the ice volume and deep-sea temperature signals in the benthic δ18O record. Lastly, the reconstruction by Naish et al. (2009), adjusted by Miller et al. (2012) to match results from geological backstripping in New Zealand, suggests a relatively stable sea level during the late Pliocene, 3.4–2.58 Myr ago, which agrees well with our results. Rohling et al. (2014), arguing that geological backstripping can provide information on relative changes in sea level but not on absolute values, added a +20 m offset to the results by Miller et al. (2012), which made those results agree better with their own reconstruction from the Mediterranean Sea. However, our own results are reasonably close to those by Miller et al. (2011) over the period between 3.4 and 2.6 Myr ago without requiring such a correction. The strong variability in global mean sea level visible in the results of Rohling et al. (2014) during this period, which suggests the repeated disappearance and reappearance of most of the East Antarctic Ice Sheet, is not visible in the data of Miller et al. (2012), nor in our own results.

Figure 6Reconstructed atmospheric CO2 for the entire simulation period, compared to the EPICA Dome C ice-core record (Bereiter et al., 2015); compared to three different reconstructions based on ice-sheet–climate models (van de Wal et al., 2011; Stap et al., 2017; Willeit et al., 2019); and compared to chemical proxies based on 11B ratios (Hönisch et al., 2009; Seki et al., 2010; Bartoli et al., 2011; Martínez-Botí et al., 2015; Stap et al., 2016; Chalk et al., 2017; Dyez et al., 2018; Sosdian et al., 2018), alkenones (Seki et al., 2010; Badger et al., 2013; Zhang et al., 2013), and stomata (Kürschner et al., 1996; Stults et al., 2011; Bai et al., 2015; Hu et al., 2015). The pre-industrial value of 280 ppmv is shown by a dashed line.


The CO2 reconstruction for the past 800 kyr is compared to the EPICA Dome C ice-core record (Bereiter et al., 2015) in Fig. 6, as well as to several different model reconstructions and proxy-based reconstructions. The three other model-based reconstructions shown were created by decoupling the benthic δ18O signal using a 1-D ice-sheet model (van de Wal et al., 2011; Stap et al., 2017), or by using an insolation-forced, fully coupled ice-sheet–climate–carbon-cycle model (Willeit et al., 2019). The geological proxies are either based on alkenones (Seki et al., 2010; Badger et al., 2013; Zhang et al., 2013) or δ11B ratios (Hönisch et al. 2009; Seki et al., 2010; Bartoli et al., 2011; Martínez-Botí et al., 2015; Stap et al., 2016; Chalk et al., 2017; Dyez et al., 2018; Sosdian et al., 2018) derived from benthic foraminifera, or based on stomata (Kürschner et al., 1996; Stults et al., 2011; Bai et al., 2015; Hu et al., 2015).

Our results broadly match those of van de Wal et al. (2011), who reconstructed CO2 using 1-D ice models, forced with a simple spatially uniform temperature offset (assuming the strength of all feedbacks in the CO2–temperature relation to be constant in time) calculated from benthic δ18O using a similar inverse forward modelling approach. Both our results and those of van de Wal et al. (2011) show values during the late Pliocene that were higher than pre-industrial but lower than present-day values (412 ppmv at the date of writing), though the values by van de Wal et al. (2011) are at the low end of the uncertainty of our own results, as can be seen in Fig. 6. Stap et al. (2017) used a 1-D ice-sheet model set-up very similar to that by van de Wal et al. (2011) but a more elaborate energy-balance model to represent the global climate. Their results differ markedly from ours and those of van de Wal et al. (2011). The reconstruction by Willeit et al. (2019), who used a coupled ice-sheet–climate–carbon cycle model, agrees with ours and that by van de Wal et al. (2011) in terms of the glacial–interglacial amplitude both before and after the MPT. However, the downward linear trend in pCO2, which Willeit et al. (2019) prescribed manually to induce the inception of Northern Hemisphere glaciation at 2.6 Myr ago, and the resulting high values during the Early Pleistocene are not visible in the other reconstructions.

Table 1Statistical comparison of the different model-based CO2 reconstructions (this study; van de Wal et al., 2011; Stap et al., 2017; Willeit et al., 2019) to the EPICA Dome C ice-core record (Bereiter et al., 2015), as well as the glacial–interglacial amplitude in CO2 both before and after the MPT for the different reconstructions.

Download Print Version | Download XLSX

None of the model-based reconstructions agrees better than the others with the available proxy data. A recent study by Badger et al. (2019) showed that the alkenone proxy has only a very low sensitivity to atmospheric CO2 at low-to-moderate CO2 concentrations, as is clearly visible in the top panel of Fig. 6, where the alkenone proxy data show a nearly constant CO2 concentration throughout the last glacial cycle. They state that the reliability of this proxy for values lower than ∼350ppmv (the entire Pleistocene and large parts of the late Pliocene) is doubtful. In addition, several studies (Indermühle et al., 1999; Jordan, 2011; Porter et al., 2019) have questioned the reliability of the stomatal proxies. However, even when we disregard both of these proxies, the level of disagreement between the different boron-isotope-based proxies is such that using them to choose one model-based reconstruction over another is not possible. The only available data record which is reliable enough to allow such a comparison is the ice-core record (Bereiter et al., 2015). In Table 1, we compare the performance of the different model-based studies at reproducing this CO2 record. All reconstructions have been prescribed an “optimized time lag” to find the highest possible correlation with the ice-core record. Our results show the best agreement in terms of the coefficient of determination R2=0.74, the root-mean-square error RMSE = 13.5 ppmv (both taken over the entire 800 kyr ice-core period), and zero time lag required to obtain these values. To put these numbers into perspective, the same numbers are listed for a simple least-squares linear fit between the LR04 benthic δ18O stack and the EPICA Dome C record. This statistical “model” yields approximately the same results as our reconstruction in terms of correlation and RMSE, although the resulting post-MPT glacial–interglacial amplitude is about 23 ppmv too small, indicating that the increased glacial–interglacial variability is not properly captured.

Figure 7Reconstructed atmospheric CO2 during the Pleistocene. The coloured regions indicate the glacial and interglacial ranges for the Early (2.6–1.2 Myr ago; blue) and Late (1.2–0.011 Myr ago; green) Pleistocene.


Figure 7 illustrates the difference between our modelled CO2 values for the Early and Late Pleistocene. Interglacial CO2 decreases from 280–300 ppmv at the beginning of the Early Pleistocene to 250–280 ppmv just before the MPT. Glacial CO2 decreases from 220–250 to 205–225 ppmv during this period, indicating a background trend of about 20 ppmv in 1.5 Myr. After the MPT, interglacial CO2 varies between 250–320 ppmv, and glacial values decrease to 170–210 ppmv, with neither range showing a clear trend. The reduced glacial–interglacial CO2 difference of 48 ppmv (averaged over the entire period) that we find in our results agrees with the value of 43 ppmv observed in the boron isotope data by Chalk et al. (2017). The values for all four model-based reconstructions, both before and after the MPT, are shown in Table 1. After the MPT, we find a glacial–interglacial difference of 87 ppmv, which agrees very well with the value of 91 ppmv over the past 800 kyr from the EPICA Dome C ice-core record (Bereiter et al., 2015).

The discussion about the nature of the Late Pleistocene glacial cycles is still ongoing. Different studies have proposed that they represent a strongly amplified response to a weak 100 kyr forcing (Ganopolski and Calov, 2011); a response to the same 40 kyr forcing as in the Early Pleistocene, with an additional threshold mechanism that causes some insolation maxima to be “skipped”, leading to 80–120 kyr cyclicity (Huybers, 2011; Tzedakis et al., 2017); or even the result of some stochastic process with a collapse threshold, containing no true periodicity (Wunsch, 2003). Since the timing of glaciations and deglaciations in our model is forced by the LR04 benthic δ18O stack, and thus its depth–age model, no meaningful conclusions regarding this discussion can be drawn from our results.

Figure 8(a) Sea level vs. atmospheric CO2 from proxy data. Blue: sea-level data from Rohling et al. (2014; last 5 Myr) vs. boron isotope CO2 data. Red: sea-level data from Elderfield et al. (2012; last 1.5 Myr) vs. boron isotope CO2 data. Gold: sea-level data from Miller et al. (2012) and Naish et al. (2009; 3.4–2.3 Myr ago) vs. boron isotope CO2 data. Purple: sea-level data from Grant et al. (2014; last 500 kyr) vs. the EPICA Dome C ice-core record (Bereiter et al., 2015). (b–f) Reconstructed sea level and atmospheric CO2 from different modelling studies (b: this study; c: van de Wal et al., 2011; d: de Boer et al., 2014; e: Stap et al., 2017; f: Willeit et al., 2019), separated into pre- and post-MPT values. De Boer et al. (2014; d) did not explicitly model CO2, calculating only a global mean temperature change. We converted this to CO2, using the logarithmic relation from van de Wal et al. (2011).


The non-linear, time-dependent relation between CO2 and sea level is visualized in Fig. 8, both for the different proxy data (Fig. 8a) and for the different model reconstructions (Fig. b–f) using scatter plots. This visualization contains several useful features; the general shape of the curve shows the sensitivity of the world's ice sheets to changes in climate, while the spread of the data points around the curve indicates the time lag between the two. An instantaneous response would result in zero spread, while a long lag would yield a wide spread. This time lag is caused by the slow response of ice sheets to climate change; the intrinsic hysteresis of global ice-sheet volume under a changing climate; and other slow-responding physical processes, such as englacial thermodynamics, glacial-isostatic adjustment, and oceanic carbon cycling.

The combination of ice-core CO2 data (Bereiter et al., 2015) and Red Sea sea-level data (Grant et al., 2014), regarded as the two most reliable proxy reconstructions for the colder-than-present worlds of the Late Pleistocene and shown by the purple data points in Fig. 8a, displays a pattern that is very similar to the model results by de Boer et al. (2014), and by this study. The reconstruction by van de Wal et al. (2011) shows a narrower range of sea levels for each CO2 value, indicating that their ice sheets respond faster to changes in climate. The reconstruction by Willeit et al. (2019) show a wider spread in the 250–350 ppmv range, implying the existence of medium-sized ice sheets (∼30m sea-level equivalent) even when CO2 is above 280 ppmv. The reconstruction by Stap et al. (2017) shows a vast spread, indicating that their ice sheets respond very slowly. The cause of this is not known (Lennert B. Stap, personal communication, 2019). While both the observational data in Fig. 8a and the results by Willeit et al. (2019) show the existence of medium-sized ice sheets under CO2 concentrations as high as 300 ppmv, the results by van de Wal et al. (2011), de Boer et al. (2014), and this study all show a “threshold” for glaciation and sea-level drop around 250 ppmv. This indicates that the effects of orbital changes on the global climate, and by extension on the surface mass balance and the timing of glacial inceptions, are underestimated in these models.

The much weaker response of sea level to changes in CO2 for concentrations higher than 280 ppmv, which is displayed by all models except Stap et al. (2017), seems to be confirmed by the sea-level data by Miller et al. (2012) and the boron isotope CO2 data. The sea-level data by Rohling et al. (2014) show much more variation, especially during the Early Pleistocene. The results by van de Wal et al. (2011) show a very clear threshold around 375 ppmv CO2, beyond which sea level rises rapidly with increasing CO2. As indicated by the colours, they only find these warm greenhouse worlds during the first half of the Miocene (20–12 Myr ago), which explains why none of the other models, nor the proxy data, display this behaviour. Although some of the boron isotope proxies suggest CO2 concentrations higher than this 375 ppmv threshold value, the sea level data do not show the associated strong rise.

4 Conclusions and discussion

We have presented a new time-continuous, self-consistent reconstruction of atmospheric CO2, ice-sheet evolution, and global mean sea level of the last 3.6 Myr, based on the benthic δ18O glacial proxy. This reconstruction was created by using a hybrid ice-sheet–climate model (Berends et al., 2018a) to decouple the contributions from ice volume and deep-water temperature to the benthic δ18O record (Lisiecki and Raymo, 2005). Our sea-level reconstruction agrees well with similar model-based results (de Boer et al., 2014), the reconstruction based on Red Sea δ18O by Grant et al. (2014), and the combined results of benthic δ18O (Naish et al., 2009) and geological backstripping (Miller et al., 2012). Our results agree less well with the reconstruction based on Mediterranean Sea δ18O by Rohling et al. (2014), which goes back further in time but has a lower signal-to-noise ratio, and with the reconstruction based on benthic δ18O and MgCa ratios by Elderfield et al. (2012). Although previous studies (van de Wal et al., 2011; de Boer et al., 2014; Stap et al., 2017) used simpler parameterized ice-sheet and climate models, the resulting sea-level reconstructions are very similar.

Our CO2 reconstruction agrees well with the EPICA Dome C ice-core record (Bereiter et al., 2015), showing the strongest correlation of all four model-based reconstructions. Our results do not agree well with different chemical proxy-based reconstructions based on foraminiferal alkenones (Seki et al., 2010; Badger et al., 2013; Zhang et al., 2013), foraminiferal δ11B ratios (Seki et al., 2010; Martínez-Botí et al., 2015; Stap et al., 2016; Chalk et al., 2017; Dyez et al., 2018; Sosdian et al., 2018), or fossil plant stomata (Kürschner et al., 1996; Stults et al., 2011; Bai et al., 201; Hu et al., 2015). However, the strong spread between those different proxies and the large uncertainty of each individual proxy record do not allow concluding whether the difference between those proxies and our model reconstructed values is significant.

Between 2.8 and 1.2 Myr ago, during the Early Pleistocene, we simulate 40 kyr glacial cycles with glacial–interglacial sea-level changes of 25–50 m. During this period, CO2 varies between 270–280 ppmv during interglacials and 210–240 ppmv during glacial maxima, with both values decreasing by about 20 ppmv over the 1.4 Myr course of the Early Pleistocene. After the Mid-Pleistocene Transition, when the glacial cycles change from 40 kyr to 80–120 kyr cyclicity, their sea-level amplitude increases to 70–120 m, with CO2 varying between 250–320 ppmv during interglacials and 170–210 ppmv during glacial maxima. This implies a glacial–interglacial contrast of about 45 ppmv pre-MPT and 85 ppmv post-MPT. The two most reliable proxy records, the sea-level reconstruction by Grant et al. (2014) and the EPICA Dome C ice-core record (Bereiter et al., 2015), display a relation between sea level and CO2 for colder-than-present worlds that matches our model results. For warmer-than-present worlds, the large spread and uncertainty in available proxy data for both sea level and CO2, as well as the much larger uncertainty in our model results during these warm periods, prevent us from drawing conclusions on the reliability of our methodology in this regard. We believe that the inverse modelling method of using benthic δ18O, combined with coupled or hybrid ice-sheet–climate models, as a proxy for past CO2 and sea level can add valuable insights into the evolution of the Earth system during these warm episodes.

During the warm late Pliocene, we find large variability in CO2, with a difference of about 100 ppmv between the cold glacial excursion 3.3 Myr ago and the warm peak around 3.2 Myr ago. The boron-isotope-based data by Martínez-Botí et al. (2015), which have the highest temporal resolution and longest temporal range of all proxy records, show a variability of about 150 ppmv during this period (albeit with an uncertainty of about 100 ppmv in either direction). While the large discrepancies between different proxy records are too large to draw any definitive conclusions, their findings do seem to support strong CO2 variability during these warm periods. Our results for this period also show a much larger uncertainty range than during the Pleistocene. We suspect that the relative sparsity of our current matrix in this warm regime biases the model towards large ice sheets in warm worlds. The high variability in our reconstructed CO2 during the late Pliocene can thus be explained as a consequence of this bias and the requirement posed by the inverse modelling approach that the prescribed δ18O record is reproduced. If a high (warm, little ice) δ18O is prescribed but little or no ice-sheet retreat occurs in the model, then this must be compensated for by an additional increase in deep-water temperature, which requires large changes in surface climate, and therefore in CO2. By adding additional GCM snapshots for worlds with >400ppmv CO2 and smaller-than-present ice sheets (particularly a reduced East Antarctic Ice Sheet), modelled ice-sheet retreat in warm climates might be increased, reducing simulated CO2 variability during these periods. While this means that the modelled sea-level rise during this period will be larger, our current results are at the low end of the uncertainties of other reconstructions, implying that a higher modelled value will not immediately lead to a mismatch.

The parameterized relation between the global climate, deep-sea temperature, and resulting contribution to the δ18O signal could be improved upon. In this study, we used the parameterization developed by Bintanja and van de Wal (2008), which is based on the fact that the LR04 stack consists predominantly of Atlantic records and on the assumption that Atlantic deep-water temperature is strongly related to annual mean surface temperature in the North Atlantic, where deep water is formed by downwelling. This constitutes a strong simplification of the complex relation between atmospheric and oceanic temperatures, and one that is no longer necessary in our model, where the globally uniform temperature offset used by Bintanja and van de Wal (2008) has been replaced with the GCM snapshots used in our matrix method. Based on output from available ocean-enabled GCM simulations of the LGM, a more elaborate parameterization could be constructed, perhaps using surface temperatures over specific downwelling regions. Instead of using a global benthic δ18O stack, it might even possible to separate the records from different ocean basins.

Although the margins of error we report for our CO2 and global mean sea-level reconstructions are relatively small, these describe only the uncertainty arising from the uncertainty in the prescribed δ18O forcing. The true uncertainty in our reconstructions depends on many other factors, including the many parameterizations included in our model (such as the matrix method and the relation between atmospheric and oceanic temperatures and benthic δ18O), and is definitely large than the margins of error reported here. However, the uncertainties arising from these factors can unfortunately not be evaluated at this stage.

Our results should not be interpreted as a realistic reconstruction of what the world looked like in terms of global climate, ice-sheet geometry, sea level, and CO2 during these periods of geological history. Rather, we believe they should be viewed as scenarios which can help in interpreting an expected new ice-core record. For example, one hypothesized mechanism behind the MPT is regolith erosion, which led to a change in basal conditions from mostly sliding ice to mostly non-sliding ice over North America and Eurasia during the MPT (Clark and Pollard, 1998; Willeit et al., 2019). Basal conditions in our model are constant over time, so that the MPT is entirely attributed to changes in CO2. A simulation identical to ours but including a change in basal conditions would result in a different pre-MPT CO2 history. Comparing the new ice-core record to these two different scenarios could help determine if regolith erosion was a crucial factor in the transition. Other proposed physical mechanisms – such as regolith erosion leading to a different relation between ice-sheet geometry and glaciogenic dust, changes in ocean circulation leading to a different dynamic equilibrium between atmospheric and oceanic carbon, or a slow background decrease in CO2 caused by weathering or reduced volcanism – could be explored in a similar manner. We believe our modelling approach is very suited for this kind of exploratory analysis.

Appendix A

Here, we briefly summarize the equations governing the matrix method. A more thorough explanation of this approach can be found in the original publication by Berends et al. (2018a).

Each GCM snapshot contains data fields describing global monthly mean 2 m air temperature, global monthly mean precipitation, and surface elevation. The interpolation routine for the temperature field uses modelled CO2 and modelled “absorbed insolation” as weighting parameters to interpolate between snapshots. This approach captures the two most important physical processes through which a change in ice-sheet geometry affects the local and global temperature: the ice–albedo feedback and the altitude–temperature feedback. Since the absorbed insolation changes not only through changes in albedo but also through changes in incoming insolation, the effects of orbital forcing are also (indirectly) accounted for. The weighting factor wCO2 is calculated as

(A1) w CO 2 = p CO 2 - p CO 2 , LGM p CO 2 , PI - p CO 2 , LGM ,

with pCO2,PI=280ppmv and pCO2,LGM=190ppmv. Multiplying the modelled surface albedo α with the prescribed insolation at the top of the atmosphere QTOA yields the absorbed insolation Iabs:

(A2) I abs ( x , y ) = 1 - α ( x , y ) Q TOA ( x , y ) .

This value is scaled between the PI and LGM reference fields to obtain the weighting parameter wins:

(A3) w ins ( x , y ) = I abs,mod ( x , y ) - I abs,LGM ( x , y ) I abs,PI ( x , y ) - I abs,LGM ( x , y ) .

This weighting field is partly smoothed with a 200 km Gaussian filter F, to account for both local and regional effects:

(A4) w ice , T ( x , y ) = 1 7 w ins ( x , y ) + 3 7 F w ins ( x , y ) + 3 7 w ins .

This is combined with the scalar pCO2 weight wCO2 to yield the final temperature weighting parameter wT:

(A5) w T ( x , y ) = w CO 2 + w ice , T ( x , y ) 2 .

This weighting field is used to interpolate both the temperature and orography fields of the GCM snapshots:


Orography is interpolated as well, so that it can be used to downscale the temperature field from the low-resolution GCM grid to the higher-resolution ice-model grid, using a simple lapse-rate approach:

(A8) T ( x , y ) = T ref ( x , y ) + λ LGM ( x , y ) h ( x , y ) - h ref ( x , y ) .

The spatially variable lapse rate λLGM(x,y) is derived from the temperature field of the GCM snapshot.

The interpolation for precipitation is based on ice thickness rather than absorbed albedo. This reflects the fact that precipitation over the ice sheet is strongly affected by altitude:


Finally, the precipitation field is downscaled from the low-resolution GCM grid to the higher-resolution ice-model grid, using the precipitation model by Roe and Lindzen (2001) and Roe (2002) to correct for the steeper slopes in the ice model:

(A12) P ( x , y ) = P ref ( x , y ) P Roe ( x , y ) P Roe ref ( x , y ) .
Code availability

The source code of ANICE2.1 is available online at (Berends et al., 2018b). Note that the model code can be compiled but cannot be run without input data describing present-day climate and topography, initial ice thickness and topography, and GCM output files constituting the climate matrix. For any questions regarding ANICE, please contact

Data availability

The reconstructed global mean sea level and atmospheric CO2 concentration over the past 3.6 Myr are available online at (Berends et al., 2020).

Author contributions

CJB, BdB, and RSWvdW designed the study. CJB created the model set-up and carried out the simulations, with support from RSWvdW. CJB drafted the paper, and all authors contributed to the final version.

Competing interests

The authors declare that they have no conflict of interest.

Special issue statement

This article is part of the special issue “Oldest Ice: finding and interpreting climate proxies in ice older than 700 000 years (TC/CP/ESSD inter-journal SI)”. It is not associated with a conference.


The Ministry of Education, Culture and Science (OCW), in the Netherlands, provided financial support for this study via the programme of the Netherlands Earth System Science Centre (NESSC). This work was sponsored by NWO Exact and Natural Sciences for the use of supercomputer facilities. Model runs were performed on the LISA Computer Cluster; we would like to acknowledge SurfSARA Computing and Networking Services for their support. Special thanks go to Bärbel Hönisch for advice on palaeo-CO2 data usage.

Financial support

This research has been supported by the Netherlands Earth System Science Centre (grant no. 024.002.001) and the NWO Earth and Life Sciences (grant no. 863.15.019).

Review statement

This paper was edited by Ed Brook and reviewed by Andrey Ganopolski and one anonymous referee.


Alley, R. B.: The Younger Dryas cold interval as viewed from central Greenland, Quaternary Sci. Rev., 19, 213–226, 2000. 

Badger, M. P. S., Schmidt, D. N., Mackensen, A., and Pancost, R. D.: High-resolution alkenone palaeobarometry indicates relatively stable pCO2 during the Pliocene (3.3–2.8 Ma), Philos. T. R. Soc. A, 371, 20130094,, 2013. 

Badger, M. P. S., Chalk, T. B., Foster, G. L., Bown, P. R., Gibbs, S. J., Sexton, P. F., Schmidt, D. N., Pälike, H., Mackensen, A., and Pancost, R. D.: Insensitivity of alkenone carbon isotopes to atmospheric CO2 at low to moderate CO2 levels, Clim. Past, 15, 539–554,, 2019. 

Bai, Y.-J., Chen, L.-Q., Ranhotra, P. S., Wang, Q., Wang, Y.-F., and Li, C.-S.: Reconstructing atmospheric CO2 during the Plio-Pleistocene transition by fossil Typha, Glob. Change Biol., 21, 874–881, 2015. 

Bartoli, G., Hönisch, B., and Zeebe, R. E.: Atmospheric CO2 decline during the Pliocene intensification of Northern Hemisphere glaciations, Paleoceanography, 26, PA4213,, 2011. 

Beerling, D. J. and Royer, D. L.: Convergent Cenozoic CO2 history, Nat. Geosci., 4, 418–420, 2011. 

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

Berends, C. J., de Boer, B., and van de Wal, R. S. W.: Application of HadCM3@Bristolv1.0 simulations of paleoclimate as forcing for an ice-sheet model, ANICE2.1: set-up and benchmark experiments, Geosci. Model Dev., 11, 4657–4675,, 2018a. 

Berends, C., de Boer, B., and van de Wal, R.: Berends_etal_2018_GMD_code, available at: (last access: 1 February 2021), Zenodo, 2018b. 

Berends, C. J., de Boer, B., Dolan, A. M., Hill, D. J., and van de Wal, R. S. W.: Modelling ice sheet evolution and atmospheric CO2 during the Late Pliocene, Clim. Past, 15, 1603–1619,, 2019. 

Berends, C., de Boer, B., and van de Wal, R.: Berends_etal_2020_CP_Supplement, Data set, Zenodo,, 2020. 

Bintanja, R., van de Wal, R., and Oerlemans, J.: Modelled atmospheric temperatures and global sea levels over the past million years, Nature, 437, 125–128, 2005. 

Bintanja, R. and van de Wal, R. S. W.: North American ice-sheet dynamics and the onset of 100,000-year glacial cycles, Nature, 454, 869–872, 2008. 

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

Brovkin, V., Ganopolski, A., Archer, D., and Munhoven, G.: Glacial CO2 cycle as a succession of key physical and biogeochemical processes, Clim. Past, 8, 251–264,, 2012. 

Bueler, E. and Brown, J.: Shallow shelf approximation as a “sliding law” in a thermomechanically coupled ice sheet model, J. Geophys. Res., 114, F03008,, 2009. 

Chalk, T. B., Hain, M. P., Foster, G. L., Rohling, E. J., Sexton, P. F., Badger, M. P. S., Cherry, S. G., Hasenfratz, A. P., Haug, G. H., Jaccard, S. L., Martínez-García, A., Pälike, H., Pancost, R. D., and Wilson, P. A.: Causes of ice age intensification across the Mid-Pleistocene Transition, P. Natl. Acad. Sci. USA, 114, 13114–13119, 2017. 

Clark, P. U. and Pollard, D.: Origin of the middle Pleistocene transition by ice sheet erosion of regolith, Paleoceanography, 13, 1–9, 1998. 

Dansgaard, W.: Stable isotopes in precipitation, Tellus, 16, 436–468, 1964. 

de Boer, B., van de Wal, R. S. W., Bintanja, R., Lourens, L. J., and Tuenter, E.: Cenozoic global ice-volume and temperature simulations with 1-D ice-sheet models forced by benthic δ18O records, Ann. Glaciol., 51, 23–33, 2010. 

de Boer, B., van de Wal, R. S. W., Lourens, L. J., and Bintanja, R.: Transient nature of the Earth's climate and the implications for the interpretation of benthic δ18O records, Palaeogeogr. Palaeocl., 335–336, 4–11, 2012. 

de Boer, B., van de Wal, R., Lourens, L. J., Bintanja, R., and Reerink, T. J.: A continuous simulation of global ice volume over the past 1 million years with 3-D ice-sheet models, Clim. Dynam., 41, 1365–1384, 2013. 

de Boer, B., Lourens, L. J., and van de Wal, R. S. W.: Persistent 400,000-year variability of Antarctic ice volume and the carbon cycle is revealed throughout the Plio-Pleistocene, Nat. Commun., 5, 2999,, 2014. 

de Boer, B., Haywood, A. M., Dolan, A. M., Hunter, S. J., and Prescott, C. L.: The Transient Response of Ice Volume to Orbital Forcing During the Warm Late Pliocene, Geophys. Res. Lett., 44, 10486–10494, 2017. 

Dolan, A. M., Haywood, A. M., Hill, D. J., Dowsett, H. J., Hunter, S. J., Lunt, D. J., and Pickering, S. J.: Sensitivity of Pliocene ice sheets to orbital forcing, Palaeogeogr. Palaeocl., 309, 98–110, 2011. 

Dolan, A. M., Haywood, A. M., Hunter, S. J., Tindall, J. C., Dowsett, H. J., Hill, D. J., and Pickering, S. J.: Modelling the enigmatic Late Pliocene Glacial Event – Marine Isotope Stage M2, Global Planet. Change, 128, 47–60, 2015. 

Dyez, K. A., Hönisch, B., and Schmidt, G. A.: Early Pleistocene Obliquity-Scale pCO2 Variability at ∼1.5 Milion Years Ago, Paleoceanography and Paleoclimatology, 33, 1270–1291, 2018. 

Elderfield, H., Ferretti, P., Greaves, M., Crowhurst, S., McCave, I. N., Hodell, D. A., and Piotrowski, A. M.: Evolution of Ocean Temperature and Ice Volume Through the Mid-Pleistocene Climate Transition, Science, 337, 704–711, 2012. 

Finsinger, W. and Wagner-Cremer, F.: Stomatal-based inference models for reconstruction of atmospheric CO2 concentration: a method assessment using a calibration and validation approach, The Holocene, 18, 757–764, 2009. 

Foster, G. L. and Rae, J. W. B.: Reconstructing Ocean pH with Boron Isotopes in Foraminifera, Annu. Rev. Earth Pl. Sc., 44, 207–237, 2016. 

Franks, P. J., Royer, D. L., Beerling, D. J., Van de Water, P. K., Cantrill, D. J., Barbour, M. M., and Berry, J. A.: New constraints on atmospheric CO2 concentration for the Phanerozoic, Geophys. Res. Lett., 41, 4685–4694, 2014. 

Ganopolski, A. and Calov, R.: The role of orbital forcing, carbon dioxide and regolith in 100 kyr glacial cycles, Clim. Past, 7, 1415–1425,, 2011. 

Gordon, C., Cooper, C., Senior, C. A., Banks, H., Gregory, J. M., Johns, T. C., Mitchell, J. F. B., and Wood, R. A.: The simulation of SST, sea ice extents and ocean heat transports in a version of the Hadley Centre coupled model without flux adjustments, Clim. Dynam., 16, 147–168, 2000. 

Grant, K. M., Rohling, E. J., Bronk Ramsey, C., Cheng, H., Edwards, R. L., Florindo, F., Heslop, D., Marra, F., Roberts, A. P., Tamisiea, M. E., and Williams, F.: Sea-level variability over five glacial cycles, Nat. Commun., 5, 5076,, 2014. 

Haywood, A. M., Hill, D. J., Dolan, A. M., Otto-Bliesner, B. L., Bragg, F., Chan, W.-L., Chandler, M. A., Contoux, C., Dowsett, H. J., Jost, A., Kamae, Y., Lohmann, G., Lunt, D. J., Abe-Ouchi, A., Pickering, S. J., Ramstein, G., Rosenbloom, N. A., Salzmann, U., Sohl, L., Stepanek, C., Ueda, H., Yan, Q., and Zhang, Z.: Large-scale features of Pliocene climate: results from the Pliocene Model Intercomparison Project, Clim. Past, 9, 191–209,, 2013. 

Haywood, A. M. and Valdes, P. J.: Modelling Pliocene warmth: contribution of atmosphere, oceans and cryosphere, Earth Planet. Sc. Lett., 218, 363–377, 2003. 

Hönisch, B., Hemming, N. G., Archer, D., Siddall, M., and McManus, J. F.: Atmospheric carbon dioxide concentration across the mid-Pleistocene transition, Science, 324, 1551–1554, 2009. 

Hu, J.-J., Xing, Y.-W., Turkington, R., Jacques, F. M. B., Su, T., Huang, Y.-J., and Zhou, Z.-K.: A new positive relationship between pCO2 and stomatal frequency in Quercus guyavifolia (Fagaceae): a potential proxy for palaeo-CO2 levels, Ann. Bot., 115, 777–788, 2015. 

Huybers, P.: Combined obliquity and precession pacing of late Pleistocene glaciations, Nature, 480, 229–232, 2011. 

Indermühle, A., Stauffer, B., Stocker, T. F., Raynaud, D., and Barnola, J.-M.: Technical Comment: Early Holocene Atmospheric CO2 Concentrations, Science, 286, 1815 pp.,, 1999. 

IPCC: Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Solomon, S., Qin, D., Manning, M., Chen, Z., Marquis, M., Averyt, K. B., Tignor, M., and Miller, H. L., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 2007. 

Jordan, G. J.: A critical framework for the assessment of biological palaeoproxies: predicting past climate and levels of atmospheric CO2 from fossil leaves, New Phytol., 192, 29–44, 2011. 

Jouzel, J., Alley, R. B., Cuffey, K. M., Dansgaard, W., Grootes, P., Hoffmann, G., Johnsen, S. J., Koster, R. D., Peel, D., Shuman, C. A., Stievenard, M., Stuiver, M., and White, J.: Validity of the temperature reconstruction from water isotopes in ice cores, J. Geophys. Res., 102, 26471–26487, 1997. 

Jouzel, J., Masson-Delmote, V., Cattani, O., Dreyfus, G., Falourd, S., Hoffmann, G., Minster, B., Nouet, J., Barnola, J. M., Chappellaz, J., Fischer, H., Gallet, J. C., Johnsen, S., Leuenberger, M., Loulergue, L., Luethi, D., Oerter, H., Parrenin, F., Raisbeck, G., Raynaud, D., Schilt, A., Schwander, J., Selmo, E., Souchez, R., Spahni, R., Stauffer, B., Steffensen, J. P., Stenni, B., Stocker, T. F., Tison, J. L., Werner, M., and Wolff, E. W.: Orbital and Millenial Antarctic Climate Variability over the Past 800,000 Years, Science, 317, 793–797, 2007. 

Kindler, P., Guillevic, M., Baumgartner, M., Schwander, J., Landais, A., and Leuenberger, M.: Temperature reconstruction from 10 to 120 kyr b2k from the NGRIP ice core, Clim. Past, 10, 887–902,, 2014. 

Kürschner, W. M., van der Burgh, J., Visscher, H., and Dilcher, D. L.: Oak leaves as biosensors of late Neogene and early Pleistocene paleoatmospheric CO2 concentrations, Mar. Micropaleontol., 27, 299–312, 1996. 

Laskar, J., Robutel, P., Gastineau, M., Correia, A. C. M., and Levrard, B.: A long-term numerical solution for the insolation quantities of the Earth, Astron. Astrophys., 428, 261–285, 2004. 

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

Martin, M. A., Winkelmann, R., Haseloff, M., Albrecht, T., Bueler, E., Khroulev, C., and Levermann, A.: The Potsdam Parallel Ice Sheet Model (PISM-PIK) – Part 2: Dynamic equilibrium simulation of the Antarctic ice sheet, The Cryosphere, 5, 727–740,, 2011. 

Martínez-Botí, M. A., Foster, G. L., Chalk, T. B., Rohling, E. J., Sexton, P. F., Lunt, D. J., Pancost, R. D., Badger, M. P. S., and Schmidt, D. N.: Plio-Pleistocene climate sensitivity evaluated using high-resolution CO2 records, Nature, 518, 49–54, 2015. 

Miller, K. G., Mountain, G. S., Wright, J. D., and Browning, J. V.: A 180-million-year record of sea level and ice volume variations from continental margin and deap-sea isotopic records, Oceanography, 24, 40–53, 2011. 

Miller, K. G., Wright, J. D., Browning, J. V., Kulpecz, A., Kominz, M., Naish, T., Cramer, B. S., Rosenthal, Y., Peltier, W. R., and Sosdian, S.: High tide of the warm Pliocene: Implications of global sea level for Antarctic deglaciation, Geology, 40, 407–410, 2012. 

Morland, L. W.: Unconfined ice-shelf flow, in: Dynamics of the West Antarctic Ice Sheet, edited by: van der Veen, C. J. and Oerlemans, J., Kluwer Acad., Dordrecht, the Netherlands, 99–116, 1987. 

Morland, L. W. and Johnson, I. R.: Steady motion of ice sheets, J. Glaciol., 25, 229–246, 1980. 

Naish, T., Powell, R., Levy, R., Wilson, G., Scherer, R., Talarico, F., Krissek, L., Niessen, F., Pompilio, M., Wilson, T., Carter, L., DeConto, R. M., Huybers, P., McKay, R. M., Pollard, D., Ross, J., Winter, D., Barrett, P., Browne, G., Cody, R., Cowan, E., Crampton, J., Dunbar, G., Dunbar, N., Florindo, F., Gebhardt, C., Graham, I., Hannah, M., Hansaraj, D., Harwood, D., Helling, D., Henrys, S., Hinnov, L., Kuhn, G., Kyle, P., Laufer, A., Maffioli, P., Magens, D., Mandernack, K., McIntosh, W., Millan, C., Morin, R., Ohneiser, C., Paulsen, T., Persico, D., Raine, I., Reed, J., Riesselman, C. R., Sagnotti, L., Schmitt, D., Sjunneskog, C., Strong, P., Taviani, M., Vogel, S., Wilch, T., and Williams, T.: Obliquity-paced Pliocene West Antarctic ice sheet oscillations, Nature, 458, 322–328, 2009. 

Pollard, D.: A retrospective look at coupled ice sheet-climate modeiing, Climatic Change, 100, 173–194, 2010. 

Pollard, D. and DeConto, R. M.: Modelling West Antarctic ice sheet growth and collapse through the past five million years, Nature, 458, 329–332, 2009. 

Porter, A. S., Evans-FitsGerald, C., Yiotis, C., Montañez, I. P., and McElwain, J. C.: Testing the accuracy of new paleoatmospheric CO2 proxies based on plant stable isotopic composition and stomatal traits in a range of simulated paleoatmospheric O2:CO2 ratios, Geochim. Cosmochim. Ac., 259, 69–90, 2019. 

Roe, G. H.: Modeling precipitation over ice sheets: an assessment using Greenland, J. Glaciol., 48, 70–80, 2002. 

Roe, G. H. and Lindzen, R. S.: The Mutual Interaction between Continental-Scale Ice Sheets and Atmospheric Stationary Waves, J. Climate, 14, 1450–1465, 2001. 

Rohling, E. J., Foster, G. L., Marino, G., Roberts, A. P., Tamisiea, M. E., and Williams, F.: Sea-level and deep-sea-temperature variability over the past 5.3 million years, Nature, 508, 477–492, 2014. 

Seki, O., Foster, G. L., Schmidt, D. N., Mackensen, A., Kawamure, K., and Pancost, R. D.: Alkenone and boron-based Pliocene pCO2 records, Earth Planet. Sc. Lett., 292, 201–211, 2010. 

Shakun, J. D., Lea, D. W., Lisiecki, L. E., and Raymo, M. E.: An 800-kyr record of global surface ocean δ18O and implications for ice volume-temperature coupling, Earth Planet. Sc. Lett., 426, 58–68, 2015. 

Singarayer, J. S. and Valdes, P. J.: High-latitude climate sensitivity to ice-sheet forcing over the last 120 kyr, Quaternary Sci. Rev., 29, 43–55, 2010. 

Sosdian, S. and Rosenthal, Y.: Deep-sea temperature and ice volume cange across the Pliocene-Pleistocene climate transitions, Science, 325, 306–310, 2009. 

Sosdian, S. M., Greenop, R., Hain, M. P., Foster, G. L., Pearson, P. N., and Lear, C. H.: Constraining the evolution of Neogene ocean carbonate chemistry using the boron isotope pH proxy, Earth Planet. Sc. Lett., 498, 362–376, 2018. 

Stap, L. B., de Boer, B., Ziegler, M., Bintanja, R., Lourens, L. J., and van de Wal, R.: CO2 over the past 5 million years: Continuous simulation and new δ11B-based proxy data, Earth Planet. Sc. Lett., 439, 1–10, 2016. 

Stap, L. B., van de Wal, R. S. W., de Boer, B., Bintanja, R., and Lourens, L. J.: The influence of ice sheets on temperature during the past 38 million years inferred from a one-dimensional ice sheet–climate model, Clim. Past, 13, 1243–1257,, 2017. 

Stults, D. Z., Wagner-Cremer, F., and Axsmith, B. J.: Atmospheric paleo-CO2 estimates based on Taxodium distichum (Cupressaceae) fossils from the Miocene and Pliocene of Eastern North America, Palaeogeogr. Palaeocl., 309, 327–332, 2011. 

Tzedakis, P. C., Crucifix, M., Mitsui, T., and Wolff, E. W.: A simple rule to determine which insolation cycles lead to interglacials, Nature, 542, 427–444, 2017. 

Valdes, P. J., Armstrong, E., Badger, M. P. S., Bradshaw, C. D., Bragg, F., Crucifix, M., Davies-Barnard, T., Day, J. J., Farnsworth, A., Gordon, C., Hopcroft, P. O., Kennedy, A. T., Lord, N. S., Lunt, D. J., Marzocchi, A., Parry, L. M., Pope, V., Roberts, W. H. G., Stone, E. J., Tourte, G. J. L., and Williams, J. H. T.: The BRIDGE HadCM3 family of climate models: HadCM3@Bristol v1.0, Geosci. Model Dev., 10, 3715–3743,, 2017. 

van de Wal, R. S. W., de Boer, B., Lourens, L. J., Köhler, P., and Bintanja, R.: Reconstruction of a continuous high-resolution CO2 record over the past 20 million years, Clim. Past, 7, 1459–1469,, 2011. 

Wagner, F., Aaby, B., and Visscher, H.: Rapid atmospheric CO2 changes associated with the 8,200-years-B.P. cooling event, P. Natl. Acad. Sci. USA, 99, 12011–12014, 2002. 

Willeit, M., Ganopolski, A., Calov, R., Robinson, A., and Maslin, M.: The role of CO2 decline for the onset of Northern Hemisphere glaciation, Quaternary Sci. Rev., 119, 22–34, 2015. 

Willeit, M., Ganopolski, A., Calov, R., and Brovkin, V.: Mid-Pleistocene transition in glacial cycles explained by declining CO2 and regolith removal, Sci. Adv., 5, eaav7337,, 2019. 

Wunsch, C.: The spectral description of climate change including the 100 ky energy, Clim. Dynam., 20, 353–363, 2003.  

Zachos, J. C., Pagani, M., Sloan, L., Thomas, E., and Billups, K.: Trends, Rhythms, and Aberrations in Global Climate 65 Ma to Present, Science, 292, 686–694, 2001. 

Zachos, J. C., Dickens, G. R., and Zeebe, R. E.: An early Cenozoic perspective on greenhouse warming and carbon-cycle dynamics, Nature, 451, 279–283, 2008. 

Zhang, Y. G., Pagani, M., Liu, Z., Bohaty, S. M., and DeConto, R. M.: A 40-million year history of atmospheric CO2, Philos. T. R. Soc. A, 371, 20130096,, 2013. 

Zwally, H. J. and Giovinetto, M. B.: Areal distribution of the oxygen-isotope ratio in Greenland, Ann. Glaciol., 25, 208–213, 1997. 

Short summary
For the past 2.6 million years, the Earth has experienced glacial cycles, where vast ice sheets periodically grew to cover large parts of North America and Eurasia. In the earlier part of this period, this happened every 40 000 years. This value changed 1.2 million years ago to 100 000 years: the Mid-Pleistocene Transition. We investigate this interesting period using an ice-sheet model, studying the interactions between ice sheets and the global climate.