The penultimate deglaciation: protocol for PMIP4 transient numerical simulations between 140 and 127 ka, version 1.0

Abstract. The penultimate deglaciation (PDG, ~ 138–128 thousand years before present, hereafter ka) is the transition from the penultimate glacial maximum to the Last Interglacial (LIG, ~ 129–116 ka). The LIG stands out as one of the warmest interglacials of the last 800 ka, with high-latitude temperature warmer than today and global sea level likely higher by at least 6 meters. Considering the transient nature of the Earth system, the LIG climate and ice-sheets evolution were certainly influenced by the changes occurring during the penultimate deglaciation. It is thus important to investigate, with coupled Atmosphere-Ocean General Circulation Models (AOGCMs), the climate and environmental response to the large changes in boundary conditions (i.e. orbital configuration, atmospheric greenhouse gas concentrations, ice-sheet geometry, and associated meltwater fluxes) occurring during the penultimate deglaciation. A deglaciation working group has recently been set up as part of the Paleoclimate Modelling Intercomparison Project (PMIP) phase 4, with a protocol to perform transient simulations of the last deglaciation (19–11 ka; although the protocol covers 26–0 ka). Similar to the last deglaciation, the disintegration of continental ice-sheets during the penultimate deglaciation led to significant changes in the oceanic circulation during Heinrich Stadial 11 (~ 136–129 ka). However, the two deglaciations bear significant differences in magnitude and temporal evolution of climate and environmental changes. Here, as part of the PAGES-PMIP working group on Quaternary Interglacials, we propose a protocol to perform transient simulations of the penultimate deglaciation under the auspices of PMIP4. This design includes time-varying changes in orbital forcing, greenhouse gas concentrations, continental ice-sheets as well as freshwater input from the disintegration of continental ice-sheets. This experiment is designed for AOGCMs to assess the coupled response of the climate system to all forcings. Additional sensitivity experiments are proposed to evaluate the response to each forcing. Finally, a selection of paleo records representing different parts of the climate system is presented, providing an appropriate benchmark for upcoming model-data comparisons across the penultimate deglaciation.


better understanding of the PDG could thus improve our knowledge of the processes that led to continental ice-mass loss during the LIG.
Recent work (Masson-Delmotte et al., 2011;Landais et al., 2013;Govin et al., 2015) also depicted a sequence of events over the PDG that contrasts with the one across the last deglaciation. Paleo-proxy records indicate that the disintegration of NH ice-sheets induced a ∼80 m sea-level rise (Grant et al., 2014) between 135 and 129 ka. This is also concomitant with Heinrich 5 Stadial 11 (HS11) (e.g. Heinrich, 1988;Oppo et al., 2006;Skinner and Shackleton, 2006;Govin et al., 2015). HS11 was characterised by weak NADW formation (Oppo et al., 2006;Böhm et al., 2015;Deaney et al., 2017), cold and dry conditions in the North Atlantic region (Drysdale et al., 2009;Martrat et al., 2014;Marino et al., 2015), and gradually warmer conditions over Antarctica (Jouzel et al., 2007), associated with a sustained atmospheric CO 2 increase of ∼60 ppm between ∼134 and 129 ka (Landais et al., 2013) (Figs. 1, 2). Increasing evidence of sub-millennial scale climate changes at high and low latitudes during 10 HS11 (e.g. Martrat et al., 2014) prompts the need to refine the sequence of events across the PDG. However, this is challenging as (i) climatic reconstructions over the PDG are still scarce and most records have insufficient resolution to allow identification of centennial-to millennial-scale climatic variability; and (ii) it is difficult to establish robust absolute and relative chronologies for most paleoclimatic records across this time interval (Govin et al., 2015).
While our knowledge of the processes and feedbacks occurring during deglaciations has significantly improved over the last 15 two decades (e.g. Cheng et al., 2009;Shakun et al., 2012;Abe-Ouchi et al., 2013;Landais et al., 2013;Cheng et al., 2016), many unknowns remain. For example, our understanding of the precise roles of atmospheric and oceanic processes in leading to the waning of glacial continental ice-sheets during deglaciations is still incomplete. It is also crucial to comprehend the subsequent impacts of continental ice-sheets disintegration on the oceanic circulation and thus the climate, the terrestrial vegetation and the carbon-cycle system. 20 Numerical simulations performed with climate models provide a dynamical framework to understand the response of the Earth system to external (i.e. insolation) and internal (e.g., albedo, GHGs) forcings that culminate in deglaciations. Atmospheric and oceanic teleconnections associated with millennial-scale variability can also be studied in detail. Model-paleoclimate proxy comparisons, including snapshot experiments at 130 ka and at 126 ka, suggest that the inclusion of freshwater forcing in the North Atlantic due to the melting of NH ice-sheets could explain the relatively cold conditions in the North Atlantic and warm 25 conditions in the Southern Ocean observed in paleodata during these time periods (Govin et al., 2012;Stone et al., 2016).
However, snapshot experiments assume that the climate state is in near-equilibrium, and because relatively rapid and large changes in both internal and external forcings occur during deglaciations, transient simulations (i.e. numerical simulations with time-varying boundary conditions) are needed. These simulations also allow a more robust paleodata-modelling comparison, thus enabling the refinement of the sequence of events. 30 Transient numerical simulations have already been performed for the last deglaciation (Liu et al., 2009;Menviel et al., 2011;Roche et al., 2011;Gregoire et al., 2012;He et al., 2013;Otto-Bliesner et al., 2014) and provide a dynamical framework to further our understanding of the climate-change drivers, teleconnections and feedbacks inherent in the Earth system. Transient simulations covering the period 135 to 115 ka have also been performed with a range of models to understand the impact of surface boundary conditions and freshwater fluxes on the LIG (Bakker et al., 2013;Loutre et al., 2014;Goelzer et al., 2016a). 35 3 Greenhouse gases GHG records are available solely from Antarctic ice cores across the time interval 140-127 ka (Fig. 2). LIG GHG records from the NEEM and other Greenland ice cores are affected by stratigraphic disturbances and in-situ CO 2 , CH 4 and N 2 O production (e.g. Tschumi and Stauffer, 2000;NEEM community members, 2013). The NGRIP ice core provides a continuous and reliable CH 4 record but it only extends back to ∼123 ka (North Greenland Ice Core project members, 2004). After a brief 5 description of existing atmospheric CO 2 , CH 4 and N 2 O records (below), we recommend using recent spline-smoothed GHG curves calculated from a selection of those records (Köhler et al., 2017). They have the benefit to provide continuous GHG records, with a temporal resolution of 1 yr on the commonly-used AICC2012 gas age scale (Bazin et al., 2013;Veres et al., 2013). This time scale is associated with an average 1σ absolute error of ∼2 ka between 140 and 127 ka.
Atmospheric CO 2 concentrations have been measured on the EDC and TALDICE ice cores (Fig. 2). The EDC records from 10 Lourantou et al. (2010) and Schneider et al. (2013) agree well overall. The Schneider et al. (2013) dataset depicts a long-term CO 2 increase starting at ∼137.8 ka and ending at ∼128.5 ka with a centennial-scale CO 2 rise above the subsequent LIG CO 2 values, also referred to as an "overshoot". The CO 2 overshoot is smaller in the Schneider et al. (2013) dataset compared to a similar feature measured in Lourantou et al. (2010): while the former displays a relatively constant CO 2 concentration of ∼275 ppm between 128 and 126 ka, the latter shows a CO 2 decrease from 280 to 265 ppm between 128 and 126 ka. The 15 offsets between CO 2 records from the same EDC core are likely related to the different air extraction techniques used in the two studies . The smoothed spline CO 2 curve we recommend using as forcing for the PDG is based on those two EDC dataset and the calculation method accounts for such potential difference in local maxima (details provided in Köhler et al. (2017)).
Atmospheric CH 4 concentration records from Vostok, EDML, EDC and TALDICE agree well within the gas-age uncer-20 tainties attached to each core (Fig. 2). They illustrate a slow rise from ∼390 to 540 ppb between ∼137 ka and 129 ka that is followed by an abrupt increase of ∼200 ppb reaching maximum LIG values at ∼128.5 ka. Because CH 4 sources are located mostly in the NH, an interpolar concentration difference (IPD) between Greenland and Antarctic CH 4 records exists.
For instance, an IPD of ∼14 ppb, ∼34 ppb and ∼43 ppb is reported during the LGM, Heinrich Stadial 1 and the Bølling warming respectively (Baumgartner et al., 2012). However, without reliable CH 4 records from Greenland ice cores, it remains 25 challenging to estimate the evolution of the IPD across the deglaciation. Hence, for the atmospheric CH 4 forcing of transient PDG simulations, we recommend using the smoothed spline CH 4 curve which is solely based on the EDC CH 4 record (Köhler et al., 2017), recognising that the values may be 1-4% lower than the actual global average. Atmospheric TALDICE, EDML and EDC N 2 O records are available between 134.5 and 127 ka (Fig. 2) (Schilt et al., 2010;Flückiger et al., 2002). From 134.5 to 128 ka, N 2 O levels increase from ∼220 to 270 ppb. Following a short decrease until ∼127 ka, N 2 O concentrations stabilise afterwards. No reliable atmospheric N 2 O concentrations are available beyond 134 ka as N 2 O concentrations measured in the air trapped in ice from the penultimate glacial maximum are affected by in-situ production related to microbial activity (Schilt et al., 2010). During the LGM (considered here as the time interval 26-21 ka), the average 5 N 2 O level was ∼201 ppb. Assuming the LGM is an analogue for the penultimate glacial maximum, we propose a 140 ka spin-up value and N 2 O transient forcing curve that starts with a 201 ppb level and then linearly increases to 218.74 ppb at 134.5 ka. From 134.5 ka, we recommend using the N 2 O smoothed spline curve calculated by Köhler et al. (2017) and which is based on the TALDICE and EDC discrete N 2 O measurements.
The CO 2 and N 2 O levels from the spline curves at 127 ka (274 ppm and 257 ppb) only differ from the values chosen 10 as boundary conditions for the PMIP4 lig127k equilibrium experiment by 1 ppm and 2 ppb respectively (Otto-Bliesner et al., 2017;Köhler et al., 2017). The comparison is less direct for CH 4 . Indeed a global CH 4 value (685 ppm) rather than an Antarctic ice core-based CH 4 value (e.g. CH 4 level of 660 ppm at 127 ka in Köhler et al. (2017)) is proposed as forcing for the lig127k simulations. However, this difference in global atmospheric CH 4 and Antarctic ice core CH 4 concentration is similar to the one observed during the mid-Holocene (23 ppb) (Otto-Bliesner et al., 2017;Köhler et al., 2017). 15

Continental ice-sheets
Changes in continental ice-sheets during the deglaciation will significantly impact the climate system through their albedo, which will directly affect the radiative balance (e.g. He et al., 2013). Changes in continental ice-sheets geometry can also significantly impact atmospheric dynamics (e.g. Zhang et al., 2014;Gong et al., 2015). Transient simulations of the PDG will thus need to be forced by the 3-dimensional and time-varying evolution of continental ice-sheets, that is currently only available 20 from numerical simulations. However, simulating the evolution of continental ice-sheets across the PDG is associated with large uncertainties, due to the climate forcing of the ice-sheet models and poorly constrained non-linearities within the ice-sheet system. Glacial geological data (e.g. glacial deposits, glacial striations...) are also available to constrain continental ice-sheet evolutions and can thus provide an estimate of the uncertainties associated with the numerical ice-sheet evolutions. In this section, we describe the available numerical ice-sheet evolutions to use as a forcing of the transient simulations of the PDG. 25 We further compare the results of these simulations with existing glacial geological constraints.

Combined ice-sheet forcing
To facilitate the transient simulations of the PDG, we are providing a combined ice-sheet forcing (available on the PMIP4 wiki), in which separate reconstructions of different ice-sheets have been merged. As the sea-level solver assumes an equilibrium initial condition, the simulations start at the previous interglacial. As is standard, the solver also requires present-day ice-30 sheet histories to bias correct against present-day observed topography. Thus, a full 240 ka ice-sheet history is required. The simulated NH ice-sheet evolution, described in Section 4.2 (Abe- Ouchi et al., 2013), is merged with the simulated Greenland and the Antarctic (Briggs et al., 2014) evolutions described in Sections 4.3 and 4.4, respectively. The resolution of the merged ice-sheet file is 1 • longitude by 0.5 • latitude.
From the LIG onward, the combined ice-sheet evolution, referred to as GLAC-1D in PMIP4, is used (e.g. Ivanovic et al., 2016). GLAC-1D includes the Greenland and Antarctic ice-sheets components described in section 4. 3 and 4.4 (Briggs et al., 2014), the North American ice-sheet simulation described in Tarasov et al. (2012) and the Eurasian ice-sheet simulation de-5 scribed in Tarasov (2014). The ice-sheet thickness from these simulations are run through a sea-level solver using the VM5a (Peltier and Drummond, 2008) Earth rheology to extract a gravitationally self-consistent topography. The surface topography is then run through a global surface drainage solver (using the algorithm described in Tarasov and Peltier, 2006) to extract the relevant surface drainage pointer field for each time-slice. This will indicate in which ocean grid cell each terrestrial grid cell will drain into.

North American and Eurasian ice-sheets
The evolution of NH ice-sheets during the PDG is given by a numerical simulation performed with the thermo-mechanically coupled shallow-ice-sheet model IcIES (Ice-sheet model for Integrated Earth system Studies) with an original resolution of 1 • by 1 • in horizontal and 26 vertical levels (Abe-Ouchi et al., 2007) (Fig. 3). IcIES uses the shallow ice approximation and computes the evolution of grounded ice but not floating ice shelves. The sliding velocity is related to the gravitational driving 15 stress according to Payne (1999) and basal sliding only occurs when the basal ice is at the pressure melting point. This icesheet model was driven by climatic changes obtained from the MIROC GCM (Abe-Ouchi et al., 2013), which was forced by changes in insolation and atmospheric CO 2 concentration. In global agreement with glacial geological constraints (Dyke et al., 2002;Svendsen et al., 2004;Curry et al., 2011;Syverson and Colgan, 2011) and other numerical simulations of NH ice-sheets evolution (Tarasov et al., 2012;Abe-Ouchi et al., 2013;Peltier et al., 2015;Colleoni et al., 2016), the simulated extent and 20 volume of the North American ice-sheet was smaller during MIS 6 than MIS 2 (Fig. 3).
In Eurasia, MIS 6 recorded the most extensive glaciation since MIS 12 (Hughes and Gibbard, 2018). The maximum extent of the Fennoscandian ice-sheet probably occurred at ∼160 ka, when it extended into central Netherlands, Germany, and the Russian Plain (Margari et al., 2010;Ehlers et al., 2011;Hughes and Gibbard, 2018). This was followed by a partial melting of the Fennoscandian ice-sheet, peaking between ∼157-154 ka, and a readvance after 150 ka (Margari et al., 2010;25 Hughes and Gibbard, 2018). The maximum extent of the NH ice-sheets probably occurred at the end of MIS 6 (Margari et al., 2014;Head and Gibbard, 2015), In Europe the late MIS 6 ice advance was less expensive than at ∼160 ka, but that was compensated by ice-sheet expansion in Russia, Siberia (Astakhov et al., 2016) andin North America (e.g. Curry et al., 2011;Syverson and Colgan, 2011). Glacial geological constraints (e.g. Astakhov, 2004;Svendsen et al., 2004) indeed suggest that the Barents-Kara ice-sheet extended further during MIS 6 than in MIS 2. The simulated Eurasian ice-sheet is in general agree-30 ment with the reconstruction of Lambeck et al. (2006), with a dome reaching 3000 m over the Kara Sea during MIS 6 that subsequently disintegrated across the deglaciation. However, the extent and volume of the simulated Eurasian ice-sheet might be underestimated since it is smaller at MIS 6 than MIS 2, whereas reconstructions suggest it should be larger at MIS 6 (Lambeck et al., 2006;Rohling et al., 2017). Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2019-41 Manuscript under review for journal Geosci. Model Dev. Discussion started: 6 March 2019 c Author(s) 2019. CC BY 4.0 License. Rohling et al. (2017) further suggest that the ice volume was almost equally distributed between Eurasia and North America at MIS 6, with a 33 to 53 m global mean sea level equivalent (sle) contribution from the Eurasian ice-sheet and 39-59 m from North America, whereas the ice-sheet simulation produces a ∼24 m sle contribution from Eurasia. Thus, the volume of the North American ice-sheet may also be overestimated.
In the ice-sheet simulation, NH ice-mass loss follows closely the boreal summer insolation and occurs mostly between ∼134 5 and 127 ka (Fig. 4), with two peaks of glacial meltwater release at ∼131 ka and 128 ka. By 132 ka, the Eurasian ice-sheet has decreased significantly and the southern and western flanks of the North American ice-sheet have disintegrated (Fig. 3).
Another significant retreat of the North American ice-sheet occurs between 132 and 128 ka, at which point it is mostly restricted to the north of the Hudson Bay. By 127 ka, the North American ice-sheet only remains over Baffin Island. 10 The Greenland model uses an updated version (GSM.G7.31.18) of the Glacial Systems Model (e.g. Tarasov et al., 2012) run at grid resolution of 0.5 • longitude by 0.25 • latitude. The model has been upgraded to hybrid shallow-ice and shallow-shelf physics, with ice dynamical core from Pollard and DeConto (2012) and includes: a 4 km deep permafrost resolving bed thermal component (Tarasov and Peltier, 2007), visco-elastic bedrock response with global ice-sheet and sea-level loading, sub-shelf melt, parametrizations for subgrid mass-balance and ice flow (Morzadec et al., 2015), and updated parametrizations for surface 15 mass-balance and ice calving.

Greenland ice-sheet
Model runs start at 240 ka with present-day ice and bedrock geometry and with an ice and bed temperature field from the end of a previous 240 ka model run. The model is then forced from 240 ka until 0 ka, with a climate forcing that is partly glacial index based, using a composite of a glaciological inversion of the GISP II regional temperature change (for the last 40 ka) and the synthetic Greenland δ 18 O curve that was deduced from the Antarctic EDC isotopic record assuming a thermal 20 bipolar seesaw pattern (Barker et al., 2011). The climate forcing also includes 2-way coupled 2D energy balance climate model (Tarasov and Peltier, 1997) to capture radiative changes.
Greenland ice-sheet model runs are scored against a large set of constraints including relative sea level (RSL), proximity to present-day ice-surface topography, present-day observed basal temperatures from various ice cores, time of deglaciation of Nares Strait, and the location of the present-day summit. The last 20 ka of the run is critical as this represents the time 25 period with most of the data constraints for Greenland. The simulation presented here (G9175) is a least misfit model from a preliminary exploratory ensemble.
This simulation suggests no significant change in Greenland ice-mass until ∼134 ka (Fig. 4), followed by a small icemass loss, mostly from floating ice, between 134 and 130 ka. In this simulation, the main phase of Greenland deglaciation occurs between 130 and 127 ka, during which Greenland loses its glacial grounded ice volume (∼2.9 m sle), but also loses 30 an additional 1.5 m sle, compared to the pre-industrial Greenland configuration. As shown in Figure 5, the extent and height of the Greenland ice-sheet is significantly smaller at 128 ka than at 132 ka. Greenland ice-mass loss is particularly evident on its western side, with a part of southwestern Greenland being ice-free. To a first order, the simulated disintegration of the Greenland ice-sheet follows the increase in boreal summer insolation and in atmospheric CO 2 (Fig. 1).
The main phase of the Greenland ice-sheet retreat in this simulation is globally in agreement with proxy records, which suggest significant runoff in the Labrador Sea at ∼130 ka and at ∼127 ka (e.g. Carlson and Winsor, 2012). However the simulated Greenland ice-sheet disintegration could be too rapid as paleoproxy records suggest significant meltwater discharge from the Greenland ice-sheet throughout the LIG (e.g. Carlson and Winsor, 2012). In addition other model simulations suggest a maximum sea-level contribution from Greenland at ∼123-121 ka (Yau et al., 2016;Bradley et al., 2018), in agreement with 5 the timing of the LIG minimum elevation at the Greenland NEEM location. This minimum elevation estimate was reconstructed from total air and water isotopic records measured on the deep ice core drilled at that site (NEEM community members, 2013).
Based on the paleoproxy records and model simulations, the protocol for the PMIP4 LIG simulation for 127ka (lig127k) recommends a pre-industrial Greenland configuration. First, SST dependence is added to the sub-shelf melt model. Second, to compensate for inadequate LIG warming, where SSTs are above present-day values, they are then given a minimum value of 3 • C (i.e. SST=MAX(SST,3.0 • C)). Even so, the Antarctic contribution to the LIG high-stand is only 1.4 m sle and is therefore inadequate given current inferences (as well as constraints on contributions from Greenland and steric effects) (e.g. Kopp et al., 2009). 20 The simulation suggests a continuous Antarctic ice-sheet discharge during the PDG, with a glacial ice-mass loss of ∼12.5 m sle between 140 and 131 ka, followed by an additional 1.4 m sle between 131 and 130 ka (Fig. 4). In this simulation, the West Antarctic ice-sheet loses significant ice-mass between 140 and 136 ka (Fig. 6), with a retreat of the grounding line over the Ross Sea as well as ice-mass loss in the Weddell Sea, on the Antarctic Peninsula and in the Amundsen Sea sector. By 132 ka, the grounding line has completely retreated over the Ross Sea and has retreated significantly over the Weddell Sea. 25

Sea-level
Direct evidence for constraining the evolution of the global sea level during the time interval 140-127 ka remains sparse.
Although the LR04 benthic δ 18 O stack (Lisiecki and Raymo, 2005) is sometimes used to approximate sea-level change on glacial-interglacial timescales, in the case of the PDG, the timing of the LR04 benthic δ 18 O stack is fixed by reference to a handful of U-series coral dates from Huon Peninsula with relatively high analytical uncertainties and questionable preservation 30 (Bard et al., 1990;Stein et al., 1993). Tying the MIS 5e peak to the average age of these coral dates results in a benthic δ 18 O minimum that is roughly centred on the main phase of coral growth during this interglacial period (122 -129 ka) (Stirling et al., Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2019-41 Manuscript under review for journal Geosci. Model Dev. Discussion started: 6 March 2019 c Author(s) 2019. CC BY 4.0 License. 1998) rather than having the onset of the interglacial aligned with the timing of the onset of the sea level high-stand at far-field sites (∼129 ka) (e.g. Stirling et al., 1998;Dutton et al., 2015).
Here, we seek to provide an improved reconstruction of sea level across the PDG by examining available RSL records.
Information on the timing and magnitude of the changes across this time interval is provided by three RSL records (Fig. 7): i) A RSL record from the Red Sea (Grant et al., 2012), that is deduced from the planktic foraminifera δ 18 O measured on Providing a robust age model for sediment records from MIS 6 to the LIG is not straightforward (e.g. Govin et al., 2015) and over time, several age models have been proposed for the Red Sea RSL record (e.g. Siddall et al., 2003;Rohling et al., 2009;Grant et al., 2012). The latest chronology is based on climatic alignment of the Red Sea RSL record to eastern Mediterranean planktic foraminifera δ 18 O records, which are in turn aligned onto the absolutely-dated Soreq cave speleothem δ 18 O record 15 (Grant et al., 2012). While the absolute ages of the speleothem record have the potential to provide a more robust age model (both in terms of accuracy and precision), the application of these dates to the Red Sea sea level reconstruction hinges on the assumption that the tie points between the Red Sea and eastern Mediterranean records have been correctly assigned, and that the intervals between these tie points can be linearly extrapolated.
The Tahiti and Huon Peninsula corals are associated with absolute radiometric dates (using U-series geochronology). For 20 the purpose of this study, all of the U-series ages have been recalculated to normalize them with the same set of decay constants for 234 U and 230 Th (Cheng et al., 2013) (Tables S3 and S4), using the methodology described by Hibbert et al. (2016). The array of data from Huon Peninsula suggest post-depositional alteration (open-system behaviour of the U-series isotopes) that complicates a precise age interpretation (Fig. 7).
The Red Sea time series published by Grant et al. (2012) depicts that, after a RSL low stand of about -100 m relative to provide bounding ages on the timing of this sea level rise pulse, with ages of corals that grew at 135.0 ka (in 0-6 m water depth) and 133.5 (±1) ka (0-25 m water depth). In between these shallower facies, there is a deeper water facies (≥20 m paleowater depth), but there are no reliable ages within this interval of the core (Thomas et al., 2009). This observation, based on changes in both the lithofacies and benthic foraminiferal assemblage, is interpreted as a pulse of sea-level rise in between about 135.0 and 133.5 ka (Fujita et al., 2010). A similar sea level oscillation has also been interpreted based on the stratigraphy as well 5 as the age and paleowater depth reconstruction at Huon Peninsula (Esat et al., 1999). The absolute timing of coral growth is only loosely constrained at this site due to open-system behaviour of the U-series isotopes (as reflected by the scatter in ages of corals collected in Aladdin's cave, ∼134 to 126 ka, Fig. 7) (Esat et al., 1999). Indeed, the corals from Terrace VII have ages (with high uncertainty) ranging from about 137 to 134.5 ka and the corals from the cave have a wide range of ages, from 134.1 to 125.9 ka (more details in Esat et al. (1999)). Given that the younger end of this age range is clearly within the MIS 5e 10 sea-level highstand (e.g. Stirling et al., 1998), it is more likely that the older end of this diagenetic array of data from Aladdin's cave is a better approximation for a the primary age (i.e. it is closer to the unaltered end member). Despite these diagenetic concerns, the agreement in the timing of this PDG sea level reversal (MWP-2A) in Tahiti and Huon Peninsula is striking (Fig.   7b).
When considering the 95% probabilistic intervals of the Red Sea RSL reconstruction on the chronology from Grant et al. 15 (2012), an overlap is observed with the coral data over the MWP-2A interval, within the stated uncertainties. Still, both coral dataset suggest that MWP-2A occurs several millennia later (i.e. ∼135-134 ka) than in the Red Sea RSL reconstruction. This mismatch is likely to be related to the difficulty to precisely anchor the dating of the current Red Sea RSL age scale over this interval (as also discussed in the supplementary information of Grant et al. (2012)). Hence, we propose a revised chronology for the Red Sea RSL record in order to provide a better agreement with the absolutely-dated corals. Given the potential ambiguities 20 of the tie point defined in Grant et al. (2012) to stretch the depth scale across this interval, we find it reasonable to adjust it such that the timing of MWP-2A is more consistent with the absolute ages provided by the Tahiti and Huon Peninsula coral data.
We note that reassigning the tie points across this interval (Table S5), where tie points are placed at the beginning and end of MWP-2A (as defined by the coral data), results in a sea-level reconstruction that more closely approximates a linear agedepth model (Fig. 7a). This revised age model for the Red Sea RSL is adopted as our preferred reconstruction for sea-level 25 change during the PDG. This reconstruction also compresses the total duration of the sea-level rise during the entirety of the PDG transition, which has implications for the freshwater forcing in the NH and for making analogies between the last and penultimate deglaciations.
This revised chronology is still attached to large uncertainties given the limits of the datasets. Also, considerable uncertainties remain with the magnitude of the sea-level pulse during MWP-2A because some of the corals cover a wide range of paleowater Glacial isostatic adjustment to the deterioration of the MIS 6 ice-sheets will also differentially affect Tahiti and Huon Peninsula, which precludes a direct comparison of the magnitude of sea-level change between these sites or a direct interpretation of global mean sea-level change in the absence of modelling. Because the changes in global mean sea level are rapid across the penultimate deglaciation, the eustatic signal is likely dominant, leading to a timing of the rapid changes that is similar between local RSL and global mean sea-level reconstructions. Still, the rate of change may be different between sites due to 5 local differences in the magnitude of sea-level change. Based on the revised chronology for the Red Sea RSL and on the coral constraints, MWP-2A starts at ∼137 ka, while MWP-2B starts at ∼133 ka ( Fig. 7, Table 4).
Finally, far-field coral data from the Seychelles and Western Australia, that have been corrected for the glacial isostatic adjustment (e.g. Dutton et al., 2015), indicate that global mean sea level passed the position of modern sea level at about 129 ka (Fig. 7). The evolution of sea level during the LIG high-stand is still debated and may have included some meter-scale 10 sea-level oscillations, but by at least some accounts, it is thought to have risen a few meters between 129 and 122 ka (e.g. Kopp et al., 2009;Dutton et al., 2015). So while the timing of peak sea level may have occurred later in the interglacial (∼122 ka), the onset of the highstand (∼129 ka) could represent an inflexion point in the rate of sea-level change coming out of the rapid deglaciation and into the interglacial. Overall, eustatic sea-level reconstructions based on paleodata and continental ice-sheets simulations (Section 4) are consis-15 tent. However, the amplitude of the eustatic sea level change across the PDG estimated from the Red Sea reconstructions is ∼10 m smaller than the combined ice-sheets simulations (Fig. 4e). In addition, the sea-level data suggest a small sea-level increase at ∼140 ka, which is not present in the ice-sheet simulations. Both suggest that the main phase of sea-level rise/continental ice-sheets disintegration initiates at ∼134 ka, even if the overall magnitude is larger in the ice-sheet simulation, but with a lower rate of change than in the Red Sea reconstructions. However, it is worth keeping in mind that both the sea level data and  If a LGM run is already available, then it is suggested to initialise the 140 ka spin-up from the climate fields and ocean state produced by the LGM equilibrium run. Starting from a LGM state may minimize the duration of the spin-up, as it should shorten the time to reach equilibrium in the deep ocean (Zhang et al., 2013). If starting from a pre-industrial set up, it is suggested to follow the PMIP4 LGM protocol (Kageyama et al., 2017) to set up the 140 ka state, but using the 140 ka boundary The model should be forced with 140 ka background conditions (Table 1), including appropriate orbital parameters, GHG concentrations as averaged over the interval 141-139 ka (191 ppm CO 2 , 385 ppb CH 4 , and 201 ppb N 2 O), as well as the NH and Antarctic ice-sheets' extent, topography and associated albedo (as described above). The forcing file describing the evolution of NH and Antarctic ice-sheets, as simulated by ice-sheets models (Abe-Ouchi et al., 2013;Briggs et al., 2014) (Section 4), is available in the Design and Data sections of the PMIP4 wiki. They include the evolution of the ice-mask, as well as surface 5 and bedrock elevations. Kageyama et al. (2017) provides guidelines for computing land-ice fraction and orography from the ice-sheet reconstruction datasets. The details of these forcings and the approach taken to compute them will ultimately depend on each model resolution.
The large glacial ice-sheets of MIS 6 impacted sea level and the land-sea mask. It would be best to modify the land-sea mask resulting from ice-sheet changes. Depending on the resolution of the model, this might not be a crucial parameter, except 10 for some bathymetry and land-sea mask features, which have particular importance for ocean circulation. For example, the Bering Strait, which is 40 to 50 m deep, exerts a significant control on NADW and North Pacific Intermediate Water formation (e.g Okazaki et al., 2010;Hu et al., 2012). We therefore recommend to close the Bering Strait during the 140 ka spin-up.
Following the recommendations for the PMIP4 LGM equilibrium run (Kageyama et al., 2017), the land-sea mask should include the exposure of the Sahul and Sunda shelves in the Indo-Australian region, as well as closure of the strait between the The model should be spun-up until near equilibrium is reached. Previous PMIP protocols recommend that the simulations are considered at equilibrium when the trend in globally averaged SST is less than 0.05 • C per century and the Atlantic Meridional Overturning Circulation (AMOC) is stable. Marzocchi and Jansen (2017) recently pointed out that the AMOC should be monitored on a centennial timescale to properly assess equilibrium. Zhang et al. (2013) further suggested that the trend in 25 zonal-mean salinity in the Southern Ocean (south of the winter sea-ice edge) should remain small (less than 0.005 psu per 100 years), especially in the Atlantic sector, to avoid potential transient characteristics in the deep ocean from impacting on AMOC strength. For models including representations of the carbon cycle or dynamic vegetation, the requirement is that the carbon uptake or release by the biosphere is less than 0.01 PgC per year. Similar to the recommendation made in Kageyama et al. (2017), the outputs of at least 100 years of the equilibrated 140 ka spin-up should be made available and fully described. The main changes in boundary conditions across the deglaciation, i.e. insolation, GHG concentrations and continental icesheets, have been described in sections 2, 3 and 4, respectively and are summarised in Table 1. For all simulations, methods should be fully documented.

Orography, bathymetry, coastlines and rivers
Disintegration of continental ice-sheets during the deglaciation affected continental topography and ocean bathymetry, and thus coastal outlines and river routing. Therefore, time-varying changes in land-ice fraction, land-sea fraction, and topography should be applied. Variations in the ice mask and topography should be updated at the same time. It is up to each group to decide the appropriate time frequency at which to update this forcing. The resolution of the files provided is 500 years, but 5 higher frequency changes obtained through linear interpolations can also be performed to avoid step changes. As mentioned for the 140 ka spin-up regarding changes in land-sea fraction, particular attention should be given to the opening of the Bering Strait and the flooding of the Sunda and Sahul shelves. When possible, these should be varied across the PDG as ice-sheets disintegrate and sea-level rises. Following the combined ice-sheet history presented here, which includes some GIA adjustment, the Bering Strait might open at ∼127.5 ka. Please note that this is a first estimate, which is associated with large uncertainties. 10 As changes in the land-sea mask could impact water delivery to the ocean through rivers, it is recommended to check that river mouths are consistent with the adjusted land-sea mask (Kageyama et al., 2017). If possible and of interest, river networks could also be remapped to take into account the ice-sheet changes.

Vegetation, land surface and other forcings
The climatic and ice-sheet changes occurring between 140 and 127 ka will significantly impact the vegetation, and thus also 15 land albedo, evapo-transpiration and the terrestrial carbon pool. The preferred option would thus be to include a dynamical vegetation model, fully coupled to the atmospheric model. Care should be taken in regions where an ice-sheet is present, as the ice-sheet and its albedo should replace any possible vegetation. If a coupled dynamical vegetation model is not available, then the experiments should be run with prescribed land-surface parameters and with fixed vegetation types and plant physiology outside of the regions covered by a continental ice-sheet, as obtained from the CMIP5 pre-industrial set up (Taylor et al., 2012;20 Ivanovic et al., 2016). In that case, care should be taken to have a consistent land-surface/vegetation forcing with the adjusted land-sea mask. Vegetation/land-surface type on a newly emerged land (for example the Sunda shelf) will thus need to be fixed based on interpolations with the nearest grid points.

Freshwater forcing
Through their impact on global salinity and ocean circulation, disintegrating ice-sheets can significantly affect the climatic 25 and biogeochemical evolution of the penultimate deglaciation (Oppo et al., 1997;Cheng et al., 2009;Hodell et al., 2009;Landais et al., 2013;Deaney et al., 2017). In particular, meltwater input in the North Atlantic region could be a significant driver of changes in NADW formation (Loutre et al., 2014;Goelzer et al., 2016b;Stone et al., 2016), including the ones associated with HS11 (further discussed in Section 8.2). It is thus strongly recommended to include a carefully designed freshwater scenario when performing transient simulations of the deglaciation.

30
As much as possible, and for all scenarios, meltwater should be added in the appropriate locations to match the evolution of the ice-sheets. Freshwater can be added over an appropriate ocean area close to the disintegrating ice-sheet, or a self-consistent paleo surface drainage forcing could also be implemented. This would involve using the provided downslope routing fields to route the water flux (fwf, cm/yr) from each grid cell: with P for precipitation (cm/yr), E for evaporation (cm/yr) and (dH/dt) ice−sheet (cm/yr, assuming an ice density of 0.91 g/cm 3 ) the change in ice-sheet thickness over time as described in the ice-sheet forcing files. Modellers may need to adjust the above to 5 account for land surface model changes in water storage. In case of negative meltwater forcing, the artificial salt flux addition should be spread globally over the ocean.
Estimates of sea-level changes across the PDG suggest a global sea-level rise of ∼100 m (Grant et al., 2012(Grant et al., , 2014 (Fig. 7) between 140 and 130 ka, due to the disintegration of NH and Antarctic ice-sheets. As the Antarctic contribution to sea level is estimated at ∼13 m (Fig. 4), this leaves a ∼87 m contribution from NH ice-sheets. Estimates of meltwater input to the North Atlantic based on changes in NH ice-sheets (Fig. 4b,e,f, black, fIC) (Abe-Ouchi et al., Meltwater input estimates derived from Red Sea sea-level records (Fig. 4f, blue, fSL) on the revised chronology (section 5) display a large (up to 0.3 Sv) 1000-year-long meltwater pulse centred at 137 ka, a broad meltwater input between 134 and 20 131 ka with a large peak at 131.7 ka, and two relatively late meltwater pulses centred at 130 ka and 128.3 ka, respectively.
The magnitude and length of the AMOC perturbation resulting from such a meltwater scenario will, of course, depend on each model's sensitivity and on the initial AMOC state. However, for most models, it is expected that NADW formation would weaken significantly between ∼137.8 and 136.5 ka as well as between ∼133.5 and 129 ka. Another small AMOC perturbation is expected between ∼129 and 127.8 ka in this scenario. 25 There are significant uncertainties associated with both the simulation of NH ice-sheets and the timing and amplitude of sea-level changes. In addition, to fully explore the potential of transient deglacial simulations, it is critical to simulate NADW changes in global agreement with proxy records. Since periods of increased IRD delivery have been associated with changes in NADW (e.g. van Kreveld et al., 2000;Rodrigues et al., 2017), we design an additional meltwater scenario, for which the timing is based on North Atlantic and Norwegian Sea IRD records (Fig. 8a). To construct the fIRD scenario, the IRD record of m originating from the Eurasian ice-sheet and 51 m from the North American one. It is expected this scenario will lead to a weakening of NADW between 139.5 and 136.5 ka as well as at ∼135 ka. The sustained meltwater pulse might induce NADW cessation between ∼133.5 and 129.4 ka, followed by a recovery sometime between 129 and 128 ka.
As can be seen in Figure 4f, the meltwater forcing scenarios based on sea-level changes (fSL) and the IRD record (fIRD) share some similarities. The main differences between these scenarios are the small meltwater pulse at ∼137 ka in fSL, which 5 is of much reduced amplitude in fIRD, and the ∼128 ka pulse in fSL not present in fIRD. Finally, the fSL scenario includes periods of significant negative meltwater forcing (i.e. artificial salt flux addition), corresponding to phases of sea-level lowering.
As described above, it is suggested to use the self-consistent drainage scheme for the location of the meltwater input. For NH scenarios fSL and fIRD, (dH/dt) ice−sheet (from equation 1) would need to be scaled to obtain the appropriate meltwater flux.
For those who wish to take part in an inter-comparison involving comparable boundary conditions, the fSL scenario is put 10 forward as the recommended option (Table 1). Depending on the sensitivity of the model to freshwater input, fIRD scenario can provide a good alternative to fSL. However, the ultimate choice of the appropriate freshwater scenario is left to each group, and sensitivity experiments to assess the climatic impact of different meltwater scenarios are encouraged (section 7).
To take into account the effect of Antarctic ice-sheet melting, freshwater should also be added close to the Antarctic coast, following the self-consistent routing scheme described above (Fig. 4d) (Briggs et al., 2014). This scenario will broadly consist

Sensitivity experiments
The penultimate deglaciation is a particularly interesting period as it provides a framework to study the impact of changes in insolation, atmospheric GHG concentrations and continental ice-sheets on climate. In addition, meltwater input associated 5 with the disintegration of continental ice-sheets will also impact the oceanic circulation and thus global climate and biogeochemistry (e.g. Liu et al., 2009;Menviel et al., 2011Schmittner and Lund, 2015;Goelzer et al., 2016a;Ivanovic et al., 2017;Menviel et al., 2018). Ultimately, transient deglacial simulations will inform on the impact of each of these processes as well as their interactions.
However, the transient simulation of the PDG proposed here might present a challenge to state-of-the-art Earth-system 10 models, as they will need to include a spin-up with 140-ka boundary conditions and be run for 13,000 years. For this reason (i.e. to avoid the necessity for additional simulations when the computational expense is prohibitive), the main experiment includes all appropriate boundary forcing as well as meltwater input due to disintegrating ice-sheets. This experiment will provide valuable information on processes occurring during the PDG and will allow for a thorough comparison of the penultimate and The main experiment is described in Table 1 and includes a comprehensive meltwater scenario. Additional transient simulations of the PDG are encouraged to, for example, assess different timing and amplitude of meltwater-input (Table 2), but also simulations with globally uniform meltwater input (fUN). As there are large uncertainties associated with meltwater input scenarios and the sensitivity of deep convection to the freshwater input, it is strongly advised to perform both the main experiment 25 and fUN to isolate the impact of the freshwater input (Goelzer et al., 2016a). In addition, the response to individual forcings (i.e. orbital parameters, GHGs and ice-sheet extent and albedo) could be assessed separately (He et al., 2013;Gregoire et al., 2015).
Even though there are some geological constraints on glacial evolution (see Section 4.1.), there remain large uncertainties associated with the reconstruction of continental ice-sheets during MIS 6, across the PDG and during MIS 5e. In addition, interpretation that can be inferred from the measured tracer. We put a special emphasis on selecting marine sediment core records that can inform on millennial-scale changes occurring in the North Atlantic sector during HS11, potentially linked 15 with changes in NADW formation. This section ends with a review of the main limitations that will be associated with the model-data comparison and recommendations.

Available surface temperature syntheses
Quantitative comparisons between paleo-reconstructions and model outputs across the time interval 140-127 ka are possible, but remain limited to a few parameters (e.g. surface-air and sea-surface temperature (SAT and SST)). Qualitative and indirect 20 comparisons are also adequate to evaluate simulations: for example, simulated AMOC strength against paleo-records (section 8.2), simulated precipitation compared to Chinese speleothem calcite δ 18 O records or the simulated vegetation patterns against pollen-based vegetation reconstructions (section 8.3, Table 5).
There are currently no paleoclimate data compilations focusing on the penultimate glacial maximum (MIS 6) or covering the full length of the PDG, but syntheses of quantitative and temporal surface temperature changes focusing on the LIG are

North Atlantic records to inform on Heinrich stadial 11 and NADW changes
Deglaciation of NH ice-sheets led to increased meltwater input into the North Atlantic, thus reducing surface water density and potentially weakening NADW formation. It has been shown that changes in NADW have a significant impact on North Atlantic  Tzedakis et al. (2018) for the high-resolution calcite δ 18 O profile) as a reference for the PDG (Text S1 and S2, Table S1).

5
As a weakened NADW transport leads to the accumulation of remineralized carbon below 2500 m in the North Atlantic (Marchal et al., 1998), changes in oceanic circulation are also inferred from benthic foraminifera δ 13 C, a tracer for the ventilation of deep-water masses (Duplessy et al., 1988;Curry et al., 1988;Menviel et al., 2015;Schmittner et al., 2017). Finally, since a NADW weakening reduces the meridional heat transport to the North Atlantic (Stouffer et al., 2007), it is generally accompanied by a sea surface cooling across this region (Kageyama et al., 2013;Ritz et al., 2013). 10 The Norwegian Sea core MD95-2010 (Risebrobakken et al., 2006) indicates an increased IRD content starting at 139.5 ka (Fig. 8b). This could suggest an initiation of the penultimate deglaciation by enhanced iceberg calving (and subsequent melting) from the Fennoscandian ice-sheet. However, this has little to no effect in the North Atlantic. Only in the western side of the North Atlantic, right in the meanders of the Gulf Stream, does core CH69-K09 display a surface cooling and a drop in δ 18 Osw starting at ∼137. 5 ka (Labeyrie et al., 1999). Atlantic (Chapman and Shackleton, 1998;Martrat et al., 2007Martrat et al., , 2014Deaney et al., 2017) and a decrease in benthic δ 13 C below 20 2000 m depth (Labeyrie et al., 1999;Shackleton et al., 2003;Hodell et al., 2008;Deaney et al., 2017). The low-resolution ǫ N d and broadly corresponds to a pause in the sea-level rise. In parallel, Pa/Th in core ODP1063 displays a significant decrease at this time (not shown), and benthic δ 13 C increases in the deepest North Atlantic cores (Fig. 8). This warming could thus be due 30 to a short lived reinvigoration of NADW formation.
In agreement with previous studies (e.g. Marino et al., 2015;Tzedakis et al., 2018), the main phase of HS11 occurs between ∼133.3 and 130.4 ka and is characterised by a large IRD peak, cold surface-water conditions, a minimum in seawater δ 18 O , and low benthic δ 13 C values in most sites of the North Atlantic (Fig. 8). Elevated rates of iceberg calving and melting of NH ice-sheets thus lead to a large meltwater pulse in the North Atlantic, a significant global seal level rise (MWP2-B) (Grant et al., 2012, 2014), NADW weakening (Böhm et al., 2015 and a surface cooling that is at least regional. The resumption of NADW formation at the end of HS11 marks the end of the penultimate deglaciation (Tzedakis et al., 2012) (∼130.4-128.5 ka). It is characterised by an increase in SST, seawater δ 18 O (Oppo et al., 2006;Mokeddem and McManus, 2016;Martrat et al., 2007Martrat et al., , 2014Deaney et al., 2017) and benthic δ 13 C from the North Atlantic region (Fig. 8). The associated 5 atmospheric warming in the North Atlantic region also leads to terrestrial regrowth and thus a sharp increase in atmospheric CH 4 occurring between ∼129 and 128.5 ka (Landais et al., 2013).  Govin et al. (2015). It is also important to recognise that the dominant drivers of speleothem δ 18 Oc may change over time and differ from one cave site to another (see Section 8.4). Therefore we strongly advise a critical assessment of these interpretations based on the most recent developments and advances in stable-isotope hydrology and the 15 original publications.

Other environmental and climate reconstructions
Tables 4 and 5 also detail the timing of the major changes (and associated uncertainties) as recorded in each timeseries. These major changes are identified using the RAMPFIT or BREAKFIT software programs (Mudelsee, 2000(Mudelsee, , 2009. Age uncertainties (1σ) reported in Tables 4 and 5 include i) the "internal" error of the event given by RAMPFIT or BREAKFIT, ii) the relative error related to the climatic alignment method for marine sediment and Ioannina records and iii) the absolute dating error 20 of the reference time scale. The selected ice-and sediment-based records are all displayed on AICC2012 or on a timescale coherent to this ice-core age scale, while terrestrial records are displayed on their own independent chronology. Among these paleoclimatic time series, several isotopic records are presented here and we strongly encourage transient simulations to be performed with oxygen and/or carbon isotope-enabled models to allow direct quantitative comparison between simulated and measured isotopic time series. 25 Little change occurs until the beginning of phase 1 at ∼136.4 ka, after which a cooling phase is identified in a few records of the North Atlantic ( Fig. 8 and Fig. 9a, major change [A]). This also corresponds to reduced monsoon activity as recorded in Chinese speleothems (Fig. 9d [M]), and the initiation of the Antarctic warming (Fig. 9h [W]). The short-lived warming event in the North Atlantic associated with phase 2 (Fig. 9a [B]) is also identified in the Chinese speleothems as a slightly wetter interval (Fig. 9d [N]). Other environmental records might not have the necessary resolution to record this multi-centennial-scale 30 event. The main phase of HS11, corresponding to phase 3, is associated with meltwater input and cold conditions in the North Atlantic ( Fig. 9a [C] and 9c [I]), dry conditions over Europe (Fig. 9b) and Asia (Fig. 9d,

interval between [O] and [P]), and
warmer conditions at high southern latitudes (Fig. 9f, h). The end of HS11 (phase 4) associated with a pause in the meltwater input ( Fig. 9c [J]) and progressively warmer conditions in the North Atlantic and southern Europe (Fig. 9a [D corresponds to a strengthening of the Asian monsoon ( Fig. 9d [P] and 9e [Q]), and maximum warmth at high southern latitudes ( Fig. 9f [R, S] and 9h [X]). Interglacial conditions in atmospheric CO 2 and CH 4 as well as North Atlantic temperatures and ventilation are attained at about 128.5 ka (Fig. 8), which is also associated with warm and wet conditions in southern Europe ( Fig. 9b [G], and 9c [L]).

Limitations and recommendations 5
One important consideration to account for when comparing simulated variables against paleodata between 140 and 127 ka, is the large uncertainties associated with both absolute and relative chronologies of most paleoclimatic records during this time interval. Dating uncertainties range from a few centuries to up to several thousand years depending on the type of archive and dating methods (Govin et al., 2015). For instance, the average absolute dating error attached to the Corchia Cave δ 18 O record is ∼0.7 ka (2σ) (Tzedakis et al., 2018). For marine sediment chronologies, which are mainly based on record alignment 10 strategies (e.g. a record on a depth scale is aligned onto a dated reference record), the age uncertainties encompass a relative dating error ("alignment error") in addition to the absolute error attached to the chronology of the dated reference record.
As a result, the overall 2σ age error associated with North Atlantic sediment core ODP976 aligned onto the Corchia record is 1.6 ka on average (Table S2). Between 140 and 127 ka, the average absolute error attached to the AICC2012 chronology used to display ice and gas records from the EDC ice core is about 4 ka (2σ) (Bazin et al., 2013;Veres et al., 2013). This large 15 AICC2012 dating uncertainty is thus attached to the GHG concentration records used to force the transient simulations (Section 3). It will therefore taint the relative timing of the changes in orbital and atmospheric CO 2 forcing that will largely drive the simulated evolution of climate and environmental changes across the PDG. The relative timing of those changes will also be affected by uncertainties attached to the temporal evolution of the ice-sheets and related meltwater forcing scenarios. However, uncertainties attached to the relative timing of changes between GHG forcing and the simulated changes in Antarctica and 20 the Atlantic Ocean basin are somewhat reduced since the ice-and sediment-based records are also displayed on AICC2012 timescale or on time scales coherent to the reference ice-core age scale (Section 8.2).
Limitations are also attached to the potential misinterpretation of climate and environmental proxies due to incomplete understanding of how some of those archives record climatic and environmental change. First, SST records highlighted here have been reconstructed using various microfossil and geochemical methods. Although the use of various tracers is known to 25 yield SST discrepancies in particular above 35 • N (MARGO project members, 2009), the extent to which these different SST reconstruction methods influence the representation of temporal climatic changes across our studied time interval is poorly known. Additional difficulties arise from the individual methods commonly used to reconstruct past SST due to, for example, the poor understanding of the modern habitat (e.g. living season and water-depth) of microfossil species (e.g. foraminifera) (Jonkers and Kučera, 2015) and alkenone producers (e.g. Rosell-Melé and Prahl, 2013). With these limitations acknowledged, ). In addition, Sime et al. (2009) suggested that the temperature at Dome C should be much higher than the one inferred from water isotopes and a constant temperature/δ 18 O slope during MIS 5e. For now, we suggest the use of simulated annual mean climatic variables for the comparison. However, it is crucial that the seasonality of paleo-records is better assessed to improve the interpretation of temperature reconstructions, and hence the model-data comparisons.
In addition, uncertainties remain regarding the dominant controlling factors (i.e. changes in temperature, rainfall amount 5 and rain sources) of speleothem calcite isotopic records (e.g. Govin et al., 2015). For instance, δ 18 Oc records throughout Asia are commonly interpreted as tracers of past changes in the intensity of the Asian monsoon. In particular, Chinese speleothem δ 18 Oc is classically interpreted as reflecting the East Asia monsoon (e.g. Cheng et al., 2009Cheng et al., , 2016. However, a water-hosing experiment performed with an oxygen-isotopes-enabled model suggests instead that δ 18 Oc variations may reflect changes in the intensity of the Indian, rather than East Asian, monsoon precipitation during Heinrich events (Pausata et al., 2011;10 Caley et al., 2014). Another recent model study also demonstrates that variations in ice core and speleothem oxygen isotope reconstructions cannot solely be attributed to climatic effects, but also reflect depleted δ 18 O of nearby oceans during glacial meltwater events (Zhu et al., 2017). Overall, we strongly encourage the use of isotope-enabled models to allow for direct and quantitative model-data comparison of isotopic tracers.
Additional multi-centennial-scale paleoclimatic reconstructions are necessary in order to further constrain the millennial 15 scale variability during the PDG. It is also crucial that comprehensive data compilation work is carried out over the entire studied time interval and covering, in particular, the penultimate glacial maximum in order to test the robustness of the initial 140 ka spin-up climate. Modelling groups running transient deglacial simulations, and/or associated sensitivity experiments, are encouraged to use multiple paleorecords for a full diagnosis of the simulations.

20
Here, we present a protocol for performing transient simulations of the PDG spanning 140 to 127 ka. Changes in boundary conditions across the PDG that will serve as a forcing are presented and discussed. This includes changes in orbital parameters, GHG concentration, NH and Antarctic ice-sheets, and associated deglacial meltwater input. While not used as a direct forcing, changes in global sea-level are also presented on a new chronology. Finally, a series of key paleoclimatic and paleoenvironmental records are suggested to perform model-data comparisons. The marine records were recovered from the North Atlantic 25 and Southern Ocean, while the continental records were retrieved from Europe, China and Antarctica. Performing transient simulations with oxygen-and/or carbon-isotope-enabled Earth-system models could significantly improve model-data comparisons by providing a more direct and quantitative comparison with paleoproxies based on measured isotopic signatures (e.g. Simulations of the penultimate deglaciation would allow a comparison with the last deglaciation, therefore highlighting 30 similarities and differences between the last two deglaciations. The evolution of insolation across the two deglaciations is different, potentially explaining the relatively more rapid disintegration of NH ice-sheets during the PDG compared to the last deglaciation. Acting both as a response to and driver of changes, atmospheric CO 2 appears to increase much more gradually EPICA community members: One-to-one coupling of glacial climate variability in Greenland and Antarctica, Nature, 444, 195-198, 2006.  (Schilt et al., 2010;Spahni et al., 2005) Ice-sheets*   Table 3. Requested 1-dimensional greenhouse gas concentrations used as forcing, as well as 2-dimensional ((X,Y) or (Y,Z)) and 3-dimensional (X,Y,Z) variables to be uploaded to the ESGF database. Most variables are requested as annual means, but * indicates that monthly outputs are requested. The CMIP6 conventions should be used for the variable names and units. See https://earthsystemcog.org/projects/wip/CMIP6DataRequest for details. Table 4: Key paleoproxy records of changes in ocean oxygen (δ 18 O) and carbon (δ 13 C) isotopic composition. The chronology of all paleoproxy records presented here is based on an alignment onto Corchia U-Th-based chronology (Table S1). Letters in brackets indicate the major changes identified in the paleoclimatic records shown in Figure 9. The dates of the major changes obtained through RAMPFIT or BREAKFIT are indicated with their associated uncertainties. Age uncertainties (1σ) reported include (1) the "internal" error of the event given by RAMPFIT (for most major changes) (Mudelsee, 2000), or BREAKFIT have been split into five phases (ϕ), with the date representing the beginning of each phase: ϕ1 is associated with the early phase of HS11 and MWP-2A, ϕ2 represents a pause within HS11, ϕ3 the main phase of HS11 and MWP-2B, ϕ4 the inception out of HS11 and ϕ5 full interglacial conditions. Implicitely the end of each phase corresponds to the beginning of the next one. 10 For the sea-level record, the new chronology discussed in section 5 is used. CS98 refers to Chapman and Shackleton (1998).

North American and
ventil. stands for ventilation.   Table 4 for key paleoproxy records of climatic and environmental changes through the PDG selected for their relatively high resolution. An indication of their climatic or environmental use is indicated in italic and the proxy used is shown in parenthesis in column 1. Age uncertainties (1σ) reported include (1) the "internal" error of the event given by RAMPFIT (for most major changes) (Mudelsee, 2000), or BREAKFIT (major change dates highlighted with **) (Mudelsee, 2009), (2) the relative error related to the climatic alignment method for marine sediment and Ioannina records and (3) the 5 absolute dating error of the reference time scale. • No dating error is provided in Yang and Ding (2014). As a result, the stated error here only encompasses the "internal" error, it thus represents only a minimal error for the timing of the stacked grain size increase. The mean age representing the beginning of each phase (ϕ) is shown in the last row and is calculated from all the available dates within each phase shown in Tables 4 and 5. FFA stands for foraminifera faunal assemblages. CS98 refers to Chapman and Shackleton (1998). Based on North Atlantic records, changes have been split into five phases (ϕ), with the date 10 representing the beginning of each phase: ϕ1 is associated with the early phase of HS11 and MWP-2A, ϕ2 represents a pause within HS11, ϕ3 the main phase of HS11 and MWP-2B, ϕ4 the inception out of HS11 and ϕ5 full interglacial conditions.

Intensity of
Implicitely the end of each phase corresponds to the beginning of the next one.   and SU90-03 (blue) (Chapman et al., 2000). b) Percentage of temperature tree pollen from the Greek Ioannina sequence (dark green) (Tzedakis et al., 2003) (Cheng et al., 2016). e) Staked grain size record of Chinese Loess (light grey) (Yang and Ding, 2014). f) SST records from Southern Ocean cores MD97-2120 (light blue) (Pahnke et al., 2003) and MD02-2488 (dark blue) (Govin et al., 2012). g) δ 13 C record from North Atlantic core MD95-2042 (light pink) (Shackleton et al., 2003) and Southern Ocean core MD02-2488 (dark blue) (Govin et al., 2012). h) Antarctic ice δD record from EPICA Dome C (Jouzel et al., 2007) on AICC2012 chronology (black) (Bazin et al., 2013;Veres et al., 2013). Letters in brackets and yellow crosses indicate the major changes identified in the paleoclimatic records, and their ±1σ dating uncertainty (horizontal grey bars, Tables 4 and 5). Chronologies and associated references are detailed in Tables 4 and 5. No details on the dating errors attached to the Chinese Loess stacked grain size record are provided in (Yang and Ding, 2014), the 1σ error should thus be treated as an underestimated dating error.