Articles | Volume 15, issue 4
Clim. Past, 15, 1603–1619, 2019
Clim. Past, 15, 1603–1619, 2019

Research article 15 Aug 2019

Research article | 15 Aug 2019

Modelling ice sheet evolution and atmospheric CO2 during the Late Pliocene

Modelling ice sheet evolution and atmospheric CO2 during the Late Pliocene
Constantijn J. Berends1, Bas de Boer2, Aisling M. Dolan3, Daniel J. Hill3, and Roderik S. W. van de Wal1 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
  • 3School of Earth and Environment, University of Leeds, Leeds, UK

Correspondence: Constantijn J. Berends (


In order to investigate the relation between ice sheets and climate in a warmer-than-present world, recent research has focussed on the Late Pliocene, 3.6 to 2.58 million years ago. It is the most recent period in Earth's history when such a warm climate state existed for a significant duration of time. Marine Isotope Stage (MIS) M2 (∼3.3 Myr ago) is a strong positive excursion in benthic oxygen records in the middle of the otherwise warm and relatively stable Late Pliocene. However, the relative contributions to the benthic δ18O signal from deep ocean cooling and growing ice sheets are still uncertain. Here, we present results from simulations of the Late Pliocene with a hybrid ice-sheet–climate model, showing a reconstruction of ice sheet geometry, sea level and atmospheric CO2. Initial experiments simulating the last four glacial cycles indicate that this model yields results which are in good agreement with proxy records in terms of global mean sea level, benthic oxygen isotope abundance, ice-core-derived surface temperature and atmospheric CO2 concentration. For the Late Pliocene, our results show an atmospheric CO2 concentration during MIS M2 of 233–249 ppmv and a drop in global mean sea level of 10 to 25 m. Uncertainties are larger during the warmer periods leading up to and following MIS M2. CO2 concentrations during the warm intervals in the Pliocene, with sea-level high stands of 8–14 m above the present day, varied between 320 and 400 ppmv, lower than indicated by some proxy records but in line with earlier model reconstructions.

1 Introduction

One of the major long-term challenges posed by anthropogenic climate change is sea-level rise due to the large-scale retreat of the Greenland and Antarctic ice sheets (e.g. Church et al., 2013). However, projecting the magnitude and especially the rate of such a retreat is limited by our understanding of the interactions between global climate and the cryosphere on centennial to multi-millennial timescales, especially in a warmer-than-present climate. In order to gain more insight into the behaviour of the Earth system in such a warmer world, numerous recent studies (Bachem et al., 2017; Bragg et al., 2012; Burke et al., 2018; de Boer et al., 2015, 2017; Dolan et al., 2011, 2015; de Schepper et al., 2014; Dowsett et al., 2016; Dwyer and Chandler, 2009; Haywood et al., 2010, 2011, 2013a, b; Hill, 2015; Lunt et al., 2009, 2010, 2012; McKay et al., 2012; Miller et al., 2012; Naish et al., 2009; Naish and Wilson, 2009; Prescott et al., 2014; Sohl et al., 2009; Swann et al., 2018; Tan et al., 2017) have focussed on the Late Pliocene, 3.6 to 2.58 million years ago, since it is the most recent period in Earth history with average global temperatures staying warmer than the present day for a significant length of time. Many of these studies, particularly those carried out as part of the Pliocene Modelling Intercomparison Project (PlioMIP; Haywood et al., 2010, 2011), focus on the mid-Pliocene Warm Period (MPWP), 3.29–2.97 Myr ago. This time slab represents a relatively stable period in Earth's climate history with warmer-than-present global temperatures, lasting longer than any of the Quaternary interglacials. Since it occurred much more recently than other warm periods, the difference in continental configuration with the present is relatively small. Modelling studies generally show a global mean annual surface temperature that was more than 3 C warmer than the present day (Bragg et al., 2012; Burke et al., 2018; Haywood et al., 2013a, b; Lunt et al., 2010, 2012). Sea surface temperatures were warmer as well, with a strongly reduced meridional gradient leading to a slight warming in the tropics and a strong warming in the polar regions (Bachem et al., 2017; Csank et al., 2011a, b; Dowsett et al., 2009, 2013, 2016; Fletcher et al., 2017; Hill, 2015; Masson-Delmotte et al., 2006; Salzmann et al., 2013). Sea-level estimates range between 10 and 30 m above present-day values (de Boer et al., 2017; Dolan et al., 2011; Dowsett et al., 2016; Dwyer and Chandler, 2009; Miller et al., 2011, 2012), caused by the almost complete deglaciation of Greenland and West Antarctica. Estimates of atmospheric CO2 concentrations during this period vary between 250 and 450 ppmv (Badger et al., 2013; Bartoli et al., 2011; Martínez-Botí et al., 2015; Seki et al., 2010; Stap et al., 2016; Zhang et al., 2013).

However, global climate during both the Pliocene in general and the MPWP in particular showed significant variability. In order to describe a more general “warm Earth” state, the MPWP is therefore treated in PlioMIP as an average of several different warm peaks that may or may not have occurred synchronously around the globe (Dowsett et al., 2016; Haywood et al., 2010, 2011). Both directly before the beginning of, and relatively shortly after, the end of the MPWP, proxy records indicate that the Earth experienced periods that were apparently colder than the present day, though neither was as cold nor as long in duration as typical Late Pleistocene glaciations.

Of particular interest is the cold excursion that occurred 3.3 Myr ago, during Marine Isotope Stage (MIS) M2, shown in Fig. 1. During the 40 000 years following the warm peak of MIS MG1 (3.315 Myr ago), benthic δ18O (LR04; Lisiecki and Raymo, 2005) increased by about 0.5 ‰ and subsequently recovered, suggesting either a global cooling, an increase in ice volume on the Northern Hemisphere and/or Southern Hemisphere or both. Sea-level records (Dwyer and Chandler, 2009; 65±25 m; Naish and Wilson, 2009; 38 m; Miller et al., 2011; 34±10 m; Miller et al., 2012; 10±10 m), as well as evidence of glacial till (Gao et al., 2012; de Schepper et al., 2014) and ice-rafted debris (de Schepper et al., 2014; Bachem et al., 2017; Smith et al., 2018) support the hypothesis of at least a partial Northern Hemisphere glaciation. De Schepper et al. (2014) and Dolan et al. (2015) provide detailed overviews of available evidence for glaciation during the Pliocene in general and MIS M2 in particular. However, because most geological fingerprints that would have been left by Pliocene ice sheets and glaciers would have been overridden or eroded by waxing and waning of the much larger Pleistocene ice sheets, evidence is limited to mostly the presence or absence of ice, providing only sparse information on geographical location and little to none on the volumes of these ice sheets.

Figure 1Benthic δ18O (LR04; Lisiecki and Raymo, 2005), 65 N summer insolation (Laskar et al., 2004) and reconstructed atmospheric pCO2 from δ18O-based model reconstructions (van de Wal et al., 2011; Stap et al., 2016) and proxy data based on alkenones (Seki et al., 2010; Badger et al., 2013; Zhang et al., 2013) and 11B ratios (Seki et al., 2010; Bartoli et al., 2011; Martínez-Botí et al., 2015; Stap et al., 2016) for the Late Pliocene. Present-day values for all variables are indicated by horizontal dashed lines; MIS M2 and KM5c are indicated by vertical dashed lines.


Dolan et al. (2015) studied MIS M2 from a climatological rather than a glaciological point of view. Using the Hadley Centre Coupled Model version 3 (HadCM3) general circulation model (GCM; Gordon et al., 2000; Valdes et al., 2017; see Sect. 2.1), they performed an ensemble of simulations of global climate during MIS M2 for different postulated and fixed ice-sheet configurations and atmospheric CO2 concentrations. By comparing the results from these different equilibrium simulations to a wide range of available climatological proxies, they attempted to constrain MIS M2 ice volume estimates through the impact such ice sheets would have on the climate. However, the available proxy records from this era have relatively large uncertainties, and where information is available, it remains difficult to use this to draw sound conclusions about Northern Hemisphere ice sheet extent. They therefore concluded that available evidence from climatological proxies was unable to constrain ice volume any further.

In this study, we adopt a different approach, combining both the glaciological and climatological viewpoints. In a recent study, Berends et al. (2018) presented and evaluated a hybrid GCM–ice-sheet model, where they proposed a matrix method of model coupling (see Sect. 2.3) to force the ANICE ice-sheet model (Bintanja and Van de Wal, 2008; de Boer et al., 2013, 2014, 2017; see Sect. 2.2) with output from the HadCM3 GCM. By using output from a simulation with HadCM3 of the Last Glacial Maximum (Singarayer and Valdes, 2010), they were able to simultaneously simulate the evolution of the ice sheets on North America, Eurasia, Greenland and Antarctica throughout the last glacial cycle and their contributions to global mean sea level and benthic δ18O. They showed that their results matched proxy-based reconstructions for ice-sheet volume, ice surface temperature, seawater δ18O, deep-water temperature and benthic δ18O. This matrix method is applied here to the Late Pliocene by using HadCM3 results from Dolan et al. (2015). The hybrid GCM–ice-sheet model presented by Berends et al. (2018) is computationally efficient enough to make large ensemble simulations feasible, opening up the opportunity to study the effects of changes in palaeotopography, pCO2 and other climatological conditions, as well as the sensitivity to ice-sheet model parameters.

However, a high-resolution, time-continuous pCO2 record needed to force the model is not available for this period. We resolve this by using the inverse modelling approach that was also used by Bintanja and van de Wal (2008), de Boer et al. (2013, 2014) and Stap et al. (2016). In this approach, we compare modelled benthic δ18O to the LR04 stack (Lisiecki and Raymo, 2005) and calculate pCO2 based on the difference between the two (see Sect. 2.4). This makes our model set-up conceptually very similar to the approach by Stap et al. (2016), who also used the LR04 stack of δ18O to force a coupled ice-sheet–climate model and thus produce a pCO2 reconstruction. However, they used a relatively simple zonally averaged energy-balance climate model coupled to a 1-D ice model, whereas we use GCM output to drive 3-D ice-sheet models, making our approach more detailed in terms of the behaviour of global climate, the ice sheets and the interactions between the two, at the expense of computational requirements. The CO2 and climate reconstructions by van de Wal et al. (2011), Stap et al. (2016) and the one presented here can be viewed as proxy-based reconstructions, based on the concept that benthic δ18O is a proxy for changes in ocean temperature and land ice volume. All three studies use a climate model describing the known relations between pCO2 and temperature and ice volume, in order to determine how pCO2 must have evolved in the past in order to produce the observed benthic δ18O signal.

2 Methodology

2.1 Climate model

HadCM3 is a coupled atmosphere–ocean general circulation model (Gordon et al., 2000; Valdes et al., 2017). It accurately reproduces the heat budget of the present-day climate (Gordon et al., 2000) and has been used for future climate projections in the Intergovernmental Panel on Climate Change Fourth Assessment Report (IPCC AR4) (e.g. Solomon et al., 2007) and paleoclimate reconstructions such as PMIP2 (Braconnot et al., 2007) and PlioMIP (Haywood and Valdes, 2003; Dolan et al., 2011, 2015; Haywood et al., 2013a). The atmosphere module of HadCM3 has 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. In the model set-up by Berends et al. (2018), the climate matrix consists of two GCM snapshots of, respectively, the pre-industrial period (PI) and the Last Glacial Maximum (LGM), produced by Singarayer and Valdes (2010) with HadCM3. Here, we include several additional snapshots focussing specifically on the Pliocene.

2.2 Ice-sheet model

To simulate the evolution of the ice sheets we use ANICE, a coupled 3-D ice-sheet–shelf model (Bintanja and Van de Wal, 2008; de Boer et al., 2013, 2014, 2017). It combines the shallow shelf approximation (SSA; Morland, 1987) for floating ice shelves with the shallow ice approximation (SIA; Morland and Johnson, 1980) for grounded ice to solve the ice flow. A Mohr–Coulomb plastic law for basal sliding is included, with basal stresses included in the SSA equations. The basal stress is calculated as a function of a till stress, which in turn depends on the local bedrock elevation (Winkelmann et al., 2011; de Boer et al., 2013). For grounded ice, the velocities resulting from both approximations are summed, resulting in a smooth transition zone between slow-flowing land ice and fast-flowing floating ice. This approach allows the grounding line to respond to changes in shelf buttressing, resulting in proper glacial–interglacial differences in Antarctic ice volume (de Boer et al., 2013; Berends et al., 2018) by the advance of grounding lines in the Filcher–Ronne and Ross basins toward the continental shelf. The surface mass balance is parameterised using an insolation–temperature scheme using monthly temperatures and precipitation, refreezing of water and a correction for orographic forcing of precipitation; a more detailed model description is provided by de Boer et al. (2013) and references therein. The horizontal resolution of ANICE for this application is 20 km for Greenland and 40 km for the other three regions (North America, Eurasia and Antarctica). The highly parameterised climate forcing and resulting computational efficiency of ANICE allow for transient simulations of multiple glacial cycles to be carried out within 10–100 h on single-core systems, making ensemble simulations feasible. Melt underneath the ice shelves is calculated using a linear relation to ocean temperature change (Pollard and DeConto, 2009; Martin et al., 2011), a parameterisation of subshelf cavity circulation based on the shortest linear distance to the open ocean (Pollard and DeConto, 2009) and the glacial–interglacial variance parameterisation by Pollard and DeConto (2009). A more detailed explanation is provided by de Boer et al. (2013), who tuned this approach to produce realistic present-day Antarctic shelves and grounding lines. A simple threshold thickness of 200 m is used to describe ice calving, whereby any shelf ice below this thickness is removed.

2.3 Matrix method

Using the definition by Pollard (2010), a climate matrix is a collection of pre-calculated output data from several steady-state GCM simulations, called “snapshots”, that differ from each other in one or more key parameters, such as prescribed atmospheric pCO2, orbital configuration or ice-sheet configuration, each creating a separate dimension of the matrix. When performing a simulation with an ice-sheet model, at every point in time during the simulation the prescribed climate forcing is determined by combining the climate states constituting the matrix according to the position of the model state within the matrix. This constitutes a middle ground between methods of offline forcing, such as a glacial index method, and fully coupled ice-sheet–climate models. When the ice-sheet model is in a state corresponding to one of the GCM snapshots, the climate from this snapshot will be prescribed, containing the effects of the altitude–temperature and albedo–temperature feedbacks of the ice sheets, the effect of ice sheet geometry on large-scale atmospheric circulation and precipitation, and possibly also the effects of changed freshwater fluxed on ocean circulation, depending on the GCM. When the ice-sheet model is in a state lying in between different GCM snapshots, the prescribed climate is a spatially variable linear interpolation of these snapshots. This means that non-linearities in different feedback processes, such as the effects of ice-sheet geometry on atmospheric circulation and precipitation, are difficult to properly account for.

In this study, we use the model set-up developed by Berends et al. (2018), who created a matrix with the HadCM3 climate states of the pre-industrial and the Last Glacial Maximum from Singarayer and Valdes (2010) and used it to force the ANICE ice-sheet model. In this set-up, temperature fields from the two climate states are combined based on a prescribed value for pCO2 and on the internally modelled ice sheets, with the feedback of the ice sheets on the climate based 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 parameterised 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 resulting plateau desert caused by the presence of a large ice sheet. Berends et al. (2018) 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, showing 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.

In this study, we extended the PI-LGM climate matrix by adding several climate states from the study by Dolan et al. (2015). The four different ice sheet configurations they used are shown in Fig. 2. The “PRISM” ice sheets are based on the PRISM3 reconstruction (Dowsett et al., 2010), which is a time-slab representation of average peak warm conditions during the MPWP. The “Small” ice sheets are present-day conditions. The “Medium” and “Large” ice sheets were based on the ICE-5G reconstruction of the last deglaciation (Peltier, 2004) at 8 and 11 kyr ago, respectively. Each of these configurations was used as boundary conditions for two simulations with HadCM3: one with 280 ppmv and one with 220 ppmv pCO2, both with 3.3 Myr orbital parameters. This adds up to eight different snapshots, plus one additional “Plio_Control” simulation with the PRISM3 ice sheet configuration and orbital parameters, and 405 ppmv pCO2. These simulations allow the climate matrix to separate effects on climate by pCO2 and ice-sheet extent and provide valuable information on climates that are both warmer and colder than the present day. Although the new climate matrix is relatively sparse for warmer-than-present worlds, containing only one snapshot (Plio_Control) for pCO2> 280 ppmv and only three snapshots for smaller-than-present ice sheets (Plio_Control, PRISM_280 and PRISM_220), we believe it is still suitable for simulating the warm Pliocene. The matrix used by Berends et al. (2018) to simulate the last glacial cycle only contained two GCM snapshots in total and still produced satisfactory results.

Figure 2The four ice-sheet configurations used by Dolan et al. (2015) as boundary conditions for their HadCM3 simulations. “PRISM” is the PRISM3 ice sheet from Dowsett et al. (2010), “Small” is the present day, “Medium” is the ICE-5G reconstruction (Peltier, 2004) at 8 kyr ago, and “Large” is ICE-5G at 11 kyr ago.

This extended matrix therefore allows for a more accurate simulation of both the warm Late Pliocene and the cold MIS M2 glaciation. For the North American and Eurasian modules of ANICE, we added the simulations of the Medium and Large ice sheets (both the 280 and 220 ppmv pCO2 versions) by Dolan et al. (2015), since those provide extra information on the effect on climate of intermediate-sized ice sheets, as well as the Plio_Control simulation for its information on pCO2 levels above 280 ppmv. In the case of North America and Eurasia, we did not use the Small and PRISM states, as there is no ice on either continent, meaning these simulations contain no additional constraints for these ice-sheet models. For the Greenland and Antarctica models, we chose to only add the simulations of the two PRISM states and the Plio_Control simulation, because they provide new information on the effect on climate of smaller-than-present-day ice sheets. The Medium and Large simulations were left out of the matrix because the ICE-5G ice sheets (Peltier, 2004) that were used to force those HadCM3 simulations have the exact same horizontal extent as the ICE-5G LGM ice sheets. Not only does this make it difficult to distinguish between these states in the interpolation routines, it also means the effect on local climate, other than through the altitude–temperature feedback, is likely to have been small.

2.4 Inverse method

The inverse forward modelling approach used to determine pCO2 based on the difference between modelled and observed benthic δ18O is very similar to that described by de Boer et al. (2013). Their method calculates how the climate at high latitudes, described by a single, spatially uniform temperature offset ΔTNH, should have evolved, such that its effect on deep ocean temperature and land ice volume reproduces the observed benthic δ18O signal. This is achieved by comparing the modelled benthic δ18O value δ18Omod at every time step to the observed value δ18Oobs. If it is too positive, then either the ocean is not cold enough or there is not enough land ice. Global mean surface temperatures are then lowered in the next time step, leading to both a cooling in the deep ocean and an increase in ice growth. This relation is quantified by de Boer et al. (2013) by the following equation:

(1) Δ T NH = Δ T NH + 20 δ 18 O mod - δ 18 O obs t + 0.1 kyr .

Here, ΔTNH is the mean surface temperature anomaly between 40 and 80 latitude at sea level over the preceding 2 kyr. The modelled benthic δ18O is calculated using ice volume, ice-sheet δ18O and deep-water temperatures relative to the present day (PD) for every 100 years. 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 global mean surface temperature anomaly. The optimum values of 2 kyr for the length of the CO2 averaging window, 3 kyr for the deep-water temperature averaging window and 20 for the scaling parameter, were determined by de Boer et al. (2013), producing a value of ΔTNH=-15 K at LGM.

Since the climate matrix used in our model determines the regional climate based on the modelled ice sheet and the scalar atmospheric CO2 concentration, we adapted this approach for our model set-up by using the difference between modelled and observed δ18O to calculate a value for pCO2, which is subsequently forwarded to the climate matrix. That way, our model reconstructs how pCO2 should have evolved in order to change global climate in such a way that the resulting changes in deep ocean temperature and land ice volume reproduce the observed benthic δ18O signal. The algorithm then becomes

(2) p CO 2 = p CO 2 + 120 δ 18 O mod - δ 18 O obs t + 0.1 kyr .

As the constrained quantity is the change in pCO2, the scaling factor changes to 120 ppmv/‰ δ18O change in order to produce a glacial–interglacial contrast of 90 ppmv pCO2. Based on the results of preliminary experiments, the length of the CO2 averaging time window was increased to 8.5 kyr, in line with the higher values given by Stap et al. (2016). For forcing the inverse routine, the LR04 stack was used. Although a few different globally distributed stacks are available (e.g. Imbrie et al., 1984; Lisiecki and Raymo, 2005; Zachos et al., 2001, 2008; Cramer et al., 2009), the differences among them are small for the Late Pliocene and the Pleistocene. In order to maintain consistency with earlier reconstructions based on inverse modelling methods (de Boer et al., 2013; Stap et al., 2016; van de Wal et al., 2011), we decided to use the LR04 stack.

A conceptual visualisation of the inverse-method-forced matrix model is shown in Fig. 3.

Figure 3A conceptual visualisation of the inverse forward modelling approach. The model is forced externally by a benthic δ18O record and an insolation reconstruction (black boxes). The inverse routine calculates pCO2 based on the difference between observed and modelled δ18O. This value is forwarded to the climate matrix, which interpolates between the GCM snapshots based on the prescribed pCO2 value and the modelled state of the cryosphere (ice thickness and albedo).


2.5 Palaeotopography reconstruction

Several recent studies have investigated the evolution of bedrock topography in the geological past. Although ANICE does include a regional solid Earth model to calculate vertical bedrock movement in response to changes in ice distribution, other processes such as erosion and plate tectonics are currently not accounted for within the model and require external forcing. For the last glacial cycle benchmark simulation, such effects are assumed to be negligibly small, but this assumption might no longer be valid when going millions of years back in time.

In our study, we use the palaeotopography reconstruction of the Barents Sea area by Butt et al. (2002), shown in Fig. 4, based on a reversal of the erosion of sediments by the Pleistocene ice sheets. In order to investigate the effect on atmospheric and oceanic circulation, Hill (2015) used HadCM3 to perform simulations both with this topography reconstruction and with present-day topography, both at 405 and 220 ppmv pCO2. By subtracting the calculated climate fields (temperature and precipitation) for the palaeotopography simulation from the present-day topography simulation, and adding the resulting “fingerprint” to the climate fields generated by our climate matrix, we take into account the effect of this change in topography (and the accompanying change in the land–ocean mask) on the global climate.

Figure 4Topography change relative to the present day for the reconstruction created by Butt et al. (2002) and used by Hill (2015).

Preliminary experiments showed that forcing the ice-sheet model with the palaeotopography reconstruction without applying this climate “fingerprint” resulted in the persistent presence of a small (∼5 m sea-level-equivalent) ice sheet over the newly exposed Barents Land. The climate fingerprint obtained from the simulations by Hill (2015) changes the local climate from an oceanic to a continental climate, with colder, dryer winters and warmer, wetter summers, resulting in more summer melt and an overall more negative mass balance, implying less ice.

Although other areas of the world might have been eroded by ice sheets (i.e. the Canadian Archipelago, Dowsett et al., 2016; Antarctica, Wilson and Luyendyk, 2009), no GCM simulations investigating the effect on global climate of reversing those changes are currently available. We have therefore chosen not to apply any of these other topography reconstructions to our model.

3 Results

3.1 Last glacial cycle benchmark

In order to assess the performance of the model when calculating pCO2 with the inverse routine instead of prescribing it directly from an ice-core record, we first performed a simulation of the last four glacial cycles, similar to the work by Berends et al. (2018). The model was calibrated by tuning the ablation parameter for the four individual ice sheets such that their volumes at LGM match the ICE-5G reconstruction. We then performed a sensitivity analysis similar to the experiment described by Berends et al. (2018), investigating the sensitivity of the modelled sea-level drop and benthic δ18O to the uncertainty in the prescribed forcing (range based on the uncertainty reported by Lisiecki and Raymo, 2005), the ablation tuning parameter (range based on the allowed values found by Berends et al., 2018) and SIA/SSA enhancement factors (increasing the ice velocities calculated for isotropic ice to more closely match those calculated for anisotropic ice according to the approach by Ma et al. 2010; ranges based on the allowed values reported by Ma et al., 2010), as well as several new model parameters involved in the inverse forcing method: the averaging time for the modelled pCO2 (range based on the values reported by de Boer et al., 2013 and Stap et al., 2016), the ratio between surface temperature anomaly and deep-sea water temperature anomaly and the tuning parameter relating pCO2 to the difference between observed and modelled δ18O (ranges based on the values reported by de Boer et al., 2010, and Bintanja and van de Wal, 2008), resulting in 17 individual simulations. The values that were used for all these parameters are listed in Table 1. The 17 ensemble members thus yield an estimate of the uncertainty related to both model parameters and forcing.

Table 1Values of the different model parameters used in the last glacial cycle (LGC) sensitivity analysis using the inverse-method-forced matrix model. All model parameters were given upper and lower bounds 10 % above and below their benchmark value, except for the SIA/SSA flow enhancement factors (values based on Ma et al., 2010), the δ18O forcing record (0.1 ‰ uncertainty stated by Lisiecki and Raymo, 2005) and the pCO2 averaging time (values of 2000 and 15 000 years given by de Boer et al., 2014 and Stap et al., 2016, respectively).

Download Print Version | Download XLSX

The simulated pCO2 record is compared to the European Project for Ice Coring in Antarctica (EPICA) Dome C ice-core record (Lüthi et al., 2008) in Fig. 5. The ranges of modelled values for pCO2 and sea-level drop at LGM for all investigated model parameters are listed in Table 2. Based on these uncertainties, the model shows that LGM pCO2 is 188–197 ppmv and that the sea-level-equivalent volume of the four continental ice sheets at LGM was 83–100 m, agreeing well with the values of 185 ppmv pCO2 from the EPICA ice core and 100 m sea-level-equivalent ice volume from the ICE-5G reconstruction (Peltier, 2004). The modelled pCO2 values match the EPICA record better than the values simulated by Stap et al. (2016), as demonstrated by the linear correlation and root mean squared error (RMSE) between the EPICA Dome C record and the reconstructions; R2=0.46 and RMSE = 23.7 ppmv for Stap et al. (2016) and R2=0.71 and RMSE = 15.2 ppmv for our simulation. The reconstruction by van de Wal et al. (2011) performs very similarly to ours (R2=0.72, RMSE = 14.7 ppmv), but since it was partly derived from the EPICA record, the comparison is not independent and therefore cannot be compared to our results.

Figure 5pCO2 throughout the last four glacial cycles (410 kyr ago – PD): observations from the EPICA Dome C ice core (Lüthi et al., 2008), reconstruction by Stap et al. (2016) and results from the inverse-method-forced matrix model (this study). Solid green line shows the benchmark run; green shaded area shows the maximum uncertainty range from the sensitivity experiment with 17 ensemble members, and the dotted line indicates the pre-industrial CO2 concentration. Linear correlation coefficients R2 are shown for the correlation between modelled pCO2 and the EPICA Dome C record.


Table 2Sensitivity of the modelled pCO2 and total eustatic sea-level contribution to the different model parameters at different points in the simulations.

Download Print Version | Download XLSX

Benthic oxygen isotope abundance and its contributions from ice volume and deep-sea water temperature are shown in Fig. 6 and compared to reconstructions by Lisiecki and Raymo (2005) and by Shakun et al. (2015), who made a proxy-based decomposition of the respective contributions to the benthic δ18O from land ice and deep-sea temperature. The simulated benthic δ18O shows a near-perfect match with the LR04 stack (Lisiecki and Raymo, 2005) that was used to force the model, as is to be expected when using the inverse forward modelling approach. The observed rapid drop in δ18Osw at the inception, between 120 and 110 kyr BP, is reproduced well, as is the drop in Tdw. Deep-water temperature between 60 ka and LGM appears to be too high, more so than for the CO2-forced model version by Berends et al. (2018).

Figure 6Benthic δ18O during the LGC for the 410 kyr simulations using the inverse-method-forced matrix model, compared to data from LR04 (a; Lisiecki and Raymo, 2005) and from Shakun et al. (2015) for both the contribution from the δ18O of seawater (b) and that from deep-water temperature (c), as well as the resulting deep-water temperature itself (d). Also shown are the results from Berends et al. (2018). Solid green line shows the benchmark run; green shaded area shows the maximum uncertainty range from the sensitivity experiment. Since the new model set-up is forced with the LR04 stack, rather than with the EPICA CO2 record used by Berends et al. (2018), the increased correlation coefficients are not a strong result.


Surface temperature anomalies over Greenland and Antarctica compared to ice-core records (EPICA Dome C; Jouzel et al., 2007; GISP2; Alley, 2000; NGRIP; Kindler et al., 2014) throughout the last glacial cycle are shown in Fig. 7. The performance of the new model version in terms of ice surface temperature is comparable to that of the model by Berends et al. (2018), as illustrated by the linear correlation coefficients and root mean square error between the modelled temperatures and the ice-core records: R2= 0.87 and RMSE = 0.86 K for Antarctica in this model versus R2= 0.84 and RMSE = 0.91 K in the old version. For Greenland, the new model produces a value of R2=0.74 and RMSE = 2.2 K versus R2=0.65 and RMSE = 2.6 K for the old model.

Figure 7Modelled versus reconstructed surface temperature anomaly ΔTs for Antarctica (EPICA Dome C; Jouzel et al., 2007) and Greenland (GISP2; Alley, 2000; NGRIP; Kindler et al., 2014) for the LGC simulation using the inverse-method-forced matrix model, compared to the direct pCO2-forced matrix model by Berends et al. (2018). Solid green line shows the benchmark run; green shaded area shows the maximum uncertainty range from the sensitivity experiment. Ice-core temperature records have been subjected to a 4 kyr running average; variance is shown by black shaded area. Linear correlation coefficients (R2) are shown for the correlation between modelled ice surface temperatures and ice-core records.


The mismatch during the inception of the glacial cycle between isotope-derived Antarctic surface temperature and ice-core CO2 on the one hand and benthic δ18O and sea level on the other hand, reported by Bintanja and van de Wal (2008), van de Wal et al. (2011), de Boer et al. (2014), Niu et al. (2017) and Berends et al. (2018), is much better in the simulations here. The linear correlation coefficient R2 between modelled and reconstructed Antarctic surface temperatures between 120 and 80 kyr ago increased from a value of 0.49 for Berends et al. (2018) to a value of 0.74 for this study. For Greenland, this value increased from 0 to 0.36. The CO2-forced model (Berends et al., 2018) produced Antarctic surface temperatures that were in good agreement with the isotope-based proxy record but failed to reproduce the strong sea-level drop. The δ18O-forced model from this study reproduces benthic δ18O and its different contributions and shows a too-strong decrease in pCO2 but is in overall agreement with proxy records of CO2, sea level and temperature, indicating that there is an added value of using the climate matrix method as applied here.

3.2 Transient simulation of the Pliocene

The comparisons between model results and (proxy) data for the simulations of the last four glacial cycles indicate that the model accurately reproduces pCO2, ice volume and general geometry (not shown), and surface temperatures. We therefore proceeded to apply the new model set-up to the Late Pliocene. We chose to start our transient simulations 3.65 Myr ago, capturing the warm period between 3.6 and 3.4 Ma. The simulations were run until 2.75 Myr ago, since the density of available pCO2 proxy data is much higher after MIS M2, allowing for a more detailed comparison of modelled pCO2 to proxy-based pCO2 reconstructions. The model was initialised with the same PRISM3 ice sheets (Dowsett et al., 2010) that were also used to force the PRISM and Plio_Control HadCM3 experiments by Dolan et al. (2015). Due to the nature of the inverse coupling method, initialising the model with present-day ice sheets quickly converges to the same result. Topography was set to the present day plus the Barents Sea erosion reversal from Butt et al. (2002) and a glacial isostatic adjustment (GIA) correction accounting for the difference in ice loading over Greenland and Antarctica according to the PRISM3 reconstruction (Dowsett et al., 2010). Insolation and benthic δ18O were prescribed according to Laskar et al. (2004) and Lisiecki and Raymo (2005), respectively. In order to estimate the uncertainty in the modelled ice volume, we performed the same sensitivity analysis as for the last glacial cycle, with the same parameter values shown in Table 1. The resulting simulated pCO2 record is shown in Fig. 8 and compared to other model reconstructions (van de Wal et al., 2011; Stap et al., 2016) and to proxy-based data derived from alkenones (Seki et al., 2010; Badger et al., 2013; Zhang et al., 2013) and δ11B ratios (Seki et al., 2010; Bartoli et al., 2011; Martínez-Botí et al., 2015; Stap et al., 2016).

Figure 8pCO2 throughout the Late Pliocene and early Pleistocene as simulated with the inverse-method-forced matrix model, compared to δ18O-based model reconstructions (van de Wal et al., 2011; Stap et al., 2016) and proxy data based on alkenones (Seki et al., 2010; Badger et al., 2013; Zhang et al., 2013) and 11B ratios (Seki et al., 2010; Bartoli et al., 2011; Martínez-Botí et al., 2015; Stap et al., 2016). Solid line shows the benchmark run; shaded area shows the maximum uncertainty range from the sensitivity experiment.


The ranges of modelled values pCO2 and sea-level change at MIS M2 and at KM5c (3.205 Myr ago) for all investigated model parameters are listed in Table 2. KM5c is used because it has been identified as a time slice representing the mid-Pliocene Warm Period (Haywood et al., 2013b). These simulations show that during MIS M2, pCO2 is 233–249 ppmv and that the sea-level-equivalent volume of the four continental ice sheets was 10–25 m bigger than the present day, with the uncertainty based on the spread in the results from the ensemble of simulations. The uncertainty in modelled pCO2 becomes much larger for warmer-than-present climates, as shown by the modelled ranges for KM5c. The sensitivity to the benthic δ18O forcing is especially high, resulting in modelled pCO2 values of 303–384 ppmv. The reason for this is that the climatological forcing resulting from the climate matrix is less constrained for warmer-than-present-day climates. Whereas the matrix contains six snapshots describing climates with more ice and/or lower pCO2, there is only one ice sheet configuration smaller than the present day (PRISM) and only one snapshot with a pCO2 higher than 280 ppmv (the Plio_Control simulation, with 405 ppmv).

The resulting modelled sea-level contributions over time are shown in Fig. 9. The modelled ice sheets over the Northern Hemisphere and Southern Hemisphere at MIS M2 and KM5c are shown in Figs. 10 and 11, respectively. In North America, MIS M2 is clearly visible as a strong peak in ice volume, which immediately disappears when pCO2 rises again. Most of the ice forms over northeastern Canada, with a smaller ice sheet developing over the northern Cordillera. In Eurasia, only small ice caps form on Svalbard and Nova Zembla (no longer islands but now small mountain areas bordering the newly exposed Barents Land), with no sizeable ice sheets forming even at the peak of MIS M2. Greenland is mostly ice-free until MIS M2, when it rapidly develops an ice sheet slightly larger than the present day. After MIS M2, the ice sheet disappears, advancing and retreating several times during the following period. Similar behaviour is observed on West Antarctica, while East Antarctica remains stable throughout the simulation. While both Greenland and Antarctica continue to show substantial variability throughout the remainder of the simulations, North America does not glaciate again until the onset of the Pleistocene glacial cycles, 2.8 Myr ago.

Figure 9Volumes of the four ice sheets over time throughout the Late Pliocene and early Pleistocene as simulated with the inverse-method-forced matrix model. Solid green line shows the benchmark run; green shaded area shows the maximum uncertainty range from the sensitivity experiment. Vertical dashed lines indicate MIS M2 (3.295 Myr ago) and KM5c (3.205 Myr ago).


Figure 10The ice sheets at the peak of MIS M2 (3.295 Myr ago), as simulated with the inverse-method-forced matrix model. Contour lines for the Northern Hemisphere (a) show ice thickness; contour lines for Antarctica (b) show surface elevation. Antarctic ice shelves are shown as light blue. Bedrock elevation where not covered by ice is shown by colours. A sizeable ice sheet exists over the present-day Hudson Bay and Baffin Island, as well as a smaller one over the northern Cordillera. Antarctic ice volume increases by 1.5–3.5 m sea-level equivalent (s.l.e.) because of the grounding of ice into the Filchner–Ronne basin.

Figure 11The ice sheets during KM5c (3.205 Myr ago), as simulated with the inverse-method-forced matrix model. Contour lines for the Northern Hemisphere (a) show ice thickness; contour lines for Antarctica (b) show surface elevation. Antarctic ice shelves are shown as light blue. Bedrock elevation where not covered by ice is shown by colours. Whereas most of the ice on Greenland has disappeared, retreat on Antarctica is limited to the Ross Sea, where the present-day ice shelf disintegrates to leave open ocean.

Global mean sea level is compared to two different reconstructions in Fig. 12. Our model results generally lie between the δ18O-based reconstruction by Miller et al. (2011) and the reconstruction based on geological backstripping from New Zealand by Miller et al. (2012). During warm periods, our model generally shows lower sea levels and less variability than Miller et al. (2011, 2012). During cold periods, our model generally shows less sea-level drop than the δ18O-based reconstruction by Miller et al. (2011) but more than the reconstruction based on geological backstripping. The δ18O-based reconstruction by Miller et al. (2011) is based on a linear relation between benthic δ18O and sea level, which is an oversimplification of separating the contributions to the benthic δ18O signal (e.g. Bintanja et al., 2005; de Boer et al., 2013). Miller et al. (2012) noted that reconstructing absolute values for local relative sea level based on geological backstripping is difficult due to the required corrections for GIA and dynamic topography. However, the relatively short duration of MIS M2 means that the reconstructed drop in sea level, of about 20 m relative to the background level, is likely to be accurate. Our model produces a value of about 24 m, in good agreement with this value.

Figure 12(a) Global mean sea level relative to the present day over time throughout the Late Pliocene as simulated with the inverse-method-forced matrix model, compared to reconstructions based on δ18O (Miller et al., 2011; blue line) and geological backstripping (Miller et al., 2012; black line with triangles). Shaded areas show maximum uncertainty ranges. Vertical dashed lines indicate MIS M2 and KM5c. (b) Peak sea-level drop during MIS M2 (3.3 Myr ago) for all three reconstructions, same vertical scale. Shaded areas show uncertainty ranges.


The evolution of the West Antarctic ice sheet agrees partially with information derived from the AND-1B sediment core, recovered from beneath the northwest part of the Ross Ice Shelf by the ANDRILL programme (Naish et al., 2009; McKay et al., 2012). Information derived from this core by de Schepper et al. (2014) is compared to model results in Fig. 13. AND-1B shows ice-free conditions in the Ross Sea up to 3.4 Myr ago, followed by glacial deposits up to 3.24 Myr ago. Our model results show ice-free conditions up to 3.32 Myr ago, just prior to MIS M2. The ice-free conditions shown in our model results around KM5c cannot be validated by AND-1B due to a lack of data. Between 3.14 and 3.04 Myr ago, AND-1B again contains glacial deposits when our model results indicate ice-free conditions. The glacial conditions between 3.04 and 2.95 Myr ago and the subsequent ice-free conditions between 2.95 and 2.90 Myr ago indicated by AND-1B match with our model results. However, since the AND-1B sediment core contains several sizeable data gaps due to geological unconformities, the possibility that observed ∼40 kyr cycles in ice-rafted debris concentration have been incorrectly matched with 40 kyr cycles in the δ18O age model cannot be precluded. We therefore conclude that the AND-1B sediment core record cannot be used to confirm or refute our model results.

Figure 13Comparison of model results to the AND-1B sediment core (de Schepper et al., 2014). Panel (a) shows the glacial conditions derived from the sediment core, classified as either “open marine” (blue), “glacial deposits” (red) or “no data” (grey). Panel (b) shows the fraction of the Ross Sea that is covered by floating (red) and grounded (green) ice in our model simulations, with shaded areas showing the maximum uncertainty range from the sensitivity experiment.


4 Discussion and conclusions

We have presented a new time-continuous, self-consistent reconstruction of pCO2, ice sheet configuration and climate for the Late Pliocene, 3.65–2.75 Myr ago. Our approach is based on the matrix method by Berends et al. (2018), where an ice-sheet model is forced with a combination of several pre-calculated GCM snapshots. We have extended their two-state climate matrix with several GCM snapshots created by Dolan et al. (2015), who simulated global climate during MIS M2 for different ice-sheet configurations and pCO2 levels. Since our initial experiment, where this model was forced with the pCO2 reconstruction by Stap et al. (2016), proved unable to constrain sea level during MIS M2 any further, we adopted the inverse forward modelling approach by de Boer et al. (2013), forcing the model with the LR04 benthic δ18O stack (Lisiecki and Raymo, 2005). By first using this δ18O-forced model set-up to simulate the last glacial cycle, we showed that it performed at least equally well to the CO2-forced set-up by Berends et al. (2018) in terms of benthic δ18O (Fig. 6) and surface temperature (Fig. 7), and better than the 1-D model set-up by Stap et al. (2016) in terms of simulated pCO2 (Fig. 5).

Our results for the Late Pliocene show a global mean sea-level drop of 10–25 m during MIS M2, with the uncertainty resulting from a sensitivity analysis investigating several key model parameters and the uncertainty in the applied δ18O forcing. This value is in good agreement with the reconstruction based on geological backstripping from New Zealand by Miller et al. (2012; 10±10 m) and the δ18O-based reconstruction by Miller et al. (2011; 34±10 m). The extra ice with respect to the present day is located mostly in eastern Canada and the northern Cordillera (9–20 m s.l.e.) and the grounded ice over the Filcher–Ronne Sea (1.5–3.5 m s.l.e.). The atmospheric CO2 concentration necessary to produce the cooling required to grow these ice sheets is shown to be 233–249 ppmv. During MIS KM5c, most of the ice on Greenland and West Antarctica disappears, raising global mean sea level to 3–10 m above the present day, caused by a pCO2 of 303–384 ppmv. The sea-level high stand of the mid-Pliocene Warm Period is achieved during MIS KM3 (3.155 Myr ago) at 8–14 m above the present day. The larger uncertainty in the modelled pCO2 during warmer periods is attributed to the fact that the climate matrix used to force our ice sheet model contains only one GCM snapshot with a pCO2 above present-day levels and only one ice sheet configuration with smaller-than-present ice sheets. Hence, the relationship between ice sheets and climate for warmer-than-present worlds is poorly constrained, which is reflected by an increased uncertainty in the simulated pCO2 and ice volume. Expanding the climate matrix with additional GCM snapshots for intermediate pCO2 levels, orbital configurations and ice-sheet geometries could help reduce this uncertainty by more accurately capturing the non-linear response of many climatological parameters to these forcings and boundary conditions.

Despite the large uncertainty, our results suggest that CO2 concentrations during this warm time interval have not been significantly higher than present-day (∼400 ppmv) values, in contrast to some of the proxy results. Comparing our Pliocene pCO2 reconstruction to those by van de Wal et al. (2011) and Stap et al. (2016), our model shows stronger variability on the 104-year timescale. In the long term, our model generally shows pCO2 levels for warm climates that are higher than those of van de Wal et al. (2011) but lower than those of Stap et al. (2016). For colder climates, our pCO2 is generally higher than Stap et al. (2016) and not clearly higher or lower than van de Wal et al. (2011). Given the level of disagreement between the different proxy-based reconstructions, it is not possible to assess the validity of the different model-based reconstructions relative to each other. However, based on the ability of the different models to reproduce the EPICA pCO2 record, assigning more confidence to the reconstruction presented here is justified.

Berends et al. (2018) provide a detailed discussion of the various advantages and disadvantages of the matrix method with respect to other methods of model forcing and coupling. Non-linear feedbacks of a growing ice sheet on the local and global climate, such as changes in atmospheric stationary waves, are not properly captured by this model set-up, although the inclusion of more GCM snapshots for intermediate-sized ice sheets should make the behaviour of the model more realistic in this respect. As a result, the inception of the last glacial cycle (100–80 kyr ago; Figs. 6, 7), is now also satisfyingly resolved in terms of temperature and sea-level drop though the decrease in pCO2 seems stronger than suggested by the ice-core record.

A drawback of the matrix method used here is that ocean temperature, required for calculating subshelf melt, is not included as a data field in the GCM snapshots. Instead, subshelf melt is calculated based on a combination of the temperature-based formulation by Martin et al. (2011) and the glacial–interglacial parameterisation by Pollard and DeConto (2009), tuned by de Boer et al. (2013) to produce realistic present-day Antarctic shelves and grounding lines. Although Berends et al. (2018) show that this set-up performs well when simulating colder-than-present climates, this is not necessarily a priori true for warmer climates, where the ice shelves are expected to retreat or even disintegrate. A more elaborate parameterisation based on GCM-calculated ocean temperatures can be expected to produce more reliable results.

Similarly, the effect of changes in insolation upon surface temperature is not well constrained. The climate matrix proposed by Berends et al. (2018) uses a parameterisation based on the locally absorbed insolation. While this allows changes in prescribed insolation to affect climate by changing the relative weights assigned to the different GCM snapshots in the climate matrix, the different GCM snapshots used in the current version of the climate matrix were all forced with the same 3.3 Ma reconstruction by Laskar et al. (2004). Expanding the climate matrix with additional GCM snapshots for different orbital parameters, along the lines of Prescott et al. (2014, 2018), would make the relation between insolation and surface temperature more explicit. We believe this could possibly lead to a further retreat of the East Antarctic ice sheet during warm periods. Another possible hindrance to significant retreat of the Antarctic ice sheet in our simulations is the lack of explicit grounding-line physics and relatively low model resolution, both of which have been shown to be required for accurate simulations of grounding-line retreat (Schoof, 2007; Gladstone et al., 2012; Leguy et al., 2014). Instead, ANICE calculates sheet and shelf ice velocities using the SIA and SSA, respectively, and add these together, without additional grounding-line parameterisations.

An additional source of uncertainty in our reconstruction is the palaeotopography of the period. Although we did include the Barents Sea erosion reversal by Butt et al. (2002) and its climate “fingerprint” as provided by Hill (2015) in our model, several other regions where ice may have existed during MIS M2 are suspected to have had a different topography – the Canadian Archipelago has been suggested to have been still one unbroken landmass which only formed later through erosion by ice during the Pleistocene glaciations (Dowsett et al., 2016), the Hudson Bay was likely not yet submerged (present today mostly due to remaining isostatic depression from the Laurentide ice sheet (Dowsett et al., 2016; Raymo et al., 2011). Similarly, based on the Eocene–Oligocene transition (34 Myr ago) palaeotopography reconstruction by Wilson and Luyendyk (2009), it is possible that, even during the Pliocene, West Antarctica was still mostly dry land (mostly submerged today due to erosion by ice and isostatic depression) and the Filchner–Ronne and Ross seas were significantly deeper (shallowed by ice-eroded sediment from West Antarctica). Although such changes in topography would likely have changed the evolution of the ice sheets, preliminary experiments for the Barents Sea showed that including the topography change without its GCM-calculated effect on climate resulted in a strong overestimation of ice volume, mostly because applying the present-day sea climate to the newly exposed high-latitude landmass resulted in a strongly positive mass balance even with pCO2 above 400 ppmv. Since no studies investigating the effects of these other topography changes on local and global climate are available yet, we did not include these changes in our study. Future work might be focussed on reinvestigating these effects once results from new GCM simulations with these topography changes become available.

Considering the results from the comparison of our model output to the available proxy data and the different uncertainties and caveats in our results, we believe our results could be of added value to future iterations of PlioMIP, to be used, for example, as boundary conditions for new GCM snapshots or even transient simulations.

Code and data availability

The reconstructed records of pCO2, global mean sea level and benthic δ18O, as well as NetCDF files containing ice thickness, bedrock topography and annual mean surface temperature and precipitation for all four ice-sheet model regions during MIS M2, MIS KM5c and MIS KM3 are available online in the supplementary information at (Berends et al., 2019).

Author contributions

CJB, BdB and RSWvdW designed the study. AMD and DJH provided data from their own studies. CJB created the model set-up and carried out the simulations, with support from BdB and 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.


Model runs were performed on the LISA computer cluster; we would like to acknowledge SurfSARA Computing and Networking Services for their support.

Financial support

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, grant no. 024.002.001). Bas de Boer was funded by NWO Earth and Life Sciences (ALW), project 863.15.019. This work was sponsored by NWO Exact and Natural Sciences for the use of supercomputer facilities. Aisling M. Dolan acknowledges funding from the European Commission under the European Union's Seventh Framework Programme (FP7/2007–2013)/ERC grant agreement no. 278636.

Review statement

This paper was edited by Ran Feng and reviewed by two anonymous referees.


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

Bachem, P. E., Risebrobakken, B., De Schepper, S., and McClymont, E. L.: Highly variable Pliocene sea surface conditions in the Norwegian Sea, Clim. Past, 13, 1153–1168,, 2017. 

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,, 2013. 

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

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,, 2018. 

Berends, C. J., de Boer, B., Dolan, A. M., Hill, D. J., and van de Wal, R. S. W.: Berends_etal_2019_GMD_supplement [Data set], Zenodo,, 2019. 

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. 

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. 

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. 

Bragg, F. J., Lunt, D. J., and Haywood, A. M.: Mid-Pliocene climate modelled using the UK Hadley Centre Model: PlioMIP Experiments 1 and 2, Geosci. Model Dev., 5, 1109–1125,, 2012. 

Burke, K. D., Williams, J. W., Chandler, M. A., Haywood, A. M., Lunt, D. J., and Otto-Bliesner, B. L.: Pliocene and Eocene provide best analogs for near-future climates, P. Natl. Acad. Sci. USA, 115, 13288–13293, 2018. 

Butt, F. A., Drange, H., Elverhøj, A., Otterå, O. H., and Solheim, A.: Modelling Late Cenozoic isostatic elevation changes in the Barents Sea and their implications for oceanic and climatic regimes: preliminary results, Quaternary Sci. Rev. 21, 1643–1660, 2002. 

Church, J. A., Clark, P. U., Cazenave, A., Gregory, J. M., Jevrejeva, S., Levermann, A., Merrifield, M. A., Milne, G. A., Nerem, R. S., Nunn, P. D., Payne, A. J., Pfeffer, W. T., Stammer, D., and Unnikrishnan, A. S.: Sea-Level Change, Climate Change 2013: The Physical Science Basis. Contribution of Working Group 1 to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, UK and New York, NY, USA, 2013. 

Cramer, B. S., Toggweiler, J. R., Wright, J. D., Katz, M. E., and Miller, K. G.: Ocean overturning since the Late Cretaceous: Inferences from a new benthic foraminiferal isotope compilation, Paleoceanography, 24, PA4216,, 2009. 

Csank, A. Z., Tripati, A. K., Patterson, W. P., Eagle, R. A., Rybczynski, N., Ballantyne, A., and Eiler, J. M.: Estimates of Arctic land surface temperatures during the early Pliocene from two novel proxies, Earth Planet. Sc. Lett., 304, 291–299, 2011a. 

Csank, A. Z., Patterson, W. P., Eglington, B. M., Rybczynski, N., and Basinger, J. F.: Climate variability in the Early Pliocene Arctic: Annually resolved evidence from stable isotope values of sub-fossil wood, Ellesmere Island, Canada, Palaeogeogr. Palaeocl., 308, 339–349, 2011b. 

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., 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., Stocchi, P., and van de Wal, R. S. W.: A fully coupled 3-D ice-sheet–sea-level model: algorithm and applications, Geosci. Model Dev., 7, 2141–2156,, 2014. 

de Boer, B., Dolan, A. M., Bernales, J., Gasson, E., Goelzer, H., Golledge, N. R., Sutter, J., Huybrechts, P., Lohmann, G., Rogozhina, I., Abe-Ouchi, A., Saito, F., and van de Wal, R. S. W.: Simulating the Antarctic ice sheet in the late-Pliocene warm period: PLISMIP-ANT, an ice-sheet model intercomparison project, The Cryosphere, 9, 881–903,, 2015. 

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. 

de Schepper, S., Gibbard, P. L., Salzmann, U., and Ehlers, J.: A global synthesis of the marine and terrestrial evidence for glaciation during the Pliocene epoch, Earth Sci. Rev., 135, 83–102, 2014. 

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. 

Dowsett, H. J., Robinson, M. M., and Foley, K. M.: Pliocene three-dimensional global ocean temperature reconstruction, Clim. Past, 5, 769–783,, 2009. 

Dowsett, H. J., Robinson, M. M., Haywood, A. M., Salzmann, U., Hill, D. J., Sohl, L. E., Chandler, M. A., Williams, M., Foley, K. M., and Stoll, D. K.: The PRISM3D paleoenvironmental reconstruction, Stratigraphy 7, 123–139, 2010. 

Dowsett, H. J., Foley, K. M., Stoll, D. K., Chandler, M. A., Sohl, L. E., Bentsen, M., Otto-Bliesner, B. L., Bragg, F. J., Chan, W.-L., Contoux, C., Dolan, A. M., Haywood, A. M., Jonas, J. A., Jost, A., Kamae, Y., Lohmann, G., Lunt, D. J., Nisancioglu, K. H., Abe-Ouchi, A., Ramstein, G., Riesselman, C. R., Robinson, M. M., Rosenbloom, N. A., Salzmann, U., Stepanek, C., Strother, S. L., Ueda, H., Yan, Q., and Zhang, Z.: Sea Surface Temperature of the mid-Piacenzian Ocean: A Data-Model Comparison, Nature Scientific Reports, 3,, 2013. 

Dowsett, H., Dolan, A., Rowley, D., Moucha, R., Forte, A. M., Mitrovica, J. X., Pound, M., Salzmann, U., Robinson, M., Chandler, M., Foley, K., and Haywood, A.: The PRISM4 (mid-Piacenzian) paleoenvironmental reconstruction, Clim. Past, 12, 1519–1538,, 2016. 

Dwyer, G. S. and Chandler, M. A.: Mid-Pliocene sea level and continental ice volume based on coupled benthic Mg/Ca palaeotemperatures and oxygen isotopes, Philos. T. R. Soc. A, 367, 157–168, 2009. 

Fletcher, T., Feng, R., Telka, A. M., Matthews, J. V. J., and Ballantyne, A.: Floral dissimilarity and the Influence of Climate in the Pliocene High Arctic: Biotic and Abiotic Influences on Five Sites on the Canadian Arctic Archipelago, Frontiers in Ecology and Evolution, 5,, 2017. 

Gao, C., McAndrew, J. H., Wang, X., Menzies, J., Turton, C. L., Wood, B. D., Pei, J., and Kodors, C.: Glaciation of North America in the James Bay Lowland, Canada, 3.5 Ma, Geology, 40, 975–978, 2012. 

Gladstone, R. M., Payne, A. J., and Cornford, S. L.: Resolution requirements for grounding-line modelling: sensitivity to basal drag and ice-shelf buttressing, Ann. Glaciol., 53, 97–105, 2012. 

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. 

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

Haywood, A. M., Dowsett, H. J., Otto-Bliesner, B., Chandler, M. A., Dolan, A. M., Hill, D. J., Lunt, D. J., Robinson, M. M., Rosenbloom, N., Salzmann, U., and Sohl, L. E.: Pliocene Model Intercomparison Project (PlioMIP): experimental design and boundary conditions (Experiment 1), Geosci. Model Dev., 3, 227–242,, 2010. 

Haywood, A. M., Dowsett, H. J., Robinson, M. M., Stoll, D. K., Dolan, A. M., Lunt, D. J., Otto-Bliesner, B., and Chandler, M. A.: Pliocene Model Intercomparison Project (PlioMIP): experimental design and boundary conditions (Experiment 2), Geosci. Model Dev., 4, 571–577,, 2011. 

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,, 2013a. 

Haywood, A. M., Dolan, A. M., Pickering, S. J., Dowsett, H. J., McClymont, E. L., Prescott, C. L., Salzmann, U., Hill, D. J., Hunter, S. J., Lunt, D. J., Pope, J. O., and Valdes, P. J.: On the identification of a Pliocene time slice for data-model comparison, Philos. T. R. Soc. A, 371,, 2013b. 

Hill, D. J.: The non-analogue nature of Pliocene temperature gradients, Earth Planet. Sc. Lett., 425, 232–241, 2015. 

Imbrie, J., Hays, J. D., Martinson, D. G., McIntyre, A., Mix, A. C., Morley, J. J., Pisias, N. G. Prell, W. L., and Shackleton, N. J.: The orbital theory of Pleistocene climate: Support from a revised chronology of the marine δ18O record, in: Milankovitch and climate: understanding the response to astronomical forcing, Part 1, 269–305, 1984. 

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 Millenioal 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. 

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. 

Leguy, G. R., Asay-Davis, X. S., and Lipscomb, W. H.: Parameterization of basal friction near grounding lines in a one-dimensional ice sheet model, The Cryosphere, 8, 1239–1259,, 2014. 

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

Lunt, D. J., Haywood, A. M., Foster, G. L., and Stone, E. J.: The Arctic cryosphere in the Mid-Pliocene and the future, Philos. T. R. Soc. A, 367, 49–67, 2009. 

Lunt, D. J., Haywood, A. M., Schmidt, G. A., Salzmann, U., Valdes, P. J., and Dowsett, H. J.: Earth system sensitivity inferred from Pliocene modelling and data, Nat. Geosci., 3, 60–64,, 2010. 

Lunt, D. J., Haywood, A. M., Schmidt, G. A., Salzmann, U., Valdes, P. J., Dowsett, H. J., and Loptson, C. A.: On the causes of mid-Pliocene warmth and polar amplification, Earth Planet. Sc. Lett., 321–322, 128–138, 2012. 

Lüthi, D., Le Floch, M., Bereiter, B., Blunier, T., Barnola, J.-M., Siegenthaler, U., Raynaud, D., Jouzel, J., Fisher, H., Kawamura, K., and Stocker, T. F.: High-resolution carbon dioxide concentration record 650,000–800,000 years before present, Nature Letters, 452, 379–382, 2008. 

Ma, Y., Gagliardini, O., Ritz, C., Gillet-Chaulet, F., Durand, G., and Montagnat, M.: Enhancement factors for grounded ice and ice shelves inferred from an anisotropic ice-flow model, J. Glaciol., 56, 805–812, 2010. 

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. 

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

McKay, R., Naish, T., Carter, L., Riesselman, C. R., Dunbar, R., Sjunneskog, C., Winter, D., Sangiorgi, F., Warren, C., Pagani, M., Schouten, S., Willmott, V., Levy, R., DeConto, R., and Powell, R. D.: Antarctic and Southern Ocean influences on Late Pliocene global cooling, P. Natl. Acad. Sci. USA, 109, 6423–6428, 2012. 

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., 99–116, Kluwer Acad., Dordrecht, the Netherlands, 1987. 

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

Naish, T. and Wilson, G.: Constraints on the amplitude of mid-Pliocene (3.6–2.4 Ma) eustatic sea-level fluctuations from the New Zealand shallow-marine sediment record, Philos. T. R. Soc. A, 367, 169–187, 2009. 

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, Ceeplusplusand 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. 

Niu, L., Lohmann, G., Hinck, S., and Gowan, E. J.: Sensitivity of atmospheric forcing on Northern Hemisphere ice sheets during the last glacial-interglacial cycle using output from PMIP3, Clim. Past Discuss.,, 2017. 

Peltier, W. R.: Global Glacial Isostasy and the Surface of the Ice-Age Earth: The ICE-5G (VM2) Model and GRACE, Annu. Rev. Earth Planet. Sc., 32, 111–149, 2004. 

Pollard, D.: A retrospective look at coupled ice sheet-climate modelling, 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. 

Prescott, C. L., Haywood, A. M., Dolan, A. M., Hunter, S. J., Pope, J. O., and Pickering, S. J.: Assessing orbitally-forced interglacial climate variability during the mid-Pliocene Warm Period, Earth Planet. Sc. Lett., 400, 261–271, 2014. 

Prescott, C. L., Dolan, A. M., Haywood, A. M., Hunter, S. J., and Tindall, J. C.: Regional climate and vegetation response to orbital forcing within the mid-Pliocene Warm Period: A study using HadCM3, Global Planet. Change, 161, 231–243, 2018. 

Raymo, M. E., Mitrovica, J. X., O'Leary, M. J., DeConto, R. M., and Hearty, P. J.: Departures from eustasy in Pliocene sea-level records, Nat. Geosci., 4, 328–332, 2011. 

Salzmann, U., Dolan, A. M., Haywood, A. M., Chan, W.-L., Voss, J., Hill, D. J., Abe-Ouchi, A., Otto-Bliesner, B. L., Bragg, F. J., Chandler, M. A., Contoux, C., Dowsett, Harry Jupiter amd Jost, A., Kamae, Y., Lohmann, G., Lunt, D. J., Pickering, S. J., Pound, M. J., Ramstein, G., Rosenbloom, N. A., Sohl, L., Stepanek, C., Ueda, H., and Zhang, Z.: Challenges in quantifying Pliocene terrestrial warming revealed by data-model discord, Nat. Clim. Change 3, 969–974, 2013. 

Schoof, C.: Ice sheet grounding line dynamics: Steady states, stability, and hysteresis, J. Geophys. Res., 112, F03S28,, 2007. 

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 sensitivty to ice-sheet forcing over the last 120 kyr, Quaternary Sci. Rev., 29, 43–55, 2010. 

Smith, Y. M., Hill, D. J., Dolan, A. M., Haywood, A. M., Dowsett, H. J., and Risebrobakken, B.: Icebergs in the Nordic Seas Throughout the Late Pliocene, Paleoceanography and Paleoclimatology, 33, 318–335, 2018. 

Sohl, L., Chandler, M. A., Schmunk, R. B., Mankoff, K., Jonas, J. A., Foley, K. M., and Dowsett, H. J.: PRISM3/GISS topographic reconstruction, U.S. Geol. Surv. Data Series, 419, 6, 2009. 

Solomon, S., Qin, D., Manning, M., Chen, Z., Marquis, M., and Averyt, K. B.: Climate Change 2007: the Physical Science Basis, Contribution Of Working Group 1 To The Fourth Assessment Report Of The Intergovernmental Panel On Climate Change, Cambridge University Press, Cambridge, UK and New York, NY, USA, 2007. 

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. 

Swann, G. E. A., Kendrick, C. P., Dickson, A. J., and Worne, S.: Late Pliocene Marine pCO2 Reconstructions From the Subarctic Pacific Ocean, Paleoceanography and Paleoclimatology, 33, 457–469, 2018.  

Tan, N., Ramstein, G., Dumas, C., Contoux, C., Ladant, J.-B., Sepulchre, P., Zhang, Z., and de Schepper, S.: Exploring the MIS M2 glaciation occurring during a warm and high atmospheric CO2 Pliocene background climate, Earth Planet. Sc. Lett., 472, 266–276, 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. 

Wilson, D. S. and Luyendyk, B. P.: West Antarctic paleotopography estimated at the Eocene-Oligocene climate transition, Geophys. Res. Lett., 36, L16302,, 2009. 

Winkelmann, R., Martin, M. A., Haseloff, M., Albrecht, T., Bueler, E., Khroulev, C., and Levermann, A.: The Potsdam Parallel Ice Sheet Model (PISM-PIK) – Part 1: Model description, The Cryosphere, 5, 715–726,, 2011. 

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,, 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
The Late Pliocene, 3.65–2.75 million years ago, is the most recent period in Earth's history that was warmer than the present. This makes it interesting for climatological research, because it provides a possible analogue for the near future. We used a coupled ice-sheet–climate model to simulate the behaviour of these systems during this period. We show that the warmest moment saw a sea-level rise of 8–14 m, with a CO2 concentration of 320–400 ppmv.