Factors controlling the last interglacial climate as simulated by LOVECLIM 1 . 3

The last interglacial (LIG), also identified to the Eemian in Europe, began at approximately 130 kyr BP and ended at about 115 kyr BP (before present). More and more proxy-based reconstructions of the LIG climate are becoming more available even though they remain sparse. The major climate forcings during the LIG are rather well known and therefore models can be tested against paleoclimatic data sets and then used to better understand the climate of the LIG. However, models are displaying a large range of responses, being sometimes contradictory between them or with the reconstructed data. Here we would like to investigate causes of these differences. We focus on a single climate model, LOVECLIM, and we perform transient simulations over the LIG, starting at 135 kyr BP and run until 115 kyr BP. With these simulations, we test the role of the surface boundary conditions (the time-evolution of the Northern Hemisphere (NH) ice sheets) on the simulated LIG climate and the importance of the parameter sets (internal to the model, such as the albedos of the ocean and sea ice), which affect the sensitivity of the model. The magnitude of the simulated climate variations through the LIG remains too low compared to reconstructions for climate variables such as surface air temperature. Moreover, in the North Atlantic, the large increase in summer sea surface temperature towards the peak of the interglacial occurs too early (at∼ 128 kyr BP) compared to the reconstructions. This feature as well as the climate simulated during the optimum of the LIG, between 131 and 121 kyr BP, does not depend on changes in surface boundary conditions and parameter sets. The additional freshwater flux (FWF) from the melting NH ice sheets is responsible for a temporary abrupt weakening of the North Atlantic meridional overturning circulation, which causes a strong global cooling in annual mean. However, the changes in the configuration (extent and albedo) of the NH ice sheets during the LIG only slightly impact the simulated climate. Together, configuration of and FWF from the NH ice sheets greatly increase the magnitude of the temperature variations over continents as well as over the ocean at the beginning of the simulation and reduce the difference between the simulated climate and the reconstructions. Lastly, we show that the contribution from the parameter sets to the climate response is actually very modest.


Introduction
The last interglacial (LIG), the previous period of low continental ice volume before the Holocene, coincides with marine isotope stage 5e (MIS5e) and is also identified, although not strictly (Sánchez Goñi et al., 1999), to the Eemian, the interval of optimal development of vegetation in Europe.Broadly speaking, it covers the time interval 130-115 kyr BP (before present).Published paleoclimatic record compilations provide geographical pattern of average temperature changes during the LIG both over the continents and the surface of the oceans.The global annual mean surface temperature during this period was 0-2 • C above the pre-industrial one (e.g., CLIMAP Project Members, 1984;Bauch and Erlenkeuser, 2003;Turney and Jones, 2010) and the global Published by Copernicus Publications on behalf of the European Geosciences Union.
annual mean sea surface temperature (SST) was 0.7 ± 0.6 • C higher than during the late Holocene (McKay et al., 2011).The Arctic in summer was about 5 • C warmer than present (CAPE-Last Interglacial Project Members, 2006;Kopp et al., 2009), and estimated summer temperatures over Europe were ∼ 2 • C higher than present (Kaspar et al., 2005).
Temporal reconstructions show that the Arctic reached its peak warmth soon after 130 kyr BP (during the first part of the LIG) (CAPE-Last Interglacial Project Members, 2006).The surface temperature at North Greenland Eemian Ice Drilling -NEEM -site (Northwest Greenland; camp position 77.45 • N 51.06 • W) peaked at 8 ± 4 • C above the mean of the past millennium at 126 kyr BP (NEEM community members, 2013).In northern Europe, the LIG can be divided into three parts: an early optimum, followed by a slight cooling and eventually a sharp drop (Brewer et al., 2008;Rioual et al., 2001).The rapid warming trend over western Europe started at the penultimate deglaciation and continued until 130 ± 2 kyr BP, with a temporary (shorter than 600 yr) cold event ca.131.4 ± 2 kyr BP (Sánchez Goñi et al., 2012).Antarctica also displays a rapid initial warming starting from 135 kyr BP (or before) until about 130 kyr BP, followed by a slow cooling (Masson-Delmotte et al., 2011).By contrast, in southern Europe, temperatures remained stable during the LIG.Depending on sites and proxies, SSTs in the Norwegian Sea and the North Atlantic peaked in the first (soon after 130 kyr BP) or second part (later than 125 kyr BP) of the LIG (Cortijo et al., 1999;Oppo et al., 2006;Martrat et al., 2007;Bauch and Kandiano, 2007;Guihou et al., 2010;Van Nieuwenhove et al., 2011;Govin et al., 2012).
There is some evidence that the global ocean circulation during the LIG was similar to the present circulation (Duplessy et al., 2007;Duplessy and Shackleton, 1985;Evans et al., 2007), but other studies suggest a high production rate of North Atlantic Deep Water (NADW) during MIS5e (Yu et al., 1996;Adkins et al., 1997).Moreover, some variations are also recorded in the strength of the ocean circulation through the LIG, indicating a maximum of the Atlantic meridional overturning circulation (AMOC) in the second part of the LIG (Oppo et al., 1997;Sánchez Goñi et al., 2012).Several studies also suggest that there was no active site of deep water formation in the Labrador Sea during the LIG (Hillaire-Marcel, 2001;Rasmussen et al., 2003a;Winsor et al., 2012).
The detailed history of sea level changes during the LIG is still under debate.Sea level was 16 to 18 m below the present one at 135 kyr BP (Gallup et al., 1994(Gallup et al., , 2002)).From 130 ± 2 kyr BP until the glacial inception at the end of the LIG, ca.116 kyr BP, the sea level reached at least 4 m above modern level (Tarasov and Peltier, 2003) but unlikely more than +9 m (Stirling et al., 1998;Henderson and Slowey, 2000;McCulloch and Esat, 2000;Gallup et al., 2002;Muhs et al., 2002;Hearty and Neumann, 2001;Overpeck et al., 2006;Kopp et al., 2009;Dutton and Lambeck, 2012;O'Leary et al., 2013).Then it fell to a low stand, as low as -19m at 113 ± 0.4 kyr BP (Gallup et al., 2002;Cutler et al., 2003).Even if the melting and growing of the Antarctic and Greenland ice sheets may have partly contributed to these sea level changes (Alley et al., 2010;Helsen et al., 2013;Koerner, 1989;Koerner and Fisher, 2002;Colville et al., 2011), other continental ice sheets have played a significant role as well (Lambeck et al., 2006;Kopp et al., 2009).Moreover, sea level higher than the present sea level stand is almost exclusively due to Antarctic and Greenland ice sheets melting, while all the ice sheets likely contributed to sea level lower than present sea level stand.
All these studies illustrate that climate changes over the LIG are complex, with variability in time and space.Although more and more climate reconstructions from proxy data become available, they are still sparse.Also, discussing sequences of climatic events is still limited by the lack of robust and coherent chronologies between records from various paleoclimatic archives despite several efforts to display data from different archives on a same chronology (Drysdale et al., 2007;Govin et al., 2012;Capron et al., 2014).
All these paleoclimatic records provide useful information for comparison to climate model simulations.In particular, few transient climate simulations have been performed over the LIG.For example, both MoBidiC, an Earth system model of intermediate complexity (Crucifix and Loutre, 2002), and the MPI/UW 3D-Earth system model (Gröger et al., 2007) simulated a cooling over continents, in particular at high latitudes, and a decrease of the high latitude SST starting from the beginning of the simulation at 126 or 128 kyr BP, respectively.None of these models were coupled with an ice sheet module.Both studies pointed out the role of the astronomical forcing in this cooling but also the importance of feedback related to sea ice, snow and vegetation.This is confirmed by Calov et al. (2005) with the CLIMBER-2 climate model coupled to the SICOPOLIS ice sheet model.Crucifix and Loutre (2002) using MoBidiC also demonstrated a slight weakening of the AMOC and an acceleration of the cooling trend between 122 and 120 kyr BP.Ritz et al. (2011) showed that the simulated AMOC weakens rapidly at terminations (the transition from a glacial to the following interglacial) in response to changes in the ocean surface density due to freshwater relocation from melting ice sheets to the ocean.Their simulation, performed with the Bern coupled three-dimensional dynamical oceanenergy balance atmosphere model taking into account ice sheet changes through time, was recently included in an intercomparison of transient surface temperature changes during the LIG simulated with seven climate models (Bakker et al., 2013).Robust features for all the models were identified such as the timing of the maximum of the summer temperature in the Northern Hemisphere (NH) (occurring between 130-122 kyr BP) and the timing of the maximum of simulated January temperature over the Arctic Ocean (occurring between 128-126 kyr BP).The role of feedback related to sea ice was also highlighted.However, other features, such as the strength of the AMOC, appeared model-dependent.
According to Bakker et al. (2013), taking into account the remnant NH continental ice sheets delays the timing of the maximum of July temperature, both in the NH and SH (Southern Hemisphere).
From a number of sensitivity equilibrium simulations, Bakker et al. (2012) concluded that the Greenland ice sheet (GIS) melting during the LIG is the major cause of shutdown of deep convection in the Labrador Sea and of a significant cooling of the northern North Atlantic Ocean.Furthermore, reduction in surface elevation and extent of the GIS may amplify the reduction in deep convection and the cooling in the Nordic Seas.In an intercomparison of snapshot simulations at 130, 128 and 125 kyr BP performed with several models of different levels of complexity, Lunt et al. (2013) noted that, although some robust features of temperature change were identified, the models generally underestimate the magnitude of the response derived from proxies.Langebroek and Nisancioglu (2014), who performed equilibrium simulations at 130, 125, 120 and 115 kyr BP with the Norwegian Earth System Model, pointed out the role of the insolation in the seasonal cycle of temperature and the importance of the greenhouse gas forcing in the actual value of temperature.They also showed that temperature is maximum during the early LIG in the NH in June-July-August (JJA) and south of 45S in JJA and December-January-February (DJF), while it reaches its maximum late in the LIG in the NH in DJF.
Here, we analyse transient simulations performed over the LIG, starting at 135 kyr BP and run until 115 kyr BP, by a single climate model, LOVECLIM.Our main purpose is to focus on the uncertainties in the simulated climate.In particular, we tackle the question of whether the uncertainties on simulated climate may explain the divergence between the reconstructed and simulated climates.Although very important as well, the issue of uncertainties on proxy data is out of the scope of this paper and will only be very quickly browsed through.Broadly speaking, the origin of the uncertainty in the simulated climate may be related to the model design and to the forcings and surface boundary conditions used.We identify how much the imperfect knowledge of the values of several physical parameters of the climate model induces changes in the simulated climate.The presence and melting of the NH ice sheets is another source of uncertainty.In the absence of accurate reconstructions of the evolution of those ice sheets during the LIG, the impact of the changes in their extent and altitude and in the amount of freshwater entering the ocean due to their waning is evaluated through sensitivity experiments.In all these experiments, we mostly focus on temperature changes, in particular their magnitude and the timing of their maximum.Nevertheless, the role of feedback is not ignored and variables other than temperature are also discussed.

Methodology
2.1 LOVECLIM1.3LOVECLIM1.3(further termed LOVECLIM) is a threedimensional Earth System Model of Intermediate Complexity (EMIC).It consists of five components representing the atmosphere (ECBilt), the ocean and sea ice (CLIO), the terrestrial biosphere (VECODE), the oceanic carbon cycle (LOCH) and the Greenland and Antarctic ice sheets (AG-ISM).The ice sheet component and the ocean carbon cycle model are not activated in the present study.The reader is referred to Goosse et al. (2010) for a more detailed description of these components, as well as to Mouchet and François (1996) for LOCH and to Huybrechts (2002) for AG-ISM.
ECBilt (Opsteegh et al., 1998) is a spectral T21 threelevel quasi-geostrophic atmospheric model that explicitly computes synoptic variability associated with weather patterns.The T21 truncation corresponds to a resolution of about 5.625 • × 5.625 • .This model includes simple parameterisations of the diabatic heating processes and an explicit representation of the hydrological cycle.Cloudiness is prescribed according to present-day climatology.A modified atmospheric balance equation is used (Sriver et al., 2013) to improve the representation of the atmospheric tropical dynamics in the most recent version of LOVECLIM.A parameterisation of katabatic winds was also implemented following Barthélemy et al. (2012).
VECODE (Brovkin et al., 2002) is a reduced-form model of vegetation dynamics, which simulates the dynamics of two plant functional types (trees and grassland) at the same resolution as that of ECBILT.
sensitivity is in the lower part of the range of the other models and its dynamical response is weak.

Insolation and greenhouse gas forcings
All the transient simulations performed in this study start at 135 kyr BP from an equilibrium state at that time and are run for 20 kyr, until 115 kyr BP.They are all forced by timedependent changes in insolation and greenhouse gas (GHG) concentrations.The model year is divided into 12 months of 30 days.The vernal equinox occurs on 21st of March.
The changes in the distribution of insolation (Fig. 1) received by the Earth are computed from the changes in orbital configuration (Berger, 1978).June insolation at all latitudes (except polar latitudes in the SH) reaches a maximum at about 127 kyr BP.Then it decreases to a minimum at 116 kyr BP.December insolation follows an almost opposite pattern for all latitudes (except northern polar latitudes).It increases starting from 127 kyr BP until ∼ 116 kyr BP.The June insolation is more than 10 % larger than the present-day value for most of the latitudes, while the December insolation is more than 10 % smaller than the present-day value for most of the latitudes.The largest anomalies occur in the polar regions during their local summer.The timing of the insolation changes is mostly driven by the changes in the climatic precession and its amplitude is large due to large values of the eccentricity.The obliquity also plays a role in the large changes in the amplitude of the insolation in particular in the polar regions.
The GHG forcing (Fig. 2) (Petit et al., 1999;Pépin et al., 2001;Raynaud et al., 2005;Loulergue et al., 2008;Spahni et al., 2005) is prescribed and the concentrations of the different gases are time-dependent.The atmospheric CO 2 , CH 4 and N 2 O concentrations of the PMIP3 protocol (https://pmip3.lsce.ipsl.fr) is used from 132 to 115 kyr BP.Before 130 kyr BP, the atmospheric CO 2 concentration experiences a fast increase by more than 50 ppm in 5 kyr.Then, it remains between 260 and 280 ppm until 115 kyr BP, except during a short time interval (128.9 to 128.3 kyr BP), when its values are higher than 280 ppm.The N 2 O concentration displays a similar behaviour, with a rapid increase from 135 to 130 kyr BP, followed by almost stable values, with a low variability.By contrast, the CH 4 concentration shows a different behaviour.It reaches a clear maximum at about 128 kyr BP.This maximum is preceded by a fast increase and followed by a slow decline.The maximum atmospheric concentrations of these GHGs during the LIG are similar to the pre-industrial ones.

Surface boundary conditions
The Antarctic ice sheet is kept fixed to its pre-industrial state and there is no additional freshwater flux in the SH.Instead, the time evolution of the NH ice sheets is prescribed to several different configurations.In a first set of simulations, the Greenland ice sheet extent and surface elevation are fixed to their pre-industrial values through the whole simulation, and without any other ice sheet in the NH (PI ice sheet in Table 1).The presence of remnant NH ice sheets at the onset of the transient experiments at 135 kyr BP makes it necessary to complement these boundary conditions.For instance, as ice sheet constraints on the penultimate glacial termination (Termination II) are sparse, the NH ice sheet retreat from the last deglaciation has been re-mapped based on other available information for Termination II.It is assumed that, if the global sea level at a time during termination II is similar to that at a given moment of termination I, then the NH ice sheet configuration is the same for both times.This method (used  (Petit et al., 1999;Pépin et al., 2001;Raynaud et al., 2005); (top right) CH 4 (Loulergue et al., 2008); and (bottom left) N 2 O (Spahni et al., 2005).e.g. by Ritz et al., 2011), which implicitly assumes that NH and SH land-ice cover evolved in phase during both periods, has been extended by relying on a recently updated sea level record (Grant et al., 2012) and by using an alternative NH ice sheet reconstruction.For the latter, a large number of geomorphological constraints on ice sheet extent have been digitised from the literature and ice sheet surface elevation has been reconstructed by assuming a plastic flow relation for the ice sheets.The resulting boundary conditions consist of a chronology of ice mask and surface elevation over the entire LIG period (called GR in Table 1), later referred to as the NH ice sheet configuration.Ice sheet albedo is also applied in this simulation.According to this reconstruction, the Fennoscandian ice sheet is fully melted at 132.1 kyr BP, while the North American ice sheet persists until 131.3 kyr BP.At the end of the interglacial, the regrowth of the ice sheet starts at 120.2 kyr BP over North America and at 117.4 kyr BP over northern Europe.
NH ice sheet melting is accompanied by an increased flow of freshwater (FW) into the surrounding oceans, which could potentially alter the strength of the AMOC, and hence the global climate.The NH FW forcing was estimated based on the same re-mapping as used for the elevation and ice mask reconstructions.The spatial distribution of produced melt water is estimated from the reconstructed elevation changes and a continental runoff routing model (Goelzer et al., 2012b) is used to identify the location and magnitude of meltwater fluxes to the ocean.The FW flux to four regions is distinguished: Pacific (PAC), Arctic Ocean (ARC), North Atlantic (NAT) and North Sea (NOS).In order to avoid too large FW input into the coastal ocean model grid points, the freshwater flux is integrated over each region and equally redistributed to the points that show a magnitude higher than 0.01 Sv any time over the period 135 to 125 kyr BP.The total freshwater flux magnitude is scaled to be in line with the total sea level contribution of the NH ice sheets (Fig. 3; GR in Table 1).
The glacial-interglacial contrast in NH sea-level contribution between penultimate glaciation and LIG is assumed to be the same as between Last Glacial Maximum (LGM) and present day, for which a value of 110 m Sea level equivalent (SLE) is taken, in line with the reference model results of Zweck and Huybrechts (2005).
Alternatively, the changes in the NH ice sheet configuration and in the additional freshwater flux in the NH are computed from the sea level record estimated from Lisiecki and Raymo (2005) (LR in Table 1).Sea level reconstruction from LR lags behind the one from GR by about 2.8 kyr during the LIG.Therefore, both reconstructed melting and glacial inception are delayed in LR compared to GR.Indeed, in LR scenario, the Fennoscandian ice sheet is fully melted at 131 kyr BP and the North American ice sheet persists until 127.5 kyr BP, while glacial inception starts at 116.5 kyr BP over North America and even later over northern Europe in LR.Moreover, the FW flux magnitude has been smoothed in time so that peaks in the FW flux are not so high in the LR reconstruction.Loutre et al. (2011) and Goosse et al. (2007) identified several parameter sets that provide modified model responses to increased atmospheric CO 2 concentration and to freshwater hosing.Amongst the parameters involved in these parameter sets, there are the albedos of the ocean and sea ice, and parameters involved in the equation of the quasi-geostrophic potential vorticity, in the formulation of the long-wave radiative scheme, in the representation of the vertical diffusivity of the ocean and in the computation of the Coriolis term in the equation of motion.A detailed description of the different parameter sets is available in Loutre et al. (2011).In addition to the reference parameter set (std), we also use parameter set 22 (as defined and used in Loutre et al. (2011), Goelzer et al. (2011) and Goosse et al. (2007); see also supplementary material) to test the response of LOVECLIM1.3.

Parameter sets
The global climate sensitivity with both parameter sets is 2.1 • C. It is computed as the global annual mean surface temperature increase after 1000 yr from the pre-industrial value in a simulation for which the atmospheric CO 2 concentration is increased by 1 % per year from the pre-industrial value until doubling is reached (after 70 yr) and held constant thereafter.The northern polar latitude sensitivity is stronger than the southern one with parameter set std, while they are almost the same for parameter set 22. The polar latitude sensitivity is computed as the annual mean surface temperature increase poleward of 75 • in either hemisphere after 1000 yr from the pre-industrial value in the same simulation as for the computation of the global climate sensitivity.
The AMOC sensitivity to a freshwater hosing is −24 % for parameter set std and −41 % for parameter set 22. The AMOC sensitivity index is computed as the percentage of decrease in the maximum value of the meridional overturning stream function below the Ekman layer in the Atlantic Ocean after 1000 yr in a water hosing experiment.In this experiment, freshwater is added in the North Atlantic (20-50 • N) with a linearly increasing rate of 2 × 10 −4 Sv yr −1 .and Table 1).The series are smoothed using a moving average over 100 yr.The black dot on the right hand side of the figures provides the corresponding simulated pre-industrial value.

General description of the simulated climate
The reference simulation (allGR) is forced by orbital and GHG changes and takes into account changes of the NH ice sheets from both their configuration and additional freshwater fluxes to the ocean (see Table 1 for a summary of the experiments).This results in a decrease in the simulated global annual mean surface temperature, global annual mean SST, Arctic (north of 69 • N) summer temperature and summer temperatures over Europe (the region bounded by 17 • W-11 • E and 36-64 • N), averaged between 135 and 115 kyr BP, by 0.2, 0.2, 0.2 and 0.3 • C, respectively, compared to the simulated pre-industrial values.
Broadly speaking, the simulation can be divided into an interval of transition from the previous glacial state (at the beginning of the simulation until ∼ 131 kyr BP), a plateau, and an interval of transition towards glacial inception (at the end of the simulation, after ∼ 121 kyr BP) (Fig. 4a).

Out of the glacial state
The input of additional freshwater in the North Atlantic before 131 kyr BP leads to a strong weakening of the AMOC (Fig. 4b), and even to its almost complete collapse for a few centuries around 132 kyr BP.However, as soon as the additional freshwater input stops the AMOC recovers.Then, although the AMOC is less intense than under pre-industrial conditions, deep convection still takes place in the Labrador Sea (Fig. 5), contrarily to some reconstructions (Hillaire-Marcel, 2001;Rasmussen et al., 2003a;Winsor et al., 2012) and convection is very active in the Norwegian and Barents seas.In response to the near collapse of the ocean circulation at ∼ 132 kyr BP, the global annual mean surface temperature (Fig. 4a) experiences a rapid decrease of 0.8 • C followed by a fast increase of 2.1 • C, all occurring in less than 1500 yr (between 132.5 and 131 kyr BP).Global annual mean surface temperature before 132 kyr BP remains lower than the pre-industrial one.Simultaneously, the simulated North Atlantic SST and western France and Greenland surface temperatures (Fig. 7a-d and g) experience a decrease of several degrees.However, at the same time, the simulated tropical South Atlantic SST (Fig. 7e) exhibits a slight increase.This is consistent with the so-called seesaw mechanism studied by Crowley (1992) and Stocker et al. (1992).Moreover, the simulated high latitude South Atlantic SST (Fig. 7f) and the simulated surface temperature at the Dome C site (Antarctica) (Fig. 7h) show almost no variations.

The plateau
The long-term evolution of the global annual mean surface temperature shows an increase from 131 kyr BP until the maximum of 16.6 • C at 128.5 kyr BP (Fig. 4a).The maximum over the LIG of the simulated global annual mean surface temperature (global annual mean SST, Arctic summer temperature, summer temperature over western France, respectively) is 0.7 • C (0.4, 2.5, 4.4 • C, respectively) higher than the corresponding pre-industrial value.These values are within the range of the reconstructed estimates for largescale changes (see section 1).The global annual mean cooling after the maximum at 128.6 kyr BP consists of three phases.First, there is a rapid cooling until ∼ 126.8 kyr BP.Then, temperatures remain stable or even slightly recover until 124.8 kyr BP.Lastly, cooling proceeds at fast rate until ∼ 121 kyr BP.However, the simulated temperature remains above the simulated PI values (15.8 • C) from 131.3 to 120.9 kyr BP.The AMOC (Fig. 4b) display no significant trend during this period, although its variability is higher at the beginning than at the end of the sub-interval.

Towards glacial inception
After 121 kyr BP, when the NH ice sheets start to regrow, the simulation shows a high variability of the AMOC (Fig. 4b).We hypothesise that the model becomes close to a bistable regime, which makes it oscillating between two modes as previously discussed in Jongma et al. (2007).Such type of bistability has been found in various versions of LOVECLIM and for different forcing conditions (Goosse et al., 2002;Friedrich et al., 2009;Jongma et al., 2007;Rahmstorf et al., 2005).It has not been formally tested for the present version of the model and LIG conditions.As previous experiments showed that this behaviour is strongly model dependent, we should analyse this feature with care.The higher variability of the AMOC is mirrored in the global annual mean temperature (Fig. 4a).The standard deviation is then ∼ 0.18 • C, while it was ∼ 0.08 • C between 125 and 122 kyr BP.This simulated high variability is also displayed in the temperature at many high northern latitude locations (Fig. 7a-d).

Presentation of the proxy data
We selected several air and sea surface temperature reconstructions from ice and marine sediment cores for comparison with the simulated temperature (Fig. 6, Table 2).The temperature time series we have considered have been reconstructed using foraminifera assemblage transfer functions for the SST time series from marine sediment cores and water isotopic time series for the surface air temperature reconstructions from the EPICA Dome C (EPICA Community Members, 2004) and NEEM (NEEM community members, 2013) ice cores.Uncertainties on each reconstructed SST record estimated from (1) the uncertainty on measurement and (2) the calibration of geochemical and microfossil proxies under modern conditions range between 0.6 and 1.9 • C depending on the SST proxies (e.g.Chapman and Shackleton, 1998;Oppo et al., 2006).It reaches up to 4 • C for the NEEM ice core precipitation-weighted temperature reconstruction (NEEM community members, 2013).
Absolute dating constraints over the LIG are sparse and this results in discrepancies (up to several thousands of years) between the various chronologies that are defined for paleoclimatic records.Those discrepancies highly depend on the strategy followed to define the original age model and the chosen reference curve (Capron et al., 2013).This is particularly critical when comparing climatic records from different paleoclimatic archives.For our model-data comparison, we take advantage of a new coherent temporal framework over the LIG between multiple ice core and marine sediment records defined by Capron et al. (2014).In par- ticular, we use several paleoclimatic reconstructions from ice core and marine cores transferred onto the recent ice core AICC2012 chronology (Bazin et al., 2013;Veres et al., 2013).The absolute uncertainty across the LIG of the ice core AICC2012 timescale is less than 2 kyr.Due to the synchronisation method of marine sediment records onto the AICC2012 chronology, a relative dating uncertainty of 2 kyr should also be considered for the marine sediment records.As a consequence, caution is still needed when discussing leads and lags of less than about 3 thousand years between model and data.We also use the pollen record from the marine core MD04-2845 in the Bay of Biscay, which provides an estimate of the mean temperature of the warmest month during the LIG (Sánchez Goñi et al., 2012).This record has not been transferred onto AICC2012 but age control points for establishing the chronology of the core are obtained from correlation with well-dated speleothem records (Drysdale et al., 2007) and methane concentration records in ices cores (Waelbroeck et al., 2008).It results in an absolute uncertainty of 2 kyr at the maximum.

Comparison with the proxy data
The temperature increase over Greenland (NEEM site; Fig. 7g) during the early LIG takes place earlier (∼ 5 kyr) in the simulation than in the reconstructions and the maximum of annual mean surface temperature is reached earlier (∼ 2 kyr) in the simulation than in the data.However, the slow cooling after the peak of the LIG is in good agreement between both the simulation and reconstructions from 126 kyr BP until about 119 kyr BP.Although the maximum of summer temperature is slightly delayed compared to the annual mean value, the general behaviour of the evolution of the simulated summer surface temperature at NEEM is not much different from the annual evolution and cannot explain the difference with the reconstructed values.In Antarctica, the model indicates a warming in annual mean of less than 2 • C for the Vostok and Dome C sites from 135 kyr BP to the peak of the LIG, while the data suggest up to 10 • C at Dome C and 7 • C at Vostok.The maximum of annual surface temperature simulated at Dome C site at ∼ 128 kyr BP (Fig. 7h), virtually simultaneous with the reconstructed one, is followed by an almost monotonous decrease, but the simulated cooling is much smaller than the reconstructed one.In other words, the magnitude of the changes in temperature in Antarctica is much smaller in the simulation than in the reconstructions.Moreover, the reconstructed temperatures are lower than the simulated ones.
The pollen data from MD04-2845 indicates a summer warming of up to 7 • C in western France from 135 kyr BP to the peak of the LIG, but it is less than 5 • C in the model (Fig. 7d).The short and abrupt cooling event at ∼ 132 kyr BP documented in the pollen data is also simulated.However, the simulated climate optimum occurs later in the model than in the data.Moreover, western France is warmer in the simulation than in the reconstruction during the warmest period of the LIG.
The magnitude of summer SST change during the LIG is smaller in the model than in the proxy data reconstructions for several regions around the globe (Fig. 7a-c; Capron et al., 2014;Lunt et al., 2013).This is particularly the case in the North Atlantic Ocean where the difference between model and reconstructions can reach several degrees.This is consistent with Langebroek and Nisancioglu (2014) who also show too small amplitudes in the simulated North Atlantic SSTs compared to the reconstructions.Many sites in the North Atlantic experience a cooling event of a few degree magnitude in the model during the termination.This feature is also identified in the reconstructions.The climate optimum in the North Atlantic is simulated earlier than the reconstructed one.The summer SST difference between model and reconstructions is smaller for many regions in the Southern Hemisphere than in the North Atlantic Ocean (not shown).
Several authors already pointed out that the amplitude in the modelled temperature change is smaller than in the proxies during the LIG (see for example Lunt et al., 2013).These authors suggested that some proxies might be biased towards warm growth-season changes.Recently, Bakker and Renssen (2014) confirmed that this effect could at least partly explain the difference in magnitude between simulated and reconstructed temperatures.On the other hand, several processes taking place in the climate system might not be fully represented in the model and might be responsible for at least part of the discrepancies.We can mention, amongst others, the representation of changes of the ocean circulation in the Southern Ocean, changes in the oceanic stratification, changes in the meridional heat exchange in the ocean, as well as the representation of clouds and radiative budget in the atmosphere.Unfortunately, it is not possible to test these hypotheses in the present framework.

Timing of maximum surface temperature
The NH and SH mean July surface temperatures reach a maximum at ∼ 128 kyr BP, simultaneously with the maximum of NH/SH June insolation.The NH and SH mean January surface temperatures exhibit only a very moderate increasing trend during the LIG, with a maximum at ∼ 119 kyr BP, in line with NH/SH October/November insolation.The timing of maximum surface temperature (defined as MWT for timing of maximum warmth; Bakker et al., 2013) at each grid cell of the model is computed here as the time of the maximum surface temperature (ts_max) in the temperature time series smoothed with a Gaussian kernel regression using a bandwidth value of 100.The significance of this result is then tested as follows: the variability of the surface temperature smoothed with bandwidth 10 (ts_10) is compared to the one with a bandwidth 100 (ts_100).Temperatures higher than ts_max minus 0.67 times the standard deviation of (ts_100-ts_10) are considered not significantly different from ts_max.If the time interval, including all these temperatures, lies within ± 0.5 kyr of the age of the temperature maximum, then MWT is considered as significant.
The annual MWT (Fig. 8a) occurs almost all over the globe at 128 kyr BP, i.e. simultaneously with the maximum of NH summer insolation, except over the tropical oceans, where it is observed much later.In particular, the MWT over the east tropical Pacific Ocean is simulated at the end of the LIG.The Sahara desert and the desert regions of the Middle East show their maximum of surface temperature late during the LIG, which is in line with the maximum grass fraction at ∼ 122 kyr BP, quickly followed by a rapid degradation towards desert conditions (strong reduction of grass fraction) in response to a reduction in precipitation (Fig. 10).Some locations over the tropical South Atlantic and the subtropical Pacific Oceans display a very early MWT related to an overshoot of temperature linked to the almost collapse of the AMOC (Fig. 7e).Indeed, this induces a short-term increase of SST in the SH during the reduction of the AMOC, which superimposes on the warming trend of the LIG.The July MWT (Fig. 8b) mirrors the annual MWT, except for the tropical oceans, where MWT occurs before 125 kyr BP.The January MWT (Fig. 8c) is characterised by a late occurrence everywhere, except in the Southern Ocean.These results are in relatively good agreement with a similar study (Bakker et al., 2013).The use of a slightly different methodology and of a different model version may explain most of the differences between both studies.Langebroek and Nisancioglu (2014) suggested that the early occurrence of the MWT in the Southern Ocean is due to the integrating and damping effect of the ocean.The January MWT occurs at the beginning of the simulation over the West Antarctic Ice Sheet and at the end of the simulation over the East Antarctic Ice Sheet, which is coherent with Langebroek and Nisancioglu (2014).

Comparison with other models
The different models included in Bakker et al. (2013) model intercomparison (Bern3D, CLIMBER-2, FAMOUS, LOVE-CLIM (version 1.2, which is similar to IGonly in this study, forced only with variations in insolation and greenhouse gases, see Table 1 and Sect.4.3), CCSM3, KCM, MPI-UW) show a similar trend, magnitudes of change and periods of maximum warmth of the LIG January and July surface temperature evolutions for most latitudinal bands.In particular the maximum warmth in January occurs towards the end of the simulation (later than 120 kyr BP) and July surface temperature is maximum in the early part of the simulation (before 128 kyr BP), for most of the models and latitudes.This is in agreement with the allGR simulation.Indeed, although the melting NH ice sheets have a strong impact on the allGR simulated climate, the maximum of warmth, strongly related to the orbital forcing and occurring later than 130 kyr BP, is hardly affected by the melting of the ice sheets.
However, there is a disagreement between the models as to the sign of the change for some latitudes and periods of the year.For example, LOVECLIM1.2simulates surface temperatures higher than the pre-industrial ones at least for a few  , 2011).The simulated series are smoothed using a moving average over 100 yr and averaged over four adjacent grid cells.The black dot on the right hand side of the figures provides the corresponding simulated pre-industrial value.The "plateau" is shaded in grey.Note that a coherent temporal framework has been recently constructed for the ENAM33, ODP980, NA87-25, SU90-08, SU90-03, V30-97, NEEM and EDC records and they are all displayed here on the recent AICC2012 chronology (Capron et al., 2014;Bazin et al., 2013).The location of the marine and ice cores is provided in Fig. 6 and in Table 2. thousand years during the LIG, in agreement with allGR in this study.But some models, such as CCSM3 in the high northern latitudes in January and in the high southern latitudes in January and July, simulate LIG surface temperature lower than the pre-industrial ones.
LOVECLIM1.2 is not the only model to show significant changes in the strength of the AMOC, without rapid change in external forcings or surface boundary conditions.FAMOUS simulates a switch around 121 kyr BP from a weak AMOC with the main site of deep convection located in the North Pacific to a strong AMOC with deep convection taking place mainly in the North Atlantic.The Bern3D model, which incorporates prescribed remnant ice sheets in the NH and the related meltwater flux to the ocean, simulates the LIG weakest AMOC between 129 and 125 kyr BP, when NH ice sheets are melting.The AMOC then remains stable before strengthening from 121 kyr BP.This is in line with the evolution of the AMOC during the LIG as simulated with allGR, except that AMOC in allGR tends to weaken after 117 kyr BP.Moreover, the weakening of the AMOC in allGR at 132 kyr BP is much larger than with Bern3D.The LIG temperature evolution simulated with allGR is strongly affected by the changes in the AMOC strength, contrarily to the Bern3D simulation.CLIMBER-2 also shows concurrent rapid changes in the LIG surface temperature of several de-grees and small amplitude changes in the AMOC strength.However, they are of opposite sign.
In the next sections, we assess the role of different surface boundary conditions described in Sect.2.3 and of the different parameter sets described in Sect.2.4.

The role of the ice sheet configuration (elevation, extent and albedo)
First, we compare the reference simulation (allGR) with a simulation that does not take into account the evolution of the NH ice sheet configuration but only includes the freshwater forcing resulting from changes in ice volume (fwfGR).
The global annual mean surface temperature (Fig. 4a) at the beginning of the simulation, until ∼ 131 kyr BP, is up to 0.9 • C (at 133.4 kyr BP) higher in fwfGR than in allGR and the cooling at the end of the simulation, from 119 kyr BP, is slightly smaller in fwfGR than in allGR (0.5 • C at 115.8 kyr BP).Many regions in the NH, both over land and over the ocean, as well as over Greenland, also show this characteristic.The major features of annual MWT are similar in fwfGR and allGR suggesting that the ice sheet configuration does not significantly impact MWT (Figs. 8a and 11b).Except during the plateau, the strength of the AMOC is most of the time slightly smaller in fwfGR than in allGR, although its variability is much larger in fwfGR than in allGR.This suggests a stabilisation effect of the NH ice sheets configuration on AMOC, or rather that fwfGR is closer to a bistable regime than allGR.Indeed, when the NH ice sheets are included (as in allGR), the sea surface temperature in their vicinity decreases (Fig. 9).This favours deepwater formation (Renssen et al., 2005) and may therefore stabilise the overall Atlantic overturning circulation although the link between the mean state and the variability of the AMOC is not necessarily straightforward.During the plateau, between 131 and 121 kyr BP, both global annual mean surface temperature and AMOC are similar in both simulations.
The conclusions drawn from the comparison of temperature between allGR and the reconstructions are not strongly modified if the changes in the configuration of the NH ice sheets during the LIG are not taken into account.The magnitude of the variations through the LIG remains too low for most of the variables.In both allGR and fwfGR, the model does not reproduce the large magnitude of the reconstructed SST variations, mostly in the NH.Moreover, the large increase in global annual mean surface temperature is simultaneous in both simulations and occurs too early compared to the reconstructions.

The role of additional freshwater flux
The reference simulation (allGR) is compared here to a simulation that does not take into account the additional freshwater flux from the ice sheets but only includes the evolution of the NH ice sheet configuration (extent, altitude, albedo) (to-poGR).Without this input in the North Atlantic, the AMOC remains more or less in its initial state until about 132 kyr BP and does not experience any collapse (Fig. 4b), although there is a rapid decrease of 3 Sv in 150 yr at 132 kyr BP.Then, its variability becomes slightly higher (the standard deviation is 1.1 Sv at the beginning of the simulation, while it rises to 1.5Sv later) and its strength becomes similar to the one simulated in allGR (∼ 19 Sv).Simultaneously, the global annual mean surface temperature (Fig. 4a) smoothly increases to reach a value similar to the one simulated in allGR at 131 kyr BP.Thus the additional freshwater flux is clearly responsible for the dramatic weakening of the AMOC of almost 18 Sv at 132 kyr BP (resulting in a collapse of the AMOC) and the cooling of 1.3 • C at the same time simulated in allGR.A decrease of several degrees in the simulated North Atlantic SST and western France and Greenland surface temperatures may also be attributed to the additional freshwater fluxes.During the plateau, between 131 kyr BP and 121 kyr BP, the differences between the simulations (to-poGR and allGR) are only minor.After 121 kyr BP, at the end of the simulation, when the NH ice sheets start to regrow, the AMOC variability is smaller in the simulation without freshwater flux (topoGR) than with the one that includes it (allGR).
The general pattern of annual MWT is very similar with (Fig. 8a) and without (Fig. 11a) additional freshwater flux due to ice sheet melting, with many regions of the globe experiencing an early MWT (∼ 128 kyr BP).MWT occurs later than 123 kyr BP in similar regions in both allGR and to-poGR.While allGR shows a rather uniform MWT for southern and eastern Asia (∼ 119 kyr BP), topoGR displays a late MWT for southern Asia (115 kyr BP) and an early MWT for eastern Asia (122 kyr BP).July MWT is mostly characterised with a later occurrence over the southern Atlantic Ocean in topoGR compared to allGR.January MWT occurred about 1 kyr later for many places in topoGR compared to allGR.However, it is characterised with an earlier occurrence (earlier than 123 kyr BP) for topoGR than for allGR in the North Atlantic, in particular over the Greenland Sea.

The orbital and greenhouse gas forcings alone
Lastly, the simulations taking into account either all (allGR) or none (IGonly) of the changes related to the growth and decay of the NH ice sheets are compared.
The global annual mean surface temperature is at maximum at ∼ 128.4 kyr BP in IGonly as in allGR (Fig. 4).Although the warming from the previous glacial maximum is fast, about 0.21 • C kyr −1 from 135 kyr BP, there is no rapid increase as simulated in allGR at about 132 kyr BP.Even if the initial conditions of allGR and IGonly at 135 kyr BP and their general behaviour until 131 kyr BP are strongly different, the simulations converge after 131 kyr BP.The AMOC intensity reaches a minimum at ∼ 128.7 kyr BP, almost simultaneous to the maximum of global annual mean surface temperature.However, it shows hardly any trend or variability before 121 kyr BP in IGonly.In particular, the   AMOC strength and temperature display a much larger variability and smaller values than before the event.
The annual MWT is very similar in IGonly (Fig. 11c) and allGR (Fig. 8a), although it occurs later in IGonly than allGR at some locations in the South Pacific and South Atlantic, which may be related to the response of the ocean in the SH to the lack of input of freshwater in the North Atlantic.In January, MWT takes place a few thousand years earlier over the North Atlantic in IGonly compared to allGR, and much earlier over the northeastern Pacific, along the coast of North America.
We can conclude from this that the evolution of the NH ice sheets greatly increases the magnitude of the temperature variations over continents as well as over the ocean at the beginning of the simulation and improves the agreement between the simulated climate and the reconstructions.

Separation factor technique -the importance of the synergism
In the previous sections the simulations were discussed within the framework of a classical methodology for sensitivity analysis, based on a single perturbation at a time.However, the synergism between the different factors is not quantified.The synergism between two factors corresponds to their joint effect that is reflected in either the enhancement or the reduction of their cumulated impact.The separation factor technique, first designed by Stein and Alpert (1993) and further extended later for palaeo-climate studies (Claussen et al., 2001;Berger, 2001;Crucifix and Loutre, 2002), allows to quantify the direct effects of the factors as well as their synergism.Previous studies applied the factor separation technique to components of the climate system or to climate forcings (Yin and Berger, 2012).Here, we apply it to investigate the role of the boundary conditions.The pure contribu-tion of the NH ice sheet configuration is given as (topoGR -IGonly).The pure contribution of the additional freshwater flux from the melting ice sheets is given as (fwfGR, IGonly).The contribution of the interaction (synergism) between these two factors is given as (allG, fwfGR, topoGR + IGonly).This analysis (Fig. 12a) confirms that, in the model, the NH ice sheet configuration is responsible for a decrease in global annual mean surface temperature of 0.9 • C from 133.5 to 131.3 kyr BP, while the additional freshwater flux is responsible for a cooling of −0.2 ± 0.1 • C over the first 2.5 kyr of the simulation, followed by a short duration (1250 yr) cooling event of up to 1.3 • C.Moreover, the synergism between the two factors, in other word the joint effect of the NH ice sheet configuration and of the freshwater forcing resulting from changes in ice volume, is responsible for a positive contribution to the global annual mean surface temperature at 135 kyr BP that becomes negative (0.1 • C) at 131.3 kyr BP.Between 131 and 121 kyr BP, during the plateau, the insolation and greenhouse gas forcings are responsible for all the changes in global annual mean surface temperature.After 121 kyr BP, at the end of the simulation, the synergism has a negative contribution (less than 0.4 • C) until 117 kyr BP.This contribution then oscillates between positive and negative values.
At the beginning of the simulation, before 131 kyr BP, the NH ice sheet configuration does not contribute to changes in the AMOC (Fig. 12b), while the largest (negative) contribution comes from the additional freshwater flux (up to 18 Sv).The synergism has a small positive contribution, leading to a strengthening of about 3 Sv, with a large variability.Similar to temperature, the insolation and greenhouse gas forcings are responsible for the whole change in AMOC, during the plateau, between 131 and 121 kyr BP.After 121 kyr BP, at the end of the simulation, the contribution of the freshwater flux and the NH ice sheet configuration are both positive (strengthening of the AMOC of up to, respectively 13 and 7 Sv), while it is negative for the synergism (reduction of the strength of the AMOC of up to 7 Sv).

Sensitivity to the reconstructed evolution of the ice sheets
Reconstructions of sea level change during the LIG (see Sect. 1) strongly suggest that ice sheets were still present over North America and Eurasia at 135 kyr BP.However, their precise extent and the exact timing of their melting are not known accurately.The impact on the simulated climate of a later NH ice sheet melting and glacial inception is analysed here through the use of two NH ice sheet reconstructions (GR and LR, see Sect.2.3) for the LIG simulation (allGR and al-lLR, respectively, see Table 1).Simulations allGR and allLR differ only from the NH ice sheet configuration and freshwater flux.Grant et al. (2012) already suggested that deep ocean temperature bias and "incorrect" orbital tuning could explain the difference in timing between both scenarios.The impact of this difference on the timing of the maximum of temperature (MWT) is strong.The annual MWT occurs about 5 kyr later in allLR (Fig. 14b) than in allGR (Fig. 8a) in southern and eastern Asia.In the southern Pacific Ocean and in the southern Indian Ocean, MWT occurs around 130 kyr BP in allLR, while it takes place between 132 and 125 kyr BP in allGR.July MWT occurs 1-2 kyr later in allLR than in allGR in the North Atlantic Ocean and southern Africa, while it is about 5 kyr earlier in allLR than in allGR in the central Pacific.
Most of the impact of the choice of the scenario for the NH ice sheets evolution during the LIG is concentrated on the beginning of the simulation.The moderate continuous input of freshwater in the LR scenario leads to a colder Greenland, western France and North Atlantic in allLR than in allGR (Fig. 13c-e); the SH is less affected (not shown).However, no strong short-lasting cooling event is simulated in allLR as is the case in allGR at about 132 kyr BP.Moreover, the freshwater flux into the North Atlantic at 128 kyr BP in allLR has a cooling effect of ∼ 2 • C over western France (0.7 • C over the NH) at a time when temperature reaches its maximum in allGR.Therefore, the maximum of surface temperature over western France is delayed with allLR compared to allGR.The warming of the North Atlantic is slower in allLR than in allGR (Fig. 13c) and therefore takes longer, from 133 to 127 kyr BP.Thus, the agreement with the reconstructions is better, although the differences are still large.Over Greenland the agreement between reconstructed and simulated (allLR) before 122 kyr BP is very good.Thereafter, the cooling proceeds differently in observations and simulations, although the general trends remain in good agreement (Fig. 13e).
Compared to the GR scenario, the LR scenario induces a delayed warming, mostly over Europe and the North Atlantic, at the beginning of the simulation.This delay induces large differences in surface temperature between the simulations using either GR or LR scenario.The difference reaches almost 10 • C over Greenland, 5 • C locally over the North Atlantic and almost 4 • C over western France.The comparison with proxy data shows that the LR scenario for the evolution of the ice sheets leads to a better agreement between modelled and reconstructed climates than the GR scenario.

The role of the parameter sets
In order to test the potential impact of the model parameter sets, we performed and compared simulations with either parameter set std or parameter set 22, which has a higher sensitivity to a freshwater flux in the North Atlantic than std.Moreover, we used both scenarios for the LIG NH ice sheet (see Table 1) configuration and freshwater flux, either based on Grant et al. (2012) or on Lisiecki and Raymo (2005) (see Table 1).In all these simulations (allGR, allLR, parGR, parLR), the model is forced with changes in insolation and atmospheric greenhouse gases.
The AMOC in parGR is weaker than in allGR throughout LIG.Moreover, the difference between the simulations decreases from 4.2 Sv at 135 kyr BP to 1.4 Sv at 115 kyr BP (Fig. 13b).Prior to 127 kyr BP, the AMOC is strongly affected by the evolution of the ice sheets.Nevertheless, for the same ice sheet scenario, the AMOC simulated with parameter set 22 is slightly weaker than with std, and remains weaker until 121 kyr BP.Starting from 120 kyr BP, the role of the parameter set in the simulated AMOC is small compared to the role of the ice sheets.
The global annual mean surface temperature is cooler in parGR than in allGR during the LIG.The average difference is 0.3 • C. It amounts to ∼ 0.6 • C at 135 kyr BP, decreasing to only ∼ 0.2 • C between 127 and 121 kyr BP, then increasing back to 0.5 • C at 115 kyr BP.Between 127.0 and 120.5 kyr BP, the temperature difference related to the choice of the parameter set may vary according to the location.For most of the locations, the major driver of temperature changes is the ice sheet changes, in particular the freshwater flow from the ice sheets.Similar conclusions can be drawn from the comparison of allLR and parLR.Therefore, we can conclude that the contribution of the parameter sets is actually very modest.However, we must take care here because the sensitivity of model configurations used is in the lower end of the full range of the tested model configurations (Goosse et al., 2007;Loutre et al., 2011;Goelzer et al., 2011).(Capron et al., 2014;Bazin et al., 2013).The black dot on the right hand side of the figures provides the corresponding simulated pre-industrial value.The "plateau" is shaded in grey.The location of the marine and ice cores is provided in Fig. 6 and in Table 2.It is also characterised by a decrease in surface temperature in the Labrador Sea and in the Barents Sea, a decrease in sea surface salinity in the Hudson Bay, Baffin Bay and Davis Strait, a change in the pattern of convection in the Labrador Sea, and an increase in winter sea ice area in the NH.

Clim
For the Holocene, palaeoceanographic data suggest that the Labrador Sea may be close to a bistable regime and could therefore potentially induce large oscillations in the climate system (Hillaire-Marcel et al., 2001).Furthermore, several model studies suggested that, during the Holocene, small freshwater forcing in the Labrador Sea may trigger the climate system into a bistable regime (Kuhlbrodt et al., 2001;Wood et al., 1999;Jongma et al., 2007).Under these conditions, small climatic perturbations may be amplified into large climatic changes.Similarly, the climate system may also be close to a bistable regime during the LIG.Indeed, abrupt climate changes are recorded in high-resolution marine proxies from the low-latitude east Atlantic margin, related to freshening and cooling in the Norwegian Sea (Maslin et al., 1998).It was suggested that they are closely connected to AMOC reorganisation (Ganopolski and Rahmstorf, 2001;Alley et al., 2001;Timmermann et al., 2003).Recently, Galaasen et al. (2014) identified episodes of brief and intermittent reduction in NADW flow lasting a few centuries between 115 and 130 kyr BP, which they associated with changes in the Nordic Seas overflow.However, the uncertainty in both the simulation and reconstructions prevents drawing any direct link between them.
Transient simulations of the LIG performed with climate models other than LOVECLIM also display such large variability of the AMOC (Bakker et al., 2013).For example, the AMOC change in FAMOUS has an opposite sign compared to IGonly and topoGR simulations and CLIMBER-2 simulated several high amplitude variations of the AMOC.Friedrich et al. (2009), using LOVECLIM, suggested that such rapid changes may be due to a flush of freshwater from the Hudson Bay to the Labrador sea due to change in wind.Therefore, they are not related to any additional freshwater flux but rather to internal variability.However, sensitivity tests performed in the framework of this work to test this hypothesis show that the closure of Hudson Bay prevents neither the occurrence of the rapid event nor the general pattern of changes in surface temperature, sea surface temperature, convection depth, and sea ice extent.In any case, as the occurrence of such events is strongly dependent on the model and of the selected forcings, they should be interpreted with care.

The role of the freshwater flux from Antarctica
The simulations performed in this work did not include the freshwater fluxes from the Antarctic ice sheet.However, Mathiot et al. (2013), using LOVECLIM with data assimilation to study the mechanisms responsible for southern high latitude cooling during the Holocene from 10 to 8 kyr BP, concluded that the Southern Ocean cooling is mainly driven by an increase in freshwater flux from Antarctica, while the continent is only slightly affected.As we may expect a similar impact during the LIG, a further step would be to quantify the evolution of the Antarctic ice sheet during the LIG (extent, altitude, amount of freshwater from its melting if any) to estimate its impact on climate during the LIG.

Uncertainties associated with the dating of paleoclimatic records
Several chronologies are involved in our comparisons between model and data.The timing of our prescribed freshwater flux is linked to the timing of the sea level curve, either from Grant et al. (2012) or Lisiecki and Raymo (2005), with the former leading the latter by about 2.8 kyr.We showed that choosing one or another may have a large impact on the timing (several thousand years) of the simulated climate (e.g. the timing of large temperature changes in the NH at the beginning of the simulation).For the model-data comparison, all paleoclimatic records have been synchronised onto the AICC2012 chronology associated with a 2-thousand-year absolute error, except core MD04-2845 (Capron et al., 2014;Sánchez Goñi et al., 2012).The synchronised marine records are associated with a total age uncertainty of about 3 thousand years.The absolute dating error for core MD04-2845 is less than 2 thousand years.Caution is thus needed when discussing leads and lags of less than about 3 thousand years between model and data.The ice below 2206.7 m from the NEEM ice core, corresponding to 108 kyr BP, is disturbed and folded (NEEM community members, 2013).Relevant information is extracted from comparison between NEEM and other ice cores for several variables.This process may also lead to uncertainties.Although some differences between climate models and reconstructions could be reduced when taking into account these uncertainties, there remain significant differences between them, in particular when amplitude is concerned.

Estimation of the contribution of the ice sheets
The evolutions of the ice sheets used here have significant uncertainties.The scenarios could certainly be improved and more precise ice sheet reconstructions for the LIG will become available in the future.However, although our scenarios of changes in the NH ice sheet configuration and FW flux have an impact on the simulated climate (at least for some climate variables at some locations), it is not the only factor explaining model-data disagreements.For instance, the changes in the NH ice sheet configuration induce an early occurrence (before 129 kyr BP) of the annual MWT in the SH, in the eastern Atlantic Ocean and in the western Pacific, in the New Zealand region, but a late (after 120 kyr BP) annual MWT in the Greenland Sea as well as a late MWT almost everywhere in the NH in January.However, the impact of the NH ice sheet configuration changes on the simulated timing of the LIG remains rather small and cannot fully explain model-data disagreements.
The contribution of the ice sheets to the change in global annual mean surface temperature during the LIG before 131 kyr BP reaches 1 • C related to changes in the ice sheet configuration only and up to 1.3 • C when the melting ice sheets release a large amount of freshwater into the ocean.The timing of the temperature rise at the beginning of the interglacial as well as its rate strongly depend on the evolution of the NH ice sheets.However, the presence of remnant ice sheets at the beginning of the LIG does not have long-term impact on the climate.Indeed, it has almost no influence between 131 and 121 kyr BP when only the Greenland ice sheet remained in the NH.After 120 kyr BP, the contribution of the different components of changes in the NH ice sheets shows a large variability (±0.5 • C).The changes in the AMOC are large (up to 18 Sv) only when the melting ice sheets deliver additional freshwater to the ocean.
The comparison between allGR and allLR suggests that the uncertainty on the surface boundary conditions can lead to a difference in the simulated global annual mean surface temperature of up to 1.4 • C before 131 kyr BP, but smaller than 0.8 • C at the end of the simulation (115 kyr BP).The maximum change in the AMOC induced by the changes in the boundary conditions varies from up to 18Sv (at the beginning of the simulation) to up to 13Sv (at the end of the simulation).This uncertainty in magnitude can be viewed, at least partly, as an uncertainty on the chronology related to the timing of the freshwater flux from the melting ice sheets.
NH ice sheets are quite instrumental in the evolution of the temperature before 132 and after 119 kyr.By contrast, the use of different parameter sets induces only small changes in the simulated climate.

Summary and conclusions
In this study we used an Earth system model of intermediate complexity to investigate the importance of the surface boundary conditions, in particular the changes during the LIG of the ice sheets configuration (extent, altitude, albedo) and of the additional freshwater flux from the resulting ice The impact of some changes in model parameters was also considered.The uncertainty is expressed here as the impact of the boundary conditions and parameter sets on a climate variable (mostly temperature and strength of the AMOC).The simulations can be broadly divided into three subintervals: the beginning of the simulation (135-131 kyr BP) potentially impacted by the remnant ice sheets, the subinterval 131-121 kyr BP when the Greenland ice sheet is the only ice sheet present in the NH (the plateau), and the end of the simulation (121-115 kyr BP) characterised by the glacial inception.At the beginning of the simulation, uncertainties in the global annual mean surface temperature can be attributed to the ice sheet configuration (0.9 • C), to additional freshwater flux (1.3 • C), to the choice of the scenario (GR versus LR) for the NH ice sheet evolution (1.5 • C) and to the choice of the parameter set (0.6 • C).Moreover, the choice of the scenario for the NH ice sheet evolution is mostly responsible for an uncertainty of several thousand years in the time response of the climate.The uncertainty remains small (less than 0.2 • C) between 131 and 121 kyr BP.At the end of the simulation, the ice sheet configuration and the choice of the parameter set are each responsible for an uncertainty of about 0.5 • C in the global annual mean surface temperature.It must be noted that these uncertainties are not additive.Moreover, an abrupt event taking place in a few centuries, such as the event discussed in Sect.7.1 is responsible for a further uncertainty of 0.4 • C.
The magnitude of summer SST change during the LIG is smaller in the model than in the proxy data reconstructions for several regions around the globe (Capron et al., 2014;Lunt et al., 2013).This is particularly the case in the North Atlantic Ocean where the difference between model and reconstructions can reach several degrees although taking into account the evolution of the NH ice sheets reduces the discrepancy.The July MWT occurs almost all over the globe at 128 kyr BP and the January MWT is characterised by a late occurrence everywhere, except in the Southern Ocean.This is in disagreement with the data showing a late warming in the NH compared to the SH.Orbital and GHG forcings are responsible for most of the climate changes between 131 and 121 kyr BP.Moreover, the evolution of the ice sheets prior 131 kyr BP has a negligible impact on the simulated climate between 131 and 121 kyr BP.In other words, there is no strong climate memory at that time.
The ice sheet configuration has only a small contribution to changes in the strength of the AMOC although it can be responsible for a high variability.However, the additional freshwater input may induce a large change in the AMOC and therefore may be responsible for an uncertainty of several Sv (up to 18 Sv in our simulations).Moreover, the choice of the scenario leads to a further uncertainty on the timing of the event.A rapid event, related to change of 4 Sv also occurred between 120 and 115 kyr BP in most of the simulations.It is thus suggested that such an event may be respon-sible for an uncertainty of several Sverdrups.In addition, the choice of the parameter set is responsible for a difference in the strength of the simulated AMOC between 1.4 and 4.2 Sv.
The evolution of the NH ice sheets greatly increases the amplitude of the climate changes before 131 kyr BP.The additional freshwater flux (FWF) from the melting NH ice sheets is responsible for the major contribution, while the changes in the configuration (extent and albedo) of the NH ice sheets only slightly impact the simulated climate.The evolution of the ice sheet is critical in the timing of the warming at the beginning of the LIG, mostly over Europe and the North Atlantic.The uncertainty in the MWT can reach several thousand years depending on the scenario of evolution of the ice sheets.Therefore, the scenario of the NH ice sheet evolution, in particular its timing, is essential for an accurate simulation of the surface temperature, in particular in the North Atlantic region.Further modifications of the climate response can be expected when models of the Greenland and Antarctic ice sheets are interactively coupled within LOVECLIM as has been done for future projections (e.g.Huybrechts et al., 2011;Goelzer et al., 2012a) and is in preparation for the LIG period in a forthcoming publication.
Clearly other sources of uncertainty, not studied here, such as uncertainties related to choices made in the representation of the model components, to their physical parameterisations and to the model resolution (Murphy et al., 2004;Stainforth et al., 2005), also affect the simulation of the climate evolution during the LIG.All these uncertainties, as well as the uncertainties on the chronology and on the values of the reconstructed variables, may significantly hamper our ability to confidently draw conclusions from model-data comparisons during the LIG.Taking into account the time evolution of the ice sheets during the LIG improves the simulation of the amplitude of the climate response during the LIG.Nevertheless, even if taking into account all the uncertainties, the model fails to reproduce the very large amplitude changes of the climate during the LIG (in particular in the NH surface temperature) found in proxy reconstructions.However, the scenario of ice sheet evolution is critical for the simulation of the time evolution of the climate.The low temperatures at the beginning of the LIG are simulated simultaneously to the large freshwater flux from the ice sheets, while the maximum temperature is mostly driven by the solar forcing.This can cause a phase shift between the simulated and reconstructed climates.
The Supplement related to this article is available online at doi:10.5194/cp-10-1541-2014-supplement.
search project (Constraining long-term climate and sea-level projections using the Last Interglacial) within its Research Programme on Science for a Sustainable Development (SD/CS/06A).
The research leading to these results has received funding from the European Union's Seventh Framework programme (FP7/2007-2013) under grant agreement no.243908, "Past4Future.Climate change -Learning from the past climate".This is Past4Future contribution 65.

Figure 1 .
Figure1.Deviations from the present-day values of the insolation at the boreal summer solstice (left) and at the boreal winter solstice (right) for the LIG(Berger, 1978).

Figure 3 .
Figure 3. Freshwater flux from melting ice sheets entering the ocean for four regions: the Pacific (PAC), Arctic Ocean (ARC), North Atlantic (NAT) and North Sea (NOS).FW flux is estimated based on sea level change reconstructions from (left; a) Grant et al. (2012) and (right; b) Lisiecki and Raymo (2005).

Figure 4 .
Figure 4. Time evolution of (a) (left) global annual mean surface temperature ( • C); (b) (right) maximum of the meridional overturning stream function in the North Atlantic from model simulations using different surface boundary conditions (see text and Table1).The series are smoothed using a moving average over 100 yr.The black dot on the right hand side of the figures provides the corresponding simulated pre-industrial value.

Figure 5 .
Figure 5. Winter (DJF) convection depth (m) simulated with allGR.(left) time evolution for four sites of the North Atlantic: Labrador Sea (black dashed line), Davis Strait (black full line), Norwegian Sea (blue dashed line) and Barents Sea (blue full line).These sites are identified with boxes on the right hand side of the figure.(right) convection depth (m) for the North Atlantic at 134 (top) and 116 kyr BP (bottom).

Figure 8 .
Figure 8. Timing (in kyr) of the simulated LIG surface temperature maximum (MWT) in annual mean (a), and for July (b) and January (c) for allGR.See text for the definition of MWT.

Figure 10 .
Figure 10.(a) Changes in land cover during the LIG over the Sahara region.(b) Time evolution of precipitation in January, July and in annual mean during the LIG over the Sahara region for simulation allGR.

Figure 12 .
Figure 12.Separation factor method applied to the global annual mean surface temperature (a) and the AMOC (b).The pure contribution of the ice sheets configuration (cf1 in blue), of the additional freshwater flux related to the melting of the ice sheets (cf2 in orange) and of the synergism between the two factors (cf12 in brown).
Figure 13.Comparison between the LIG climate simulated using different parameter sets (see text and Table1) and different Northern Hemisphere ice sheet evolutions (see text and Table1), and the reconstructed climate.Time evolution of (a) global annual mean surface temperature ( • C); (b) maximum of the meridional overturning stream function in the North Atlantic from model simulations using different surface boundary conditions; (c) NE Atlantic summer SST ( • C) (ODP 980; 55.8 • N, 14.11 • W, Oppo et al., 2006; NA87-25; 55.57 • N, 14.75 • W, Cortijo et al., 1994); (d) July surface temperature over western France (MD04-2845; 45 • N, 5 • W; Sánchez Goñi et al., 2012); and (e) Precipitation-weighted temperature reconstruction corrected for elevation change at the NEEM site, Greenland (77.45 • N, 51.06 • W)(NEEM community members; 2013).The simulated series are smoothed using a moving average over 100 yr and averaged over four adjacent grid cells.NA87-25, ODP980 and NEEM records are displayed here on the recent AICC2012 chronology(Capron et al., 2014;Bazin et al., 2013).The black dot on the right hand side of the figures provides the corresponding simulated pre-industrial value.The "plateau" is shaded in grey.The location of the marine and ice cores is provided in Fig.6and in Table2.
on the uncertainty in the simulated climate.

Table 1 .
Description of the transient simulations.All the simulations are forced by changes in insolation and atmospheric GHG concentrations (see text).The reference parameter set (std)

is de- fined and used in Goosse et al. (2010), while parameter set 22 is defined and used in Loutre et al. (2011) and Goosse et al. (2007) (see also supplementary material). PI stands for pre-industrial; GR (respectively LR) means that the LIG ice sheet configuration and/or freshwater flux is reconstructed based on Grant et al. (2012) (re- spectively
Lisiecki and Raymo, 2005)).

Table 2 .
Location of the study sites, including marine and ice core sites.Site (12) corresponds to an unpublished marine core for which only the simulated temperature is available.