# Dynamic Greenland ice sheet driven by pCO2 variations across the Pliocene Pleistocene transition

Ning Tan, Jean-Baptiste Ladant, Gilles Ramstein, Christophe Dumas, Paul Bachem, Eystein Jansen

► To cite this version:

Ning Tan, Jean-Baptiste Ladant, Gilles Ramstein, Christophe Dumas, Paul Bachem, et al.. Dynamic Greenland ice sheet driven by pCO2 variations across the Pliocene Pleistocene transition. Nature Communications, 2018, 9, pp.4755. 10.1038/s41467-018-07206-w. hal-01945992

HAL Id: hal-01945992

https://hal.sorbonne-universite.fr/hal-01945992

Submitted on 5 Dec 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ARTICLE

DOI: 10.1038/s41467-018-07206-w

OPEN

# Dynamic Greenland ice sheet driven by pCO₂ variations across the Pliocene Pleistocene transition

Ning Tan<fnref n="1" />,², Jean-Baptiste Ladant<fnref n="1" />,³,⁴, Gilles Ramstein<fnref n="1" />, Christophe Dumas<fnref n="1" />, Paul Bachem⁵ & Eystein Jansen⁶

It is generally considered that the perennial glaciation of Greenland lasting several orbital cycles began around 2.7 Ma along with the intensification of Northern Hemisphere glaciation (NHG). Both data and model studies have demonstrated that a decline in atmospheric pCO₂ was instrumental in establishing a perennial Greenland ice sheet (GrIS), yet models have generally used simplistic pCO₂ constraints rather than data-inferred pCO₂ evolution. Here, using a method designed for the long-term coupling of climate and cryosphere models and pCO₂ scenarios from different studies, we highlight the pivotal role of pCO₂ on the GrIS expansion across the Plio-Pleistocene Transition (PPT, 3.0-2.5 Ma), in particular in the range between 280 and 320 ppm. Good qualitative agreement is obtained between various IRD reconstructions and some of the possible evolutions of the GrIS simulated by our model. Our results underline the dynamism of the GrIS waxing and waning under pCO₂ levels similar to or lower than today, which supports recent evidence of a dynamic GrIS during the Plio-Pleistocene.

<footnote n="1">¹Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, 91191 Gif-sur-Yvette, France. ²Key Laboratory of Cenozoic Geology and Environment, Institute of Geology and Geophysics, Chinese Academy of Sciences, 100029 Beijing, China. ³Département de Géosciences, École Normale Supérieure, Paris 75005, France. ⁴LMD/IPSL, CNRS/ENS/UPMC, Ecole Polytechnique, Paris 75005, France. ⁵Uni Research Climate, Bjerknes Centre for Climate Research, Jahnebakken 5, 5007 Bergen, Norway. ⁶Dep. of Earth Science, University of Bergen, Bjerknes Center for Climate Research, Jahnebakken 5, 5007 Bergen, Norway. Correspondence and requests for materials should be addressed to N.T. (email: ning.tan@mail.iggcas.ac.cn)</footnote>
The long-term trend that led to the initiation of the cyclic Northern Hemisphere glaciations (NHG) can be dated from 3.6 Ma¹. A set of reconstructed SST records²⁻⁴ spread across the globe describes a gradual cooling trend for the middle to late Pliocene and early Pleistocene (3.6–2.2 Ma), associated with a progressive increase in high δ¹⁸O peaks of benthic foraminifera⁵, reflecting both lower deep-water temperatures and increased ice volume. Terrestrial data from Lake El'gygytgyn are also consistent with long-term ocean cooling⁶. During this period, the first marked peak of IRD was found around 3.4–3.3 Ma (MIS MG2 and MIS M2) East of Greenland (ODP site 907⁷) and is interpreted as evidence of an important glacial event occurring before the establishment of the cyclic NHG around 2.7 Ma⁷,⁸. Significant amounts of IRD were recovered in different ocean drilling sites in the North Atlantic during the PPT⁹. These deposited IRD likely originated from Greenland, Iceland, North America and Scandinavia, suggesting large land ice expansions to the coast during this interval⁷,⁹. In addition, records based on volcanic and sedimentary facies from Iceland also indicate the onset of a large glaciation after 2.6 Ma¹⁰,¹¹. The onset of a perennial GrIS and of the cyclic NHG at the end of the PPT stands out as a tipping point in the climate evolution of the Earth. Indeed, it marks the beginning of a low pCO₂ world with perennial ice sheets in both hemispheres, an infrequent occurrence in the Earth's history¹², thereby creating specific geologic and climatic conditions allowing the development of glacial/interglacial cycles.

Greenland, however, may have experienced waxing and waning of ice before the intensification of NHG, as suggested by Eocene, late Miocene and early Pliocene IRD records¹³⁻¹⁵. In particular, the last large Northern Hemisphere glaciation prior to the major intensification at 2.7 Ma occurred during MIS-M2 (3.312–3.264 Ma)¹,¹⁶⁻¹⁸. This 50 kyrs glaciation was followed by the well-established Mid-Pliocene Warm Period (MPWP) from 3.3 Ma to 3.0 Ma for which numerous observational and modelling studies performed in the framework of the PlioMIP project¹⁹ have demonstrated that the warmer conditions led to a reduced GrIS²⁰⁻²² (Supplementary Fig. 1). The subsequent build-up of a perennial GrIS across the PPT remains poorly constrained from a spatio-temporal point of view due to its minor contribution to the signal in global benthic foraminiferal δ¹⁸O records and the paucity of direct geological data associated with the GrIS expansion. Models have the potential to provide valuable insights into the GrIS evolution but its transient nature remains complicated to explore with fully coupled Global Climate Models (GCM). Most climate/cryosphere studies on the GrIS onset during the PPT using GCM have generally produced snapshot climatic simulations subsequently used to force an ice sheet model (e.g. ref. ²³). For instance, Lunt et al.²³ demonstrated, based on a series of fully coupled GCM snapshot experiments with different forcing factors, that pCO₂ decline was the major driver of the GrIS glaciation. However, this result was obtained from equilibrium simulations at 2.7 Ma and the underlying GCM simulations included a pre-existent Pliocene GrIS²¹. A more recent study, also using snapshot simulations at 2.7 Ma, demonstrated that pCO₂ values had to remain low to counterbalance the increasing summer insolation in order to maintain the GrIS after the initial onset²⁴. In contrast, low resolution, conceptual and intermediate complexity models have been used to perform transient long-term experiments²⁵⁻²⁷, but these models remain simplified with respect to many processes and do not have the spatial resolution to focus specifically on Greenland. Willeit et al.²⁸ have recently worked on increasing the resolution but their simulations of the GrIS across the PPT still used pre-defined pCO₂ forcings. Therefore, although important insights into the dynamics of the cryosphere across the PPT have emerged from

these studies, simulation of the GrIS evolution across PPT with a forward physically-based model driven by realistic forcings remains a major challenge. Here, we use a recent numerical interpolation method²⁹ (see Methods), which couples climate simulations obtained with the fully coupled IPSL-CM5A model³⁰ and the 15-km resolution version of the ice sheet model GRISLI³¹, in order to investigate the transient evolution of the GrIS across the PPT. The major advantage of our approach is the ability to directly test the response of the climate-ice sheet system to different pCO₂ evolution scenarios (see Methods) in order to define plausible pCO₂ scenarios that led to the GrIS inception and variability. We thus apply various reconstructions of pCO₂ evolution across the PPT from both proxy records and inverse modelling studies (e.g. ref. ²⁷,²⁸,³²⁻³⁵), as well as constant pCO₂ evolutions.

## Results

**GrIS sensitivity to constant pCO₂ scenarios.** We first present six ice sheet experiments with constant pCO₂ values ranging from 220 to 405 ppmv using realistic orbital variations³⁶ (Supplementary Fig. 2b). As shown in Fig. 1a, the evolution of the GrIS clearly shows that below 280 ppmv of pCO₂, it is possible to trigger and maintain a large perennial ice sheet over Greenland even during intervals of strong summer insolation, in particular around 2.6 Ma, whereas pCO₂ higher than 320 ppmv prevent the GrIS from remaining perennial across the whole PPT period. Between 240 and 360 ppmv, the sensitivity of the GrIS volume to orbital variations is considerable, illustrating the role of the complex interplay between pCO₂ and orbital variations in the dynamics of the GrIS, especially after 2.7 Ma when insolation becomes highly variable. Finally for pCO₂ levels above 320 ppm, the ice sheet extent remains limited across the whole 3–2.5 Ma interval. The results of these constant simulations are in good agreement but go one step beyond previous transient experiments at constant pCO₂³⁷ in that we demonstrate that the GrIS possesses a dynamism on orbital timescales across the PPT for only a narrow range of atmospheric pCO₂ concentrations.

**GrIS sensitivity to the initial boundary conditions.** Existing knowledge indicates that the GrIS largely retreated during the MPWP, but the exact configuration of the Greenland ice sheet and its topography remain highly uncertain¹⁷. In order to investigate the impact of different initial states of the GrIS on its subsequent evolution, we test two different initial boundary conditions: an ice-free Greenland (Supplementary Fig. 1c) and the modern GrIS (Supplementary Fig. 1a). We perform constant pCO₂ simulations at three different levels (220 ppmv, 320 ppmv, 405 ppmv). Our results show that for extreme pCO₂ values (220 or 405 ppm) the simulated ice sheet evolution does not depend upon the initial configuration of the ice sheet at 3.0 Ma (Fig. 1b). For pCO₂ levels closer to the modelled threshold for glaciation (320 ppmv), orbital and pCO₂ forcings have comparable influence on the simulated GrIS volume but these simulations demonstrate that the initial configuration of the ice sheet also controls the GrIS response. Indeed, bedrock vertical motion driven by the presence or absence of ice can alter topography for tens or hundreds of meters, which, at pCO₂ levels close to the threshold, may lead to differences in ice volume and extent. However, the general evolution of the GrIS volume remains similar in the scenario starting with an ice-free Greenland and in that starting from a full ice sheet. More importantly, during major deglaciation episodes driven by high summer insolation, the modelled GrIS in both cases recedes to identical configurations. As there is evidence for a much reduced, perhaps even non-existent, GrIS during interglacials of the MPWP (e.g., ref. ³⁸,³⁹), we initialize the ice sheet
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Sensitivity tests with constant pCO₂ scenarios. a Simulated GrIS volume (ISV on the figure) with various constant pCO₂ concentrations (220 ppmv, 240 ppmv, 280 ppmv, 320 ppmv, 360 ppmv, 405 ppmv; so that ISV_220 indicates the simulated GrIS with constant pCO₂ of 220 ppmv). b Simulated GrIS volume with the three pCO₂ concentrations (220 ppmv, 320 ppmv and 405 ppmv under two different initial GrIS configurations: the present-day full GrIS or an ice-free Greenland. Light blue dash lines represent the summer mean insolation at 65 °N</figcaption></figure>

model with an ice-free Greenland. Indeed, the starting date of our simulation (3 Ma) is preceded by a very high summer insolation maximum (Supplementary Fig. 2b).

**GrIS evolution using pCO₂ reconstructions.** In the following, we force our set-up with published reconstructions of atmospheric pCO₂, derived from recent proxy reconstructions³²⁻³⁴ or from inverse modelling²⁷,²⁸,³⁵, in order to compare the modelled GrIS evolution to previous work and to IRD records across the PPT.

We employ recent pCO₂ reconstructions based on alkenones and boron measurements³²⁻³⁴ (teal lines on Fig. 2). These records depict different pCO₂ evolution scenarios from a mostly linear and low-resolution decrease for the older record³⁴ to high-resolution and high pCO₂ variations in the most recent record³³. When forced by the estimated high pCO₂ levels of Seki et al.³⁴, the GrIS evolution confirms the major role of orbital variations when pCO₂ levels are confined to between ~280 ppmv and ~320 ppmv (Fig. 2a). While pCO₂ levels remain above 350 ppm, even large orbital variations do not significantly affect the extension of the GrIS (3.0 to 2.9 Ma). During the 2.9–2.7 Ma interval, the low variability in insolation prevents the GrIS from growing in spite of pCO₂ levels dropping to 300 ppm. The simulated GrIS evolution then shows a significant increase at 2.7 Ma and pronounced orbital scale variability between 2.7 and 2.5 Ma, but importantly, there is no full retreat of the GrIS during the insolation maximum after 2.6 Ma because this maximum is associated with pCO₂ levels below 300 ppm (Fig. 2a). In contrast, the Seki et al.³⁴ estimates of low pCO₂ levels generate a large perennial GrIS as early as 2.9 Ma triggered by the combination of an insolation minimum and low pCO₂. A similar modelled evolution of the GrIS is obtained when forced by the pCO₂ record of Bartoli et al.³². Regardless of the uncertainties on the absolute values of pCO₂ concentration, the low levels in this record force an early onset of a perennial GrIS (Fig. 2b). However, because of their low resolution, these two records do not show any pCO₂ variability on the ~10 kyrs timescale contrary to that of Martinez-Boti et al.³³, which demonstrates that the PPT pCO₂ evolution is in fact much more variable than previously thought. The uncertainties associated with the estimates of the Martinez-Boti et al.³³ record show that a completely different evolution of the GrIS can be simulated (Fig. 2c). The high (low) estimates depict

an evolution close to that driven by the high (low) estimates of Seki et al.³⁴, because the pCO₂ levels during times of insolation extremes (3.0–2.9 Ma and 2.7–2.6 Ma) are similar. However, the simulation forced with the mean pCO₂ estimates of Martinez-Boti et al.³³ shows an early attempt of GrIS expansion (~2.98–2.94 Ma) before a progressive onset from 2.8 Ma onwards, with a significant increase in GrIS volume around 2.72 Ma and 2.6 Ma (Fig. 3) because of the combination of low insolation and low pCO₂ levels. After 2.6 Ma, the GrIS rapidly melts down to small ice caps in the southern and southeastern margins (Figs 2c and 3) because of a simultaneous increase in summer insolation and pCO₂. The ice sheet expansion then resumes at ~2.55 Ma (A complete GrIS evolution with the mean pCO₂ estimates of Martinez-Boti et al.³³ is shown in the Supplementary Movie 1). Interestingly, the pCO₂ records of Martinez-Boti et al.³³ generate a sharp GrIS volume decrease after 2.6 Ma, regardless of the uncertainties and earlier shape of the expansion of the ice sheet because of an increase in pCO₂ values during the insolation maximum. This seems contradictory to the presence of large and perennial ice sheets from 2.7 Ma onwards but comparison with IRD records⁷,⁹ does not invalidate a scenario with a widescale deglaciation after 2.6 Ma.

When forced by the monotonously decreasing and obliquity modulated pCO₂ scenario defined by Willeit et al.²⁸, our simulated GrIS shows a similar evolution to that of Willeit et al.²⁸ (Fig. 4), in particular in the large increases in GrIS volume reproduced at 2.7 Ma and 2.55 Ma, as well as the large decrease at 2.6 Ma, although our GrIS volume displays less variability because of a lower sensitivity to the orbital forcing in our fully coupled climate-ice sheet model. Inverse modelling studies²⁷,³⁵ have also provided potential pCO₂ reconstructions for this interval (Supplementary Fig. 3). However, the low pCO₂ concentrations throughout most of these reconstructions lead to a large perennial GrIS even during the well-established MPWP warming period¹⁹, in contrast to evidence for a much reduced GrIS³⁸,³⁹.

## Discussion

In this work, we restrict high-resolution ice sheet modelling to scenarios of the GrIS evolution across the PPT although much evidence for ice growth outside Greenland has been reported during this period⁸,⁴⁰,⁴¹. This decision is warranted because Greenland is one of the primary location for ice nucleation in the
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Simulated GrIS volume evolution based on different pCO₂ records. **a** Seki et al.³⁴; **b** Bartoli et al.³²; **c** Martinez et al.³³. Dashed light blue line represents the boreal summer insolation at 65°N. Orange lines are the simulated GrIS volumes based on pCO₂ records and their uncertainties, represented by the solid teal lines and shaded teal areas. **a** Alkenone-based pCO₂ record of Seki et al.³⁴, determined using size-corrected ε₃₇,₂ values for the modern range of b-values³⁴. The solid (dashed) orange line is the simulated GrIS volume obtained with the solid (dashed) pCO₂ line. **b** Boron isotopes-based pCO₂ record of Bartoli et al.³², with a 2 σ error range. The solid orange line is the simulated GrIS volume with mean pCO₂ values and shaded orange lines are the GrIS volumes obtained with low and high extremes of the pCO₂ uncertainties. **c** As in **b** but for the boron isotopes-based pCO₂ record of Martinez-Boti et al.³³, with an uncertainty range corresponding to 68% of 10,000 Monte Carlo simulations with full propagation of uncertainties³³. Circled numbers at 2.9 Ma, -2.8 Ma, -2.72 Ma and 2.6 Ma correspond to the ice sheet snapshot displayed in Fig. 3</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Greenland ice sheet thickness snapshots in the modelled GrIS evolution. These snapshots are taken from the simulated GrIS evolution based on the pCO₂ record of Martinez-Boti et al.³³ (mean) at 2.9 Ma (**a**), 2.8 Ma (**b**), 2.72 Ma (**c**) and 2.6 Ma (**d**), as shown on Fig. 2c</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Simulated GrIS volume evolution with the pre-defined best scenario of pCO₂ from Willeit et al.²⁸. The brown dash line represents the GrIS evolution of Willeit et al.²⁸ under the same pCO₂ scenario. The dashed light blue line represents the boreal summer insolation at 65°N</figcaption></figure>

North Hemisphere¹⁵ and because there is currently no agreement on the possible extent and configuration of possible ice sheets outside Greenland. However, this restriction hampers comparisons with δ¹⁸O and sea level records because the maximum ice volume that can be accommodated over Greenland represents roughly 7 m of equivalent sea level, i.e. in the range of experimental and/or instrumental errors. Significant insights can however be gained from adjacent IRD records⁹ in order to further constrain the possible evolution of the GrIS and of pCO₂ across the PPT. It should be kept in mind that attempting to constrain the GrIS geometry and volume from IRD records remains speculative because the absence of IRDs does not necessarily correlate with the absence of GrIS and because IRD peaks may contain material derived from other sources than the melting Greenland icebergs. In addition, iceberg trajectories may change and the ambient temperature along the iceberg trajectory may influence the melt-out rates of the IRD contained in the icebergs. For instance, variations in North Atlantic SST and/or currents or changes in the sediment contents of the calved icebergs could explain the absence of IRD even if the GrIS maintains a significant iceberg discharge⁴¹,⁴². In addition, IRD peaks may represent an ice growth phase, during which increasing GrIS volume could lead to increased iceberg discharge, as well as an ice melt phase, during which the warmer climatic conditions could lead to enhanced ice sheet melting and, consequently, enhanced ice flux at the margin. Despite these uncertainties, IRD records offer a first-order insight into ice sheet dynamics and their relationship to orbital variations, which may aid in defining the most probable scenario for GrIS evolution. We used ODP Site 907⁷, IODP Site U1307⁴³ and ODP Site 611⁴¹ (Fig. 5. More details about these data can be found in Supplementary Note 1) because the IRD records from these sites can be confidently assumed to originate primarily from the GrIS⁴¹. Sites U1307 and 907 are located offshore from Greenland margins and show small but continuous IRD deposition from as early as 3 Ma (Fig. 5b), with the exception of a single peak at 2.92 Ma at Site 907. Around 2.7 Ma, several IRD peaks in both records suggest intensification of the iceberg discharge, broadly correlated to North Atlantic sea surface temperatures (SST) cooling events (Fig. 5c). Accordingly, the record from Site 611, located further south from Greenland, shows an absence of IRD deposition before 2.72 Ma and several peaks afterwards.

There is variable agreement between our reconstructed GrIS volume scenarios and inferred GrIS extension as implied by IRD records. For instance, the GrIS evolution forced by the pCO₂ record of Bartoli et al.³² agrees poorly with the IRD records

because it suggests a large to near-complete perennial GrIS as early as 2.98–2.95 Ma (Fig. 2b). Around 2.9 Ma, a large ice sheet covering the totality of the southern Greenland margin is present in the results forced by the pCO₂ records of Bartoli et al.³². However, there are not any particular changes above the 3.0–2.7 Ma background IRD values at Site U1307. At Site 907, an early increase of a large GrIS could explain the 2.92 IRD peak, but in this case the IRD deposition should have continued with higher levels relative to the pre-glaciation interval. In addition, the absence of any IRD at Site 611 at that time suggests that the GrIS was still of relatively limited size. In contrast to the pCO₂ record of Bartoli et al.³², the mean and high pCO₂ estimates of Martinez-Boti et al.³³ and the high estimate of Seki et al.³⁴ generate GrIS evolutions more consistent with the IRD records. The modelled GrIS using these reconstructions is limited to small ice caps on the southern and southeastern margins of Greenland during the interval 3.0–2.7 Ma interval, which agrees well with small but continuous IRD inputs at Site 907 and U1307 (Fig. 2a, c, Fig. 3 and Fig. 5a, b). The sharp increase in IRD deposition at 2.7 Ma is also accounted for by these GrIS evolution scenarios. Furthermore, the IRD records of Site 611⁴¹ displays four IRD peaks (at ~2.7 Ma, ~2.64 Ma, ~2.6 Ma and ~2.52 Ma) that are relatively well correlated in time with decreases in North Atlantic SSTs and large modelled GrIS variations (Fig. 5b, c), in particular for the Martinez-Boti et al.³³ mean and high pCO₂ scenario. Importantly, the large GrIS volume decrease after 2.6 Ma in these scenarios is not at odds with the IRD records, in that the minimal GrIS state still reaches the southern and southeastern margins of Greenland, allowing small but uninterrupted iceberg discharge at Site 907. These three pCO₂ scenarios (Martinez et al.³³ high, mean and Seki et al.³⁴ high pCO₂ estimates) therefore seem to correlate relatively well with IRD records because the GrIS expands preferentially from the south and east regions and intervals of large GrIS variations are synchronous with IRD peaks.

There is evidence however showing that after 2.7 Ma, ice sheets were not restricted to Greenland (e.g., ref. ⁴⁰,⁴¹). Although by model design ice growth in other regions of the Northern Hemisphere is not simulated, simple inferences of the impact of North American and Scandinavian ice sheets after 2.7 Ma would tend to suggest that the large GrIS retreat around 2.6 Ma might have been more limited. Indeed, the cooling effects of such ice masses on regional and global climate would presumably have partly counterbalanced the combined increase in pCO₂ and summer insolation. Still, the evidence for other NH ice sheets can tentatively be used to suggest that the GrIS evolution driven by the Martinez-Boti et al.³³ mean pCO₂ estimates is the most
<figure><img src="image_6.png" /><figcaption>**Fig. 5** The comparison between the simulated GrIS and adjacent IRD records. a Map showing simplified main ocean currents in NATL and the ocean deep drilling sites referenced in this study (Ocean Data View, Schlitzer, R., Ocean Data View, odv.awi.de, 2018). b Related IRD records in the North Atlantic regions: site 907⁷; site U1307⁴³; site 611⁴¹ and simulated GrIS volume with the mean and high pCO₂ estimates of Martinez-Boti et al.³³ and the high pCO₂ estimate of Seki et al.³⁴; c Reconstructed sea surface temperature at ODP site 982 by Lawrence et al.² and IODP site U1313 by Naafs et al.⁶³</figcaption></figure>

plausible pCO₂ scenario among those that are the most consistent with IRD records (i.e. Martinez-Boti et al.³³ mean and high, Seki et al.³⁴ high estimates), because it leads to a much larger GrIS between 2.8 and 2.6 Ma (Fig. 2c). The presence of ice sheets over North America and Scandinavia after 2.7 Ma is thus likely to be more consistent with a largely or fully glaciated Greenland because ice preferentially nucleates on Greenland before potentially expanding over other regions of the Northern Hemisphere.

In conclusion, modelling the GrIS evolution across the PPT using a recent method specifically designed for transient ice sheet experiments shows that the long-lasting paradigm explaining the large GrIS onset occurring at 2.7 Ma due to a minimum of summer insolation at 65°N is, by far, too simplistic. In our simulations, the GrIS volume appears very sensitive to pCO₂ changes during the PPT interval. Our experiments demonstrate that pCO₂ levels have to remain below 320 ppmv in order to develop and maintain perennial large GrIS across the whole PPT interval, although the absolute pCO₂ values may be model dependent. When forced by the most recent pCO₂ reconstruction available³³ for this interval, the simulated GrIS shows a good agreement with IRD records from adjacent locations⁷,41,43. The evolution of the GrIS across the PPT under this pCO₂ scenario shows large variations of ice volume of continent-sized amplitude

but indicates a continuous presence of ice even during intervals of higher pCO₂ and stronger summer insolation. This result supports recent work suggesting a very dynamic, yet persistent GrIS for the last millions of years⁴⁴, possibly receding to very small ice centres⁴⁵. Ultimately, our work emphasizes the crucial role of pCO₂ in shaping the evolution of the cryosphere across the PPT and provides numerical arguments for a vigorously dynamic GrIS under pCO₂ levels similar to, or lower than modern values.

## Methods

**Model description.** The climate model used in this study is the IPSL-CM5A GCM³⁰. The atmosphere component is the LMDZ5A version of the LMDz model (including the ORCHIDEE land-surface model) with a resolution of 3.75°×1.875° and 39 vertical layers. More details about the physical parameterization can be found in³⁰,⁴⁶. The ocean model is NEMOv3.2⁴⁷, which integrates the dynamical model OPA, the LIM2 sea ice model and the PISCES biogeochemical model. NEMO runs on a tri-polar grid (one pole in the Southern Hemisphere under Antarctica and two poles in the Northern Hemisphere under North America and Asia) in order to improve the representation of ocean dynamics in the northern high-latitudes. There are 31 unequally spaced vertical levels and a nominal resolution of 2° that is refined up to 0.5° in the equatorial area. The atmosphere and ocean models are linked through the coupler OASIS⁴⁸, ensuring energy and water conservation. Additional details about the IPSL-CM5A model can be found in Dufresne et al.³⁰.

The ice sheet model used in this study is the GRenoble Ice-Shelf and Land-Ice model (GRISLI). GRISLI is a three-dimensional thermo-mechanical model that
simulates the evolution of ice sheet geometry (extension and thickness) and the coupled temperature–velocity fields in response to climate forcing. A comprehensive description of the model can be found in ref. 31 and ref. 49. In this study, the GRISLI model is nested to the Greenland region on a cartesian grid of 15 km × 15 km. Over the grounded part of the ice sheet, the ice flow resulting from internal deformation is governed by the shallow-ice approximation⁵⁰. The model also deals with ice flow through ice shelves using the shallow-shelf approximation⁵¹ and predict the large-scale characteristics of ice streams using criteria based on the effective pressure and hydraulic load. At each time step, the velocity and vertical profiles of temperature in the ice are computed, as well as the new geometry of the ice sheet. The isostatic adjustment of the bedrock in response to the ice load is governed by the flow of the asthenosphere, with a characteristic time constant of 3000 years, and by the rigidity of the lithosphere. The temperature field is computed both in the ice and in the bedrock by solving a time-dependent heat equation. Here, as the ice sheet model GRISLI is not synchronously coupled with the IPSL-CM5A model, temperature and precipitation fields are asynchronously passed from IPSL-CM5A to GRISLI. The surface mass balance is defined as the sum between accumulation and ablation computed by the positive degree-day (PDD) method⁵².

**Interpolation method.** In the absence of synchronous coupling between GRISLI and IPSL-CM5A, we utilise an interpolation that asynchronously couples both models and offers the possibility of carrying long-term numerical integration of the ice sheet model while accounting for the time evolution of the main climate forcings. This method consists of building a matrix of possible climate states that are generated by IPSL-CM5A under various combinations of forcings, which are chosen from the range of possible values taken by forcings. The ice sheet model GRISLI can then be continuously forced by temperature and precipitation fields obtained by interpolation between the different IPSL-CM5A climatic states, interpolation based on the time evolution of the forcings. The principle of this method has been described in details in Pollard (2010) and an early version of it has been applied to the Eocene-Oligocene Transition (EOT) in Antarctica⁵³. In this work, we use an improved version of this method, which has first been applied to the EO glaciation²⁹ and that we have adapted to Greenland. Specifically here, we build a three-dimensional matrix to account for the three main drivers of an ice sheet evolution that are: (1) realistic insolation variations, (2) the atmospheric pCO₂ evolution and (3) the ice sheet feedbacks on itself. The matrix hence comprises temperature and precipitation fields that are obtained from reference IPSL-CM5A runs initialized with different combinations of orbital parameters, pCO₂ and ice sheet size. In a second step, the temperature and precipitation fields that force GRISLI can be computed by interpolating between the reference IPSL-CM5A climatic states, based on the current value taken by the 65°N summer insolation and by the atmospheric pCO₂ and on the instantaneous size of the ice sheet in GRISLI. This ensures that the climatic fields passed to GRISLI are appropriately updated at each time step to follow the evolution of both external (insolation and pCO₂) and internal (ice sheet geometry) forcings.

**Experiment design.** In this study, we have chosen two orbital configurations that produce respectively the maximal (warm orbit) and the minimal (cold orbit) mean summer insolation at 65°N between 3.0 and 2.5 Ma, according to the astronomical solution calculated by the model of ref. 36. At each time step of the ice sheet model simulations for the 3–2.5 Ma period, the impact of the summer insolation can be included by appropriately interpolating between reference IPSL-CM5A runs with warm or cold orbits.

The reference IPSL-CM5A runs are initialised with four different pCO₂ concentrations: 220 ppmv, 280 ppmv, 360 ppmv and 405 ppmv. Late Pliocene pCO₂ records indeed document a range of variation of the atmospheric pCO₂ comprised between ~200 ppmv and ~400 ppmv (See in Supplementary Fig. 2). Similarly, the instantaneous value of pCO₂ over the course of the simulation for the 3–2.5 Ma periods can be interpolated between reference runs with aforementioned pCO₂ values. As pCO₂ and temperature are linked via a logarithmic relationship, we prescribe a logarithmic interpolation between fixed pCO₂ reference runs. Conversely, the interpolation is kept linear for the insolation.

To obtain the reference Greenland ice sheet sizes that are prescribed in the reference IPSL-CM5A runs, we carry out a preliminary experiment in which we model the ice sheet development in an offline, one-way regrowth experiment (see ref. 24 for details). We initialise IPSL-CM5A with standard PlioMIP phase 1 conditions⁵⁴, which are then modified to start from an ice-free Greenland and extremely favourable conditions for glacial inception represented by the cold orbit described above and a pCO₂ concentration of 220 ppmv. GRISLI is then force with constant temperature and precipitation fields from the IPSL-CM5A simulation until the ice sheet reaches equilibrium. The ice sheet geometry is then fed back to IPSL-CM5A while pCO₂ and orbital parameters are kept identical. New climatic fields are thus obtained and GRISLI further simulates the ice sheet gain until a new equilibrium is reached. The geometry of the ice sheet can then be passed again to IPSL-CM5A. This process is repeated until a new iteration does not markedly increase the ice volume. Here, seven iterations allow us to obtain seven Greenland ice sheet sizes ranging from very small to nearly full size ice sheet (Supplementary Fig. 1c-i).

The matrix of reference IPSL-CM5A climatic states is then built using the forcings described above. A total of 2 (for the insolation) × 4 (for pCO₂) × 7 (for Greenland ice sheet sizes) simulations are run in parallel to provide reference T and P fields that cover the range of possible variations of the three main drivers that are insolation, pCO₂ and ice sheet size. A continuous T and P forcing can then be calculated based on the prescribed (for insolation and pCO₂) and emerging (for ice sheet) evolutions of these drivers. Although complex and fastidious to implement, one particular advantage of this method is that it allows virtually any pCO₂ scenario to be tested without additional GCM runs.

It should be noted that vegetation feedbacks linked to vegetation changes under Late Pliocene conditions are taken into account in our IPSL-CM5A simulations, since the tundra-taiga feedback has been shown to play a role in the onset of NH glaciation²⁰. We divide the IPSL-CM5A boundary conditions into three types relative to their presumed impact on vegetation: cold, intermediate and warm. The cold conditions are defined by the combination of the cold orbit and either 220 ppmv or 280 ppmv of pCO₂ concentration. In the reference IPSL-CM5A simulations whose boundary conditions fall under the cold criterion, we modify the PlioMIP vegetation map⁵⁵ by specifying tundra north of 50°N. The intermediate conditions are defined by the combination of the cold orbit and either 360 ppmv or 405 ppmv of pCO₂ concentration and the PlioMIP vegetation map is modified by specifying tundra north of 65°N. Finally, the warm conditions are defined by the warm orbit, regardless of the pCO₂ concentration. Under these conditions, we keep the PlioMIP vegetation map unchanged. All the reference AOGCM experiments are summarized in Supplementary Table. 1.

**Model evaluation.** We assess our model performances and our modelling strategy by presenting first a comparison between the simulated surface mass balance (SMB) obtained from the CMIP5 historical IPSL-CM5A experiment³⁰ and averaged over the period 1981–2005 and the SMB computed with the state-of-the-art regional polar model MAR⁵⁶ forced by reanalyses over the same interval. Second, we present a sensitivity experiment of the Eemian deglaciation (150–110 Ka) in which the pCO₂ evolution is very well constrained.

Supplementary Fig. 4 shows that the simulated SMB with our climate model IPSL-CM5A is consistent with the reconstructed SMB from MAR, except in some regions of the northwest and southeast Greenland. These differences can be attributed: (1) to the coarser resolution of IPSL-CM5A compared to MAR, which provide some limitations in the model's ability to accurately simulate precipitation in mountainous regions (a common problem in fully coupled climate models, see e.g. ref. 57), such as South Greenland here; (2) to biases in simulated temperatures and precipitations in the Greenland region in IPSL-CM5A³⁰ and/or MAR. Yet overall, the simulated SMB using IPSL-CM5A is reasonably close to that of MAR.

In order to simulate the Eemian deglaciation, we have performed four additional climate simulations with a LGM ice sheet configuration, two low pCO₂ scenarios (185 ppmv and 220 ppmv) and our two extreme orbital configurations. We also restrict the ice sheet simulation to Greenland. Because the climate states added in the matrix have extensive ice sheet and sea ice covers over the Northern Hemisphere, we allow the ice sheet model to use variable basal melting factors following the multi-proxy index of Quiquet et al.⁵⁸. Basal melting factors thus vary linearly between 0 and 5 m/year in regions where depth are lower than 1300 m and are fixed to 10 m/year where depth greater than 1300 m. In addition, we use the sea level curve of Waelbroeck et al.⁵⁹ as a forcing to determine which ice model grid points are floating or grounded. These two processes are added to take into account the climatic influence of the LIG deglaciation outside Greenland because we do not explicitly model this ice sheet evolution. Finally, the ice sheet model is initialized with the LGM GrIS.

The results from this sensitivity experiment are shown on Supplementary Fig. 5. The combined increase in summer insolation and pCO₂ leads to the waning of the GrIS between 133 and 123 Ka. The timing of the deglaciation is consistent with results from other studies (e.g., ref. 60–62). Stone et al.⁶⁰ “best score” simulation places the start of the GrIS deglaciation at 135 Ka (their Fig. 6b, black lines) and the end around 123 Ka. Goelzer et al.⁶² restrict their figures to the interval 130–115 Ka so that the start of the deglaciation in their simulations is unclear but the lowest GrIS volume is also found around 123 Ka. Finally, Bradley et al.⁶¹ also place the beginning and end of the GrIS deglaciation at 133 and 123 Ka, respectively (their Fig. 2). The contribution of the GrIS to sea level rise during the LIG is poorly constrained and estimated between 0.6 and 3.5 m⁶¹. Our simulated minimal GrIS reaches a volume of 2.846 × 10¹⁵ km³ (Supplementary Fig. 6), which is slightly less than the equilibrium GrIS volume (2.87 × 10¹⁵ k m³) simulated by GRISLI under preindustrial IPSL-CM5A forcings³⁰. Compared to other studies, our simulated LIG GrIS does not reach a volume low enough. However, considering that we just adapted a method that is designed for the PPT and not specifically for the Eemian, we argue that the model does a reasonable job in reproducing the Eemian deglaciation. In summary, the reasonable agreement between our modelled deglaciation and the results from other transient LIG simulations provides confidence in our modelling strategy.

**Code availability.** The codes of the IPSLCM5 and GRISLI models are available on request to the authors.
Data availability

Data that support the results of this study are available on request to the authors.

Received: 17 January 2018 Accepted: 10 October 2018

Published online: 12 November 2018

References

1. Mudelsee, M. & Raymo, M. E. Slow dynamics of the Northern Hemisphere glaciation. *Paleoceanography* **20**, 4022–4036 (2005).

2. Lawrence, K. T., Herbert, T. D., Brown, C. M., Raymo, M. E. & Haywood, A. M. High-amplitude variations in north atlantic sea surface temperature during the early Pliocene warm period. *Paleoceanography* **24**, 1–15 (2009).

3. Herbert, T. D., Peterson, L. C., Lawrence, K. T. & Liu, Z. Tropical ocean temperatures over the past 3.5 million years. *Science* **328**, 1530–1534 (2010).

4. Venti, N. L., Billups, K. & Herbert, T. D. Increased sensitivity of the Plio-Pleistocene northwest Pacific to obliquity forcing. *Earth Planet. Sci. Lett.* **384**, 121–131 (2013).

5. Lisiecki, L. E. & Raymo, M. E. A Pliocene-Pleistocene stack of 57 globally distributed benthic δ¹⁸O records. *Paleoceanography* **20**, 1003–1020 (2005).

6. Brigham-Grette, J. et al. Pliocene warmth, polar amplification, and stepped Pleistocene cooling recorded in NE Arctic Russia. *Science* **340**, 1421–1427 (2013).

7. Jansen, E., Fronval, T., Rack, F. & Channell, J. E. T. Pliocene-Pleistocene ice rafting history and cyclicity in the Nordic Seas during the last 3.5 Myr. *Paleoceanography* **15**, 709–721 (2000).

8. Maslin, Ma, Li, X. S., Loutre, M. F. & Berger, A. The contribution of orbital forcing to the progressive intensification of Northern Hemisphere glaciation. *Quat. Sci. Rev.* **17**, 411–426 (1998).

9. Flesche Kleiven, H., Jansen, E., Fronval, T. & Smith, T. M. Intensification of Northern Hemisphere glaciations in the circum Atlantic region (3.5-2.4 Ma) - Ice-rafted detritus evidence. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **184**, 213–223 (2002).

10. Eiriksson, J. & Geirsdóttir, Á. A record of Pliocene and Pleistocene glaciations and climatic changes in the North Atlantic based on variations in volcanic and sedimentary facies in Iceland. *Mar. Geol.* **101**, 147–159 (1991).

11. GEIRSDÓTTIR, Á. & EIRIKSSON, J. O. N. Sedimentary facies and environmental history of the Late-glacial glaciomarine Fossvogur sediments in Reykjavik, Iceland. *Boreas* **23**, 164–176 (1994).

12. Ramstein, G. Climates of the Earth and cryosphere evolution. *Surv. Geophys.* **32**, 329 (2011).

13. Thiede, J. et al. Million years of Greenland Ice Sheet history recorded in ocean sediments. *Polarforschung* **80**, 141–149 (2011).

14. Jansen, E. & Sjöholm, J. Reconstruction of glaciation over the past 6 Myr from ice-borne deposits in the Norwegian Sea. *Nature* **349**, 600–603 (1991).

15. Larsen, H. C. et al. Seven million years of glaciation in Greenland. *Science* **264**, 952–955 (1994).

16. De Schepper, S. et al. Northern hemisphere glaciation during the globally warm early Late Pliocene. *PLoS ONE* **8**, e81508 (2013).

17. Dolan, A. M. et al. Modelling the enigmatic Late Pliocene Glacial Event - Marine Isotope Stage M2. *Glob. Planet. Change* **128**, 47–60 (2015).

18. Tan, N. et al. Exploring the MIS M2 glaciation occurring during a warm and high atmospheric CO₂Pliocene background climate. *Earth Planet. Sci. Lett.* **472**, 266–276 (2017).

19. Haywood, A. M. et al. On the identification of a Pliocene time slice for data-model comparison. *Philos. Trans. A Math. Phys. Eng. Sci.* **371**, 1–21 (2013).

20. Koenig, S. J., DeConto, R. M. & Pollard, D. Late Pliocene to Pleistocene sensitivity of the Greenland Ice Sheet in response to external forcing and internal feedbacks. *Clim. Dyn.* **37**, 1247–1268 (2011).

21. Hill, D. J., Haywood, A. M., Hindmarsh, R. C. A. & Valdes, P. J. in *Deep-Time Perspectives on Climate Change: Marrying the Signal From Computer Models and Biological Proxysex* (eds Williams, M. E et al.) 517-538 (Geological Society of London, 2007).

22. Koenig, S. J. et al. Ice sheet model dependency of the simulated Greenland Ice Sheet in the mid-Pliocene. *Clim. Past.* **11**, 369–381 (2015).

23. Lunt, D. J., Foster, G. L., Haywood, A. M. & Stone, E. J. Late Pliocene Greenland glaciation controlled by a decline in atmospheric CO₂ levels. *Nature* **454**, 1102–1105 (2008).

24. Contoux, C., Dumas, C., Ramstein, G., Jost, A. & Dolan, A. M. Modelling Greenland ice sheet inception and sustainability during the Late Pliocene. *Earth Planet. Sci. Lett.* **424**, 295–305 (2015).

25. Berger, A., Li, X. S. & Loutre, M.-F. Modelling northern hemisphere ice volume over the last 3 Ma. *Quat. Sci. Rev.* **18**, 1–11 (1999).

26. De Boer, B., de Wal, R. S. W., Bintanja, R., Lourens, L. J. & Tuenter, E. Cenozoic global ice-volume and temperature simulations with 1-D ice-sheet models forced by benthic δ¹⁸O records. *Ann. Glaciol.* **51**, 23–33 (2010).

27. Stap, L. B., van de Wal, R. S. W., de Boer, B., Bintanja, R. & 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).

28. Willeit, M., Ganopolski, A., Calov, R., Robinson, A. & Maslin, M. The role of CO₂ decline for the onset of Northern Hemisphere glaciation. *Quat. Sci. Rev.* **119**, 22–34 (2015).

29. Ladant, J.-B., Donnadieu, Y., Lefebvre, V. & Dumas, C. The respective role of atmospheric carbon dioxide and orbital parameters on ice sheet evolution at the Eocene-Oligocene transition. *Paleoceanography* **29**, 810–823 (2014).

30. Dufresne, J.-L. et al. Climate change projections using the IPSL-CM5 Earth System Model: from CMIP3 to CMIP5. *Climate Dynamics* **40**, 2123–2165 (2013).

31. Ritz, C., Rommelaere, V. & Dumas, C. Modeling the evolution of Antarctic ice sheet over the last 420,000 years: Implications for altitude changes in the Vostok region. *J. Geophys. Res.* **106**, 31943 (2001).

32. Bartoli, G., Hönisch, B. & Zeebe, R. E. Atmospheric CO₂ decline during the Pliocene intensification of Northern Hemisphere glaciations. *Paleoceanography* **26**, 4213–4227 (2011).

33. Martínez-Botí, M. A. et al. Plio-Pleistocene climate sensitivity evaluated using high-resolution CO₂ records. *Nature* **518**, 49–54 (2015).

34. Seki, O. et al. Alkenone and boron-based Pliocene pCO₂ records. *Earth Planet. Sci. Lett.* **292**, 201–211 (2010).

35. Van De Wal, R. S. W., De Boer, B., Lourens, L. J., Köhler, P. & Bintanja, R. Reconstruction of a continuous high-resolution CO₂ record over the past 20 million years. *Clim. Past.* **7**, 1459–1469 (2011).

36. Laskar, J. et al. A long-term numerical solution for the insolation quantities of the Earth. *Astron. Astrophys.* **428**, 261–285 (2004).

37. DeConto, R. M. et al. Thresholds for Cenozoic bipolar glaciation. *Nature* **455**, 652–656 (2008).

38. Alley, R. B. et al. History of the Greenland Ice Sheet: paleoclimatic insights. *Quat. Sci. Rev.* **29**, 1728–1756 (2010).

39. Dolan, A. M. et al. Using results from the PlioMIP ensemble to investigate the Greenland Ice Sheet during the mid-Pliocene Warm Period. *Clim. Past.* **11**, 403–424 (2015).

40. Knies, J. et al. The Plio-Pleistocene glaciation of the Barents Sea - Svalbard region: a new model based on revised chronostratigraphy. *Quat. Sci. Rev.* **28**, 812–829(2009).

41. Bailey, I. et al. An alternative suggestion for the Pliocene onset of major northern hemisphere glaciation based on the geochemical provenance of North Atlantic Ocean ice-rafted debris. *Quat. Sci. Rev.* **75**, 181–194 (2013).

42. Andrews, J. T. Icebergs and iceberg rafted detritus (IRD) in the North Atlantic: facts and assumptions. *Oceanography* **13**, 100–108 (2000).

43. Sarntlein, M. et al. Mid-Pliocene shifts in ocean overturning circulation and the onset of Quaternary-style climates. *Clim. Past.* **5**, 269–283 (2009).

44. Bierman, P. R., Shakun, J. D., Corbett, L. B., Zimmerman, S. R. & Rood, D. H. A persistent and dynamic East Greenland Ice Sheet over the past 7.5 million years. *Nature* **540**, 256 (2016).

45. Schaefer, J. M. et al. Greenland was nearly ice-free for extended periods during the Pleistocene. *Nature* **540**, 252 (2016).

46. Hourdin, F. et al. The LMDZ4 general circulation model: climate performance and sensitivity to parametrized physics with emphasis on tropical convection. *Clim. Dyn.* **27**, 787–813 (2006).

47. Madec, G. & Imbard, M. A global ocean mesh to overcome the North Pole singularity. *Clim. Dyn.* **12**, 381–388 (1996).

48. Valcke, S. PRISM: An Infrastructure Project for Climate Research in Europe. Rep. No 3 (CERFACS, 2006).

49. Peyaud, V., Ritz, C. & Krinner, G. Modelling the Early Weichselian Eurasian Ice Sheets: role of ice shelves and influence of ice-dammed lakes. *Clim. Past.* **3**, 375–386 (2007).

50. Morland, L. W. Thermomechanical balances of ice sheet flows. *Geophys. Astrophys. Fluid Dyn.* **29**, 237–266 (1984).

51. MacAyeal, D. R. Large-scale ice flow over a viscous basal sediment: theory and application to ice stream B, Antarctica. *J. Geophys. Res.* **94**, 4071 (1989).

52. Fausto, R. S., Ahlström, A. P., Van As, D., Bøggild, C. E. & Johnsen, S. J. A new present-day temperature parameterization for Greenland. *J. Glaciol.* **55**, 95–105 (2009).

53. DeConto, R. M. & Pollard, D. Rapid Cenozoic glaciation of Antarctica induced by declining atmospheric CO₂. *Nature* **421**, 245–249 (2003).

54. Contoux, C., Ramstein, G. & Jost, A. Modelling the mid-pliocene warm period climate with the IPSL coupled model and its atmospheric component LMDZ5A. *Geosci. Model Dev.* **5**, 903–917 (2012).

55. Salzmann, U., Haywood, A. M., Lunt, D. J., Valdes, P. J. & Hill, D. J. A new global biome reconstruction and data-model comparison for the Middle Pliocene. *Glob. Ecol. Biogeogr.* **17**, 432–447 (2008).
56. Fettweis, X. et al. Important role of the mid-tropospheric atmospheric circulation in the recent surface melt increase over the Greenland ice sheet. *Cryosph* 7, 241–248 (2013).

57. Sun, Y., Solomon, S., Dai, A. & Portmann, R. W. How often does it rain? *J. Clim.* **19**, 916–934 (2006).

58. Quiquet, A., Ritz, C., Punge, H. J. & Y Mélia, D. Greenland ice sheet contribution to sea level rise during the last interglacial period: a modelling study driven and constrained by ice core data. *Clim. Past.* **9**, 353–366 (2013).

59. Waelbroeck, C. et al. Sea-level and deep water temperature changes derived from benthic foraminifera isotopic records. *Quat. Sci. Rev.* **21**, 295–305 (2002).

60. Stone, E. J., Lunt, D. J., Annan, J. D. & Hargreaves, J. C. Quantification of the Greenland ice sheet contribution to Last Interglacial sea level rise. *Clim. Past.* **9**, 621–639 (2013).

61. Bradley, S. L., Reerink, T. J., van de Wal, R. S. W. & Helsen, M. M. Simulation of the Greenland Ice Sheet over two glacial-interglacial cycles: investigating a sub-iceshelf melt parameterization and relative sea level forcing in an ice-sheet-ice-shelf model. *Clim. Past* **14**, 619–635 (2018).

62. Goelzer, H., Huybrechts, P., Loutre, M.-F. & Fichefet, T. Impact of ice sheet meltwater fluxes on the climate evolution at the onset of the Last Interglacial. *Clim. Past.* **12**, 1721–1737 (2016).

63. Naafs, B. D. A. et al. Late Pliocene changes in the North Atlantic Current. *Earth Planet. Sci. Lett.* **298**, 434–442 (2010).

## Acknowledgements

We thank Camille Contoux for providing the initial ice-free Greenland bedrock configuration and Zhongshi Zhang, Alan Haywood and Florence Colleoni for discussions. We also thank Matteo Willeit for providing their model results, Masa Kageyama and Aurélien Quiquet for their help with the set-up of the Eemian sensitivity simulation and Mary Minnock for her help on improving the manuscript. This study was performed using HPC resources from GENCI-TGCC (Grant 2016-GENCI t2016012212) and supported by the French project LEFE "ComPreNdrE" (A2016-992936), French State Program Investissements d'Avenir (managed by ANR), ANR HADOC project, grant ANR-17-CE31-0010 of the French National Research Agency and the Norwegian Project "OCCP" (NFR project number 221712).

## Author contributions

N.T. carried out the modelling experiments and prepared for the first manuscript, J.-B.L. and C.D. contributed equally to the experiment design and the analysis of the results. G.R. and J.-B.L. helped to improve the paper. P.B. and E.J. helped on the data-model comparison work.

## Additional information

**Supplementary Information** accompanies this paper at https://doi.org/10.1038/s41467-018-07206-w.

**Competing interests:** The authors declare no competing interests.

**Reprints and permission** information is available online at http://npg.nature.com/reprintsandpermissions/

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_7.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2018