Last interglacial temperature evolution – a model inter-comparison

. There is a growing number of proxy-based reconstructions detailing the climatic changes that occurred during the last interglacial period (LIG). This period is of special interest, because large parts of the globe were characterized by a warmer-than-present-day


Introduction
To strengthen our confidence in climate models, it is important to assess their ability to realistically simulate a climate different from the present-day climate (Braconnot et al., 2012).The last interglacial period (LIG; ∼ 130 000-115 000 yr BP) provides an interesting period, because many proxy-based reconstructions show temperatures up to several degrees higher than present day (CAPE-members, 2006;Turney and Jones, 2010;McKay et al., 2011).However, to date, the evolution of the climate during the LIG is still under debate.This is especially true for the establishment of peak interglacial warmth in different regions.For instance, proxy-based reconstructions of surface temperatures from the Norwegian Sea and the North Atlantic are inconclusive on whether peak interglacial warmth occurred in the first or in the second part of the LIG (Bauch and Kandiano, 2007;Nieuwenhove et al., 2011;Govin et al., 2012).The main cause of this uncertainty is the difficulty in establishing a coherent stratigraphic framework for the LIG period, not only between different regions (e.g. the Norwegian Sea and the North Atlantic) but also between different types of proxyarchives (e.g.speleothems, ice cores, deep-sea cores and lake sediments, e.g.Waelbroeck et al., 2008).Deciphering the evolution of LIG surface temperatures is further complicated by the fact that different types of proxies record different parts of the climatic signal: for instance maximum summer warmth, the number of days above a threshold temperature, the seasonal temperature contrast or average summer temperatures (Jones and Mann, 2002;Sirocko et al., 2006).Climate simulations covering the LIG period can be used to facilitate the interpretation of proxy-based temperature reconstructions by providing information on the timing of peak interglacial warmth for different months and on possible spatial differences in the evolution of temperatures.
For the LIG period a large number of equilibrium simulations have been analysed (Montoya, 2007;Lunt et al., 2012, and references therein).However, to investigate the evolution of temperatures throughout this period and the timing of maximum warmth (MWT), the transient nature of two of the major forcings, changes in the astronomical configuration and changes in the concentrations of the major greenhouse-gases (GHGs), has to be incorporated.A small number of transient climate simulations have previously been performed for the LIG (e.g.Calov et al., 2005;Gröger et al., 2007;Ritz et al., 2011a).Amongst other things they indicate the importance of changes in the overturning circulation and the sea-ice cover.These are however known to be highly model-dependent stressing the need for a larger model intercomparison.
In this study we present the first investigation of the common LIG temperature signal in long, > 10 000 yr, transient simulations performed with seven different climate models.Included are both published LIG transient simulations (Gröger et al., 2007) as well as ones recently performed within the PMIP3 framework (Paleoclimate Modelling Intercomparison Project).The climate models used in this inter-comparison study differ in complexity from 2.5-D atmosphere-ocean-vegetation models to general circulation models (GCMs).Some also differ in terms of the climatic forcing and in the components of the climate system which are included (Table 1).
The objectives of this model inter-comparison are the following: (1) to determine the transient temperature response to LIG forcings which is common to the different models, (2) to analyse the simulated spatio-temporal response of temperatures during the LIG and (3) to indicate in which regions climatic feedbacks likely played a crucial role in shaping the LIG temperature evolution.This study provides an important step towards a future comparison of LIG proxy-based reconstructions and transient model simulations.

Model simulations
We performed transient LIG climate simulations with a total of seven different climate models of different complexity.In the following model descriptions, we focus only on the relevant differences between the models and the simulation design.For a complete description, see Table 1.In the second part of this section, an overview of the evolution of the main climate forcings of the LIG period is given in terms of changes in the insolation received by Earth and the changes in the GHG concentrations.

Bern3D
The Bern3D Earth system model of intermediate complexity (EMIC) consists of a two-dimensional atmospheric energy and moisture balance model that is coupled to a threedimensional sea-ice-ocean model.In the atmospheric component, heat is transported horizontally by diffusion only while moisture is transported by both diffusion and prescribed advection (Edwards and Marsh, 2005;Müller et al., 2006;Ritz et al., 2011a,b).This means that, compared to other models, the spatial and temporal changes in surface temperatures simulated by the Bern3D model are more directly linked to local changes in the radiative forcing.This simulation includes prescribed changes in the extent of the Northern Hemisphere (NH) continental ice sheets (the Antarctic ice sheet is fixed to present-day configuration because of the coarse resolution of the model at high latitudes).The extent of the NH continental ice sheets is calculated using the benthic δ 18 O stack (a proxy for global ice volume) of Lisiecki and Raymo (2005) in order to scale the ice sheet between the modern and the Last Glacial Maximum extent (Ritz et al., 2011a).Consequently, until ∼ 125 ka BP remnants of the North American and Eurasian ice sheets from the preceding glacial period are prescribed next to the Greenland Table 1.List of the main features of the climate models involved in this model inter-comparison.Unless stated otherwise, sea-level, vegetation cover and ice-sheet configuration are fixed to pre-industrial values.The equilibrium GHG values used by CCSM3 and KCM are respectively 272 and 280 ppm for CO 2 , 622 and 760 ppb for CH 4 and 259 and 270 ppb for N 2 O.In column six information about the applied spin-up procedure is given by stating if the spin-up was an equilibrium (eq.) or transient (trans.)simulation, what the corresponding start year (ka BP) of the corresponding forcing scenario is and, in the case of an equilibrium spin-up simulation, what the length (ka) of the spin-up period is (between brackets).The other used acronyms are Earth system model of intermediate complexity (EMIC), general circulation model (GCM), astronomical configuration (orb), astronomical acceleration with a factor of 10 (acc), greenhouse-gas concentrations (ghg), and prescribed changes in ice sheet configuration (ice).

CCSM3
The CCSM3 (Community Climate System Model, version 3) is a GCM which is composed of four components representing atmosphere (CAM3), ocean (POP), land, and sea ice (Collins et al., 2006).For the simulation in this study, the low-resolution version of CCSM3 is used, which is described in detail by Yeager et al. (2006).The transient simulation has been carried out with 10 times accelerated astronomical forcing (see Sect. 2.3 for details).

CLIMBER-2
The CLIMBER-2 EMIC is a 2.5-D atmosphere-oceanvegetation model of intermediate complexity (Petoukhov et al., 2000).The atmospheric component is a low-resolution 2.5-D statistical-dynamic model.The oceanic component is a zonally averaged multi-basin (Atlantic, Indian and Pacific) model which resolves these basins only in the latitudinal direction.CLIMBER-2 includes a thermodynamic seaice model that computes the evolution of sea-ice coverage and thickness.Note that, in contrast to most of the simulations in this inter-comparison, vegetation is actively simulated.This transient simulation is part of a longer simulation covering the last 4 glacial-interglacial cycles (420-0 ka BP) resulting in the initial conditions of this simulation being different from the other simulations.

FAMOUS
The FAMOUS GCM (Jones et al., 2005;Smith and Gregory, 2012;Smith, 2012) is a low-resolution version of the HadCM3 GCM (Gordon et al., 2000) with roughly half the horizontal resolution in both the atmosphere and ocean and a longer time step.

Kiel Climate Model
The Kiel Climate Model (KCM) GCM consists of the ECHAM5 atmospheric GCM coupled to the Nucleus for European Modeling of the Ocean (NEMO) ocean-sea-ice GCM (Park et al., 2009).The simulation runs from 126 to 115 ka BP and has been performed with 10 times accelerated astronomical forcing (see Sect. 2.3 for details).

LOVECLIM
The LOVECLIM EMIC includes a simplified atmospheric component and a low-resolution ocean GCM (Goosse et al., 2010).

MPI-UW
The MPI-UW (Max Planck Institute for Meteorology and University of Wisconsin-Madison) Earth system model (Mikolajewicz et al., 2007) is a GCM which consists of the ECHAM3 atmospheric GCM, the Large Scale Geostrophic Ocean (LSG2) GCM including a simple sea-ice model, the Hamburg Ocean Carbon Cycle model (HAMOCC) and the Lund Potsdam Jena dynamical terrestrial vegetation model (LPJ).The latter two components encompass the entire carbon cycle, which allows for the prognostic calculation of atmospheric pCO 2 from the fluxes of HAMOCC and LPJ.
In turn, the prognostic pCO 2 is then used in the radiative calculations resulting in the GHG forcing in this simulation to be different.The atmosphere in this MPI-UW simulation was integrated in periodically synchronous mode.In this method the slow components such as the ocean and vegetation are integrated for the entire time span, but the fast (and computationally expensive) atmospheric component is integrated only part of the time.In the meantime the ocean model is driven with fluxes from previous synchronous integration periods in combination with an EBM-type damping for small sea surface temperature anomalies.The main underlying assumption is that the atmosphere is in statistical equilibrium with the underlying sea-surface temperature and sea ice distribution.The LIG MPI-UW simulation runs from 128 to 115 ka BP (Schurgers et al., 2007;Gröger et al., 2007).Note that in this simulation the vegetation is not fixed at preindustrial values but actively simulated.

Evolution of the main climatic forcings of the LIG period
Having an overview of the changes in main climate forcings will allow us to identify if the simulations show a linear relation between changes in temperature and the climatic forcings or if feedback mechanisms within the climate system are important.The two climate forcings discussed here are the amount of insolation received by Earth and atmospheric GHG concentrations (Figs. 1 and 3; values according to the PMIP3 protocol, http://pmip3.lsce.ipsl.fr/).The changes in the amount of insolation received by Earth result from changes in the astronomical configuration.Globally averaged, the anomalies are close to zero for the LIG period.However, changes in the distribution over the different latitudes and seasons are not (Berger, 1978).During the first part of the LIG, the NH received more insolation in summer compared to the present-day period.Differences in insolation between both periods exceed 60 W m −2 for June at 65 • N (∼ 130-122 ka BP with a peak around ∼ 127 ka BP; Fig. 1).We do note that maximum insolation anomaly did not occur during the same period for all summer months (Fig. 3; Berger, 2001).For the Southern Hemisphere (SH) summer, the changes in insolation show a strong minimum between 126-123 ka BP while insolation values are just above the present-day ones during the late LIG (< 118 ka BP).The evolution of insolation during NH (SH) winter shows a maximum during the late (early) LIG.However, it is important to remember that, at high latitudes (>∼ 67 • ), absolute insolation and the insolation anomalies are close to zero in winter.
While all simulations in this model inter-comparison include changes in the astronomical configuration, only some include changes in GHG concentrations (Bern3D, CLIMBER-2, FAMOUS, LOVECLIM and MPI-UW).These are kept constant in the other simulations (CCSM3 and KCM).In accordance with the PMIP3 protocol, the FA-MOUS and LOVECLIM simulations include changes in CO 2 , CH 4 and N 2 O; the Bern3D simulation only includes changes in CO 2 and CH 4 ; and the CLIMBER-2 simulation only includes changes in CO 2 .The GHG concentration values for the LIG period in the PMIP3 protocol are based on ice-core data of Luthi et al. (2008), Loulergue et al. (2008) and Schilt et al. (2010) for respectively CO 2 , CH 4 and N 2 O.A linear interpolation was applied to these data in order to get a 1 yr time resolution.The GHG concentration changes result in a radiative forcing which increases from low values around 130 ka BP towards maximum values, just above pre-industrial, around 128.5 ka BP.This peak around 128.5 ka BP is sharp and short-lived (i.e.centennial time scale), after which the radiative forcing related to the GHG concentration changes remains approximately constant at values just below pre-industrial (Fig. 1; radiative forcing calculated after Houghton et al., 2001).Note that the radiative forcing provided by the changes in the three major GHGs is small, < 0.2 W m −2 , compared to the forcing provided by  2010) for CO 2 , CH 4 and N 2 O respectively.The corresponding pre-industrial values are given by the dotted lines.In constructing the forcing scenarios, a linear interpolation was applied to the data to get a 1 yr resolution.The bold black line presents the combined radiative forcing of the three main GHG concentration changes (W m −2 ; concentrations according to the PMIP3 protocol; formulation of radiative forcing after Houghton et al., 2001).The fixed average LIG GHG concentrations applied in the CCSM3 simulation are depicted by the arrows on the left-hand side of the upper panel.In the MPI-UW simulation, the radiative forcing of the GHGs is prognostically calculated from simulated changes in the carbon cycle.The resulting CO 2 changes are depicted in light green.In the lower part of the figure, the June insolation anomaly for 65 • N is given (W m −2 relative to pre-industrial values; Berger, 1978).
the insolation changes (Fig. 1).In the simulation performed with the MPI-UW model the pCO 2 concentration is a prognostic variable.Therefore the evolution of the GHG forcing is different from the other simulations (Fig. 1; note that CH 4 and N 2 O are neglected).The simulated GHG evolution in the MPI-UW simulation is characterized by a slow increase between 128-122 ka BP from ∼ 270 ppm towards more stable values of around 285 ppm between 122-115 ka BP.The KCM and CCSM3 simulations have GHG concentrations fixed at pre-industrial and average LIG values respectively (see Table 1

for details).
In all the simulations described in this study, a fixed-day calendar is used.Even though this is common practice in palaeoclimate simulations, it does have a non-negligible effect on the monthly temperature anomalies and has to be considered when comparing the results with proxy-based reconstructions (Joussaume and Braconnot, 1997;Chen et al., 2011).Another important uncertainty in this study is the age determination of the LIG (Kukla et al., 2002) and the related uncertainty of the applied forcings.However, since all simulations in this inter-comparison applied the same fixed-day calendar method and the same definition of the LIG period, these uncertainties have no impact on the results when we consider the robustness of the simulated temporal evolution of LIG temperatures and the differences between the various simulations (Joussaume and Braconnot, 1997).

Data processing
All the simulated temperature fields were averaged into 50 yr averages for every month.For the CCSM3 and KCM simulations, which are performed with a 10-fold acceleration of the changes in the insolation forcing, averaging over 50 astronomical years effectively means an average over only five model years.Therefore, in the CCSM3 and KCM simulations only sub-decadal climate variability is filtered out while in the other simulations the variability on multi-decadal time scales is also filtered out.Furthermore, the accelerated orbital forcings might distort the evolution of the deep ocean circulation due to its long response time (see also Sect. 3.5).In order to smooth out the artificial noise resulting from the fact that the MPI-UW model was integrated in periodically synchronous mode the temperature series have been averaged using a 200 yr window and afterwards linearly interpolated to obtain 50 yr averages.All temperatures of the different simulations were linearly re-gridded onto a common rectangular 1 • × 1 • grid.Throughout this manuscript, when dealing with "temperatures" we refer to differences between simulated LIG and pre-industrial near surface air temperature.The pre-industrial temperatures were obtained by averaging over the last 30-100 yr of long (> 500 yr) equilibrium simulations with pre-industrial values for the orbital parameters and GHG concentrations.

Results and discussion
In the first part of this section, we focus on the robustness among the models of the simulated LIG temperature evolution and the corresponding LIG temperature maximum.In the second part these results are then related to changes in insolation and GHG concentrations.In part three we discuss the importance of model complexity and resolution and in parts four to seven the impact of climate feedbacks on the simulated LIG temperature evolutions.An assessment of the inter-model differences and spatial patterns in the simulated temperature evolutions can provide valuable information about the climate system and the important internal feedback mechanisms operating during the LIG period.However, we acknowledge that temperature changes provide www.clim-past.net/9/605/2013/Clim.Past, 9, 605-619, 2013 only indirect indications, but a more in-depth investigation of the underlying causes of the inter-model differences is outside the scope of this manuscript.Nonetheless, we deem that highlighting potentially important climate feedbacks is an essential first step towards future sensitivity experiments.

Simulated robust LIG temperature evolution
Because of the seasonal and latitudinal differences in insolation, we investigate the robustness of the simulated LIG temperature evolution for both January and July in five different latitudinal bands: high and mid-northern latitudes, the tropical latitudes and the mid-and high southern latitudes (respectively 60 • N-90 • N; 30 • N-60 • N; 30 • S-30 • N; 60 • S-30 • S and 90 • S-60 • S).These specific latitudinal bands were chosen, because many important feedback processes in the climate system are roughly confined to these regions, e.g. the albedo and sea-ice feedback in the higher latitudes and the monsoon system in the tropical latitudes.With this approach we explicitly assume that longitudinal differences in temperature anomalies are small compared to latitudinal differences, an assumption of which we will discuss the validity of in the final part of this section.We focus on January and July temperatures to investigate changes in either the cold or warm season.However, we do note that these months do not always represent the warmest or coldest months for a given location.
For most latitudinal bands and for both January and July, the different models show LIG temperature evolutions which are comparable in their trend, magnitude of the temperature change and period of maximum warmth (Fig. 2).The robustness of the simulated temperatures is depicted by the multimodel mean (MMM) evolution, the corresponding standard deviation (STDEV) and the spread in the MMM period of maximum warmth (Fig. 3).In the high and mid-latitudes of the NH, we find robust peak July MMM temperature anomalies of respectively 0.3-3.7 K and 0.7-5.3K compared to preindustrial during the time intervals of respectively 128.4-125.1 ka BP and 129.4-126.3ka BP (temperature ranges and time intervals are given as the one standard deviation interval around the MMM; see caption of Fig. 3 for more details).Simulated January temperature anomalies for the NH midlatitudes also show a robust trend with a −0.8 to 2.1 K peak between 121.3-117.2ka BP.Finally, all simulations which include changes in GHG concentrations according to the PMIP3 protocol simulate a period of maximum July warmth for the SH high latitudes of −1.3 to 2.3 K between 129.2-126 ka BP.The period of maximum January and July LIG warmth and the corresponding temperature anomalies for all latitudinal bands are listed in Table 2.

Temperature evolution and forcings
The simulated evolution of MMM January and July temperature anomalies shows in most cases a clear correspondence with respectively December and June insolation (Fig. 3).Furthermore it is apparent that, when the rate of change of the radiative forcing is large, the model results are robust, but without such a strong trend the resulting temperature evolutions tend to differ largely.The latter is especially true for the simulated temperature evolutions in the winter months in high latitude regions which are characterized by a very small insolation forcing.Simulated July temperature anomalies for the SH exemplify another finding.We see that the models that include changes in GHG concentrations according to the PMIP3 protocol (Bern3D, CLIMBER-2, FAMOUS and LOVECLIM) tend to simulate a more distinct temperature evolution compared to the simulations that include fixed GHG concentrations (CCSM3 and KCM).The different GHG forcing evolution applied in the MPI-UW simulation mainly impacts the simulated temperatures of the SH high latitudes.These findings highlight the importance of the radiative forcing provided by the changes in GHG concentrations, though mainly when the magnitude and trend in the insolation forcing is small.

Model resolution and complexity
The presented model inter-comparison shows that, on long timescales (> 1 ka) and on large spatial scales (latitudinal bands for instance), the differences between EMICs and GCMs in model complexity and model resolution are of minor importance for the simulated temperature evolution.The results do show that the simulated high-frequency climate variability (50-100 yr timescales in this study because of the applied time-averaging) is much larger in the models of higher complexity and resolution.This is not only caused by the higher complexity of the equations of motion in the GCMs but also by the fact that in the simulations, which include accelerated orbital forcings, the 50 yr averages are actually resulting from only five model years instead of 50.Similarly, because in the MPI-UW simulation the atmosphere is integrated in a periodically synchronous way, the 50 yr averages include only a smaller number of years of atmospheric output, which results in more high-frequency variability.Finally, two of the models (CLIMBER-2 and MPI-UW) include a dynamical terrestrial biosphere component.However, its impact is not easily identified.Schurgers et al. (2007) describe that the vegetation-albedo feedback in MPI-UW causes higher temperatures in the mid-to high latitudes of the NH in the period 128-121 ka BP.Indeed, the MPI-UW temperatures in this latitude band are higher compared to most other models, but in the CLIMBER-2 temperature evolution this is not the case.

Sea ice and the LIG temperature evolution
While the astronomical and GHG forcings are longitudinally homogeneous, the net effect of changes in the radiative forcings and feedback mechanisms within the climate system can We define maximum warmth when the temperature is between the upper limit (which is the absolute temperature maximum) and the lower limit (which is the upper limit minus 10 % of the differences between the LIG maximum and minimum values).Note that these calculations are not performed on the 50 yr average-data but on the values from a polynomial which is fitted to the temperature-series in order to get a robust, multi-millennial signal.The individual period of maximum warmth is only shown if the maximum temperature along the polynomial is outside the 1 σ range of the LIG mean.Note that in all simulations shown in this figure a fixed-day calendar was used.
Table 2. List of the MMM timing of the period of maximum LIG warmth (ka BP) and the corresponding temperature anomalies (K) for every latitude band and for both January and July.The time interval of maximum warmth and the corresponding temperature anomalies are calculated as follows: the middle of the period is the time for which the maximum MMM temperature is found, which is in turn the middle temperature given in column seven; the maximum and minimum of the period of maximum warmth are the mid-point ±1 σ , which gives a measure of the spread in the periods of maximum warmth of the individual simulations (for details see captions of Figs. 2 and 3); the temperature range is found by taking the minimum and maximum of the MMM temperatures ±1σ occurring during the period of maximum warmth.In the last column it is given if the mid-point of the MMM period of maximum warmth is outside the 1 σ range of the LIG MMM average temperature.Note that, in all simulations underlying these results, a fixed-day calendar was used.result in longitudinal differences in the temperature evolution.To investigate the importance of feedback mechanisms, we focus not only on the temperature evolution of the individual models and the MMM but also on spatial patterns in the MWT for January and July (Figs. 4 and 5).
For the Arctic Ocean, we find an overall agreement on an early MWT in January (before approximately 123 ka BP) which we relate to a strong feedback involving sea-ice changes.It appears that, in line with findings for the Holocene (Renssen et al., 2005), the early LIG June insolation maximum results in a decline in the summer sea-ice cover and thickness and a subsequent decline in the winter sea-ice cover and thickness.This could in turn enhance the heat flux from the ocean to the atmosphere leading to higher atmospheric winter temperatures.Changes in winter insolation or in the GHG forcing are not likely the cause since the pattern is clearly confined to the Arctic.Moreover winter insolation changes in the region are close to zero, and the MPI-UW simulation does show an early January temperature optimum over the Arctic while it does not include an early LIG GHG forcing maximum.These findings provide strong indications that the sea-ice feedback plays an important role in determining the LIG winter temperature evolution in the Arctic region.We note, however, that the January MMM STDEV in parts of the Arctic is larger than it is over the surrounding regions, which indicates the model dependency of the seaice-temperature feedback (Fig. 5).
Interestingly, the situation in the sea-ice covered regions surrounding Antarctica is rather different.In contrast to the Arctic Ocean, in most of the simulations the MWTs in January and July for the sea-ice covered areas of the SH summer do not coincide.Rather, the July MWT in these areas is ∼ 127 ka BP, while the January MWT is after ∼ 118 ka BP.Also note that no robust January MWT is found over most of the Southern Ocean (Fig. 5) and that the early LIG July MWT over the high latitudes of the SH is only found in 4 of the models (see also Sect. 3.2).
These findings highlight the need for a more thorough investigation to establish the importance of sea-ice cover changes for the evolution of LIG temperatures and the possible differences between the Arctic and Antarctic.

The AMOC and the LIG temperature evolution
Several of the simulations in this inter-comparison study show large, abrupt changes in surface temperature anomalies that cannot easily be related to changes in the forcings.However, we are able to relate some of them to changes in the strength of Atlantic meridional overturning circulation (AMOC) by comparing the timing of these temperature changes with the simulated AMOC evolution (Fig. 6).The strength of the AMOC has a large impact on the exchange of heat between the SH and NH and on the exchange of heat between the atmosphere and the ocean.However, the strength of the AMOC can change rapidly because of its potentially bi-stable behaviour (Stommel, 1961).Moreover, its strength, stability and locations of the main convective regions are highly model-dependent.As a result, the simulated LIG evolution of characteristics of the AMOC differs largely between the models.In the following we will discuss the apparent relation between the simulated LIG temperatures and the strength of the AMOC.Note, however, that these are indirect indications only and additional sensitivity experiments would be necessary to confirm the connection and understand the underlying mechanisms.Note that the number of models used in calculating the MMM and STDEV differs per latitude band and per time interval, because Bern3D has latitudinal grid limits at 76 • N and 76 • S and the MPI-UW and KCM simulations run only from respectively 128-115 ka BP and 126-115 ka BP.Furthermore, the applied orbital acceleration in the CCSM3 and KCM and the applied periodically synchronous integration in MPI-UW impact the size of the STDEV (see also Sect.2.3).The vertical grey lines give the MMM mid-points of the period of maximum warmth.The width of the corresponding grey box gives the spread in the periods of maximum warmth of the individual simulations (1 σ ) while its height gives an indication of the number of simulations used in the calculation (see Fig. 2).This model spread is further illustrated by the depicted midpoints of the individual periods of maximum warmth (see Fig. 2 for the corresponding colour coding).The period of maximum warmth is only shown if for at least four individual simulations we find a maximum temperature along the polynomial which is outside the 1 σ range of the model individual LIG mean.Also shown in this figure are the corresponding 130-115 ka BP insolation anomalies for the mid-latitude of the chosen latitude bands (W m −2 ).The left-hand column gives the December (yellow), January (orange) and February (red) insolation anomalies, while in the right-hand column insolation anomalies for June (yellow), July (orange) and August (red) are given.All values are anomalies relative to pre-industrial values.The FAMOUS model simulates an abrupt decline in January NH temperatures around 121 ka BP (Fig. 2), a signal which originates largely from the northern Pacific (not shown here).This abrupt decline is followed within centuries by an abrupt increase of both January and July temperatures at mid-southern latitudes.Furthermore, the FAMOUS simulation shows an anomalous January MWT pattern over the Southern Ocean and over the northern and north-eastern Pacific (Fig. 4).The LIG simulation performed with the FA-MOUS model is characterized by two different modes of the AMOC.Between 130-121 ka BP a weak AMOC is simulated (∼ 70 % weaker compared to the strong mode of the AMOC; Fig. 6) and the main site of deep convection is located in the North Pacific (not shown).Around 121 ka BP the AMOC switches to a strong mode which is characterized by deep convection mainly taking place in the North Atlantic.The changes in the global overturning circulation around 121 ka BP are likely related to the concurrent temperature changes and the anomalous MWT in parts of the North Pacific and Southern Ocean, and the mechanisms behind it should be investigated more thoroughly in a future study.
The second simulation showing abrupt temperature changes which can be related to changes in the AMOC strength is LOVECLIM.The simulation performed with LOVECLIM shows an abrupt January cooling at around 120 ka BP at high northern latitudes (Fig. 2).Accompanying this temperature shift is a region of anomalous timing of winter warmth in the Labrador Sea (Fig. 4).The evolution of the AMOC strength in the LOVECLIM simulation is characterized by an abrupt weakening by about 15 % around 120 ka BP (compared to the period before), after which the AMOC remains unstable (Fig. 6).The simulated changes in AMOC strength in the LOVECLIM model are related to periods of weakened convection in the Irminger Sea (not shown here).Two quasi-stable AMOC modes have already been described for the LOVECLIM model by Schulz et al. (2007).An interesting question is whether the mode transitions in the FAMOUS and LOVECLIM simulations are determined by an external forcing or do they occur by chance (stochastically) both at around 120 ka BP.Running the ensemble simulations necessary to answer this question is however outside the scope of the manuscript.
The Bern3D model shows a anomalous spatial pattern over the Labrador Sea with a timing of maximum winter warmth clearly offset from the surrounding areas.As mentioned before, the Bern3D model incorporates both prescribed remnant ice sheets and a related meltwater flux entering into the ocean between roughly 130-125 ka BP (see Sect. 2.1.1).Concurrent with this melt flux, a somewhat weakened AMOC is simulated for the period 129-125 ka BP (20-30 % compared to the period 125-121 ka BP during which no freshwater flux is applied; Fig. 6).After 121 ka BP, growth of the NH continental ice sheets is prescribed in the Bern3D simulation.To compensate for this, a volume of freshwater is removed globally from the ocean surface.Therefore the surface ocean becomes more dense (not shown), and as a result an increase of the AMOC strength is simulated.The timing of maximum winter warmth over the Labrador region, after 117 ka BP, corresponds to the period of maximum AMOC strength.Note however that, compared to the simulations performed with FAMOUS and LOVECLIM, the simulated changes in the AMOC strength in Bern3D do not seem to have such a clear impact on the simulated LIG temperature evolution (Fig. 2).
The last simulation showing abrupt changes in the LIG temperature evolution which can be related to changes in the AMOC strength is CLIMBER-2.In the period between 123-120 ka BP, we see large temperature fluctuations with a duration of about 1 ka (Fig. 2).Similar shifts are apparent in the simulated strength of the AMOC in the CLIMBER-2 simulation (Fig. 6).Interestingly and in contrast to the three simulations described above, while the changes in AMOC strength are relatively small, they cause major, ∼ 4 K shifts in the NH mid-latitude temperatures.
Of the remaining three models, the KCM and CCSM3 simulations show a more stable AMOC throughout the LIG.But note that the KCM and CCSM3 simulations have been performed with an acceleration technique, and therefore this analysis does not allow us to determine if the apparent stability of the AMOC strength is an actual model characteristic or an artefact of the applied acceleration technique.Finally, the MPI-UW simulation shows fluctuations in the AMOC strength (Fig. 6) and corresponding high-frequency temperature variability in the SH (Fig. 2) which are related to multi-centennial variability in Southern Ocean deep convection (not shown).
Regardless of the exact mechanisms causing the changes in the AMOC in the different simulations, these results show the importance of the evolution of the configuration and strength of the AMOC for the simulated evolution of LIG temperature anomalies.This appears to be especially true for regions like the North Atlantic, the Labrador Sea, the Norwegian and the Barents seas, the north-eastern Pacific and possibly the Southern Ocean.In order to simulate a more robust LIG temperature evolution for these regions, stronger constraints are needed on how the configuration of the AMOC evolved and if mode-switches occurred during the LIG, which is, to date, still inconclusive (Nieuwenhove et al., 2011).In the transient LIG simulation performed with the Bern3D model, changes in the size of the continental NH ice sheets were prescribed.According to the method applied by Ritz et al. (2011a) to reconstruct remnant ice sheets, parts of the NH continents remained glaciated until ∼ 125 ka BP (see Sect. 2.1.1).Therefore, this simulation provides the possibility to investigate the impact of remnant ice sheets on the LIG temperature evolution.The impact of remnant ice is most clearly visible in the simulated July temperature evolution in high northern latitudes.The peak warmth simulated by the Bern3D model occurs several thousands of years later than in most of the other simulations (Fig. 2).Note, however, that it is not easy to distinguish between the impact of remnant ice and the weakening of the AMOC.It is to be expected that remnant ice has a larger impact on July temperatures, while the changes in the strength of the AMOC more likely result in changes in January temperatures.However, this is not clear from our results.Renssen et al. (2009) have shown for the present interglacial period that remnant ice sheets had a profound influence on the surrounding regions, delaying the thermal maximum to several thousands of years after the insolation optimum.The simulation performed with the Bern3D model, even though of lower resolution than the model used by Renssen et al. (2009), indicates a similar delay of maximum warmth.However, a thorough comparison of several simulations including remnant ice sheets is needed to retrieve a more robust signal of the impact of remnant ice on the LIG temperature evolution.Such a model inter-comparison is however further complicated by the fact that different reconstructions of the changes in altitude and extent of the ice sheets during the LIG are still inconclusive.According to Kopp et al. (2009), maximum LIG global sea level was at least 6 m above present day.A small part of this can be accounted for through thermal expansion of the ocean waters (McKay et al., 2011) and a loss of mountain glaciers.But the relative contributions of the Greenland and Antarctic ice sheets are still heavily debated (Dutton and Lambeck, 2012, and references therein).

The monsoon and the simulated LIG temperature evolution
Finally we examine the anomalous pattern in simulated LIG temperatures in the regions centered on the Sahel and on India.The simulated July MMM MWT over these regions is clearly later than the surrounding regions (Fig. 5).However, the large STDEV value indicates that this anomalous pattern is only simulated by some of the models.Maximum July temperatures over the regions centered on the Sahel and India are reached before 126 ka BP in Bern3D and LOVECLIM, around 124 ka BP in CLIMBER-2 and after 120 ka BP in CCSM3, KCM and MPI-UW.Again another pattern is simulated by the FAMOUS model, showing an early (> 127 ka BP) MWT over the Sahel but a late (< 118 ka BP) MWT over parts of India (Fig. 4).Solely based on the geographical distribution of these anomalous MWT patterns, we relate them to changes in the monsoon system.There are several important differences between the different models which can at least partly explain the large STDEV in these areas.Most of the GCMs in this model inter-comparison (CCSM3, KCM and MPI-UW) simulate a MWT in the Indian and African monsoon regions which is delayed with respect to the insolation forcing.This can be explained by strong feedbacks between the insolation, land evapotranspiration and cloudiness as described by Mulitza et al. (2008).In the MPI-UW simulation this negative feedback related to clouds and precipitation is partly compensated by the positive vegetation-albedo feedback in these regions as savanna is partly replaced by tropical and temperate forests (not shown).The EMICs in this model intercomparison have difficulty to realistically simulate changes in the monsoon system because of their low resolution and simplified atmospheric physics and dynamics.For instance in the Bern3D model, the moisture transport is driven by fixed, zonally averaged winds which inhibit any changes in the monsoon system to be simulated.In LOVECLIM, the simulated changes in the tropical regions should also be treated with care since the results are strongly affected by the quasi-geostrophic nature of its atmosphere and the fixed cloud cover.This model inter-comparison shows the importance of changes in monsoon systems for the evolution of LIG temperatures in regions centered on the Sahel and India even though the exact climate change mechanisms at work for these specific regions are likely different in the different simulations.Furthermore, the model inter-comparison has shown that the LIG temperature evolution in monsoon regions tends to be very different between EMICs and GCMs.

Summary
In this manuscript we have presented the first model intercomparison study of long, > 10 ka transient simulations covering the LIG period with a total of seven different climate models.The aim has been to determine the common transient temperature response to the LIG forcings and to indicate which feedbacks appear to play a crucial role.Despite large differences between the incorporated climate models and the forcings of the different simulations, the results show for large parts of the globe a robust evolution of LIG January and July surface air temperature anomalies compared to pre-industrial values.The main findings are outlined below.
-Simulated July temperature anomalies for the NH show a robust temperature maximum between 130-125 ka BP with a magnitude ranging from 0.3-5.3K.
-Peak January warmth is found for the NH mid-latitudes (30 • N-60 • N) between 121-117 ka BP with anomalies of −0.8 to 2.1 K, but for the NH high latitudes (60 • N-90 • N) the results of the simulations are inconclusive.
-Over the Arctic Ocean a robust 128-126 ka BP timing of peak January warmth is found.
-The simulated January temperature evolution for the SH shows a robust temperature maximum after 121 ka BP with corresponding temperature anomalies of −1 to 1.2 K.
-The July temperature maximum for the SH is found between 130-124 ka BP but only in the four simulations which include prescribed changes in GHG concentrations.
-The simulations show a close relation between changes in insolation and temperature in both summer and winter for almost all latitudes.There are two main exceptions: NH high-latitude winter temperature changes result largely from sea-ice related feedbacks and are therefore highly model-dependent; and SH mid-to highlatitude winter temperature changes appear strongly affected by changes in GHG concentrations.
-The impact of model complexity and or resolution is of minor importance for the simulated temperature evolution over longer (> 1 ka) timescales and large spatial scales.
Based on the differences between the simulations and the investigation of regional patterns in the timing of maximum warmth, we found that several climate feedbacks are likely to be of major importance when simulating the evolution of LIG temperatures: the sea-ice feedback, changes in the strength and configuration of the AMOC, remnants of NH continental ice sheets and changes in the monsoon system.
Our results yield important information for future modeldata comparison studies by detailing the following: (1) in which regions the model results are robust and can thus potentially help interpret temperature reconstructions; (2) for which regions more constraints on the LIG climate are needed before a successful model-data comparison can be achieved, because the simulated temperature changes are closely linked to model-dependent climate feedbacks.Moreover, our findings will serve as the starting point for a number of sensitivity studies that will be performed to determine exactly how important these feedbacks are and to investigate the mechanisms behind them.

Fig. 1 .
Fig. 1.In the upper part the changes in the main GHG concentrations over the period 130-115 ka BP are depicted.Concentrations are according to the PMIP3 protocol (http://pmip3.lsce.ipsl.fr/),which is based on the ice-core data of Luthi et al. (2008), Loulergue et al. (2008) and Schilt et al. (2010) for CO 2 , CH 4 and N 2 Orespectively.The corresponding pre-industrial values are given by the dotted lines.In constructing the forcing scenarios, a linear interpolation was applied to the data to get a 1 yr resolution.The bold black line presents the combined radiative forcing of the three main GHG concentration changes (W m −2 ; concentrations according to the PMIP3 protocol; formulation of radiative forcing afterHoughton et al., 2001).The fixed average LIG GHG concentrations applied in the CCSM3 simulation are depicted by the arrows on the left-hand side of the upper panel.In the MPI-UW simulation, the radiative forcing of the GHGs is prognostically calculated from simulated changes in the carbon cycle.The resulting CO 2 changes are depicted in light green.In the lower part of the figure, the June insolation anomaly for 65 • N is given (W m −2 relative to pre-industrial values;Berger, 1978).

Fig. 2 .
Fig. 2. Simulated 130-115 ka BP January and July surface air temperature anomalies (K) by the seven climate models for five different latitude bands.All temperatures are anomalies relative to pre-industrial values.The temperature series are 50 yr averages.The horizontal bars accompanying the surface temperature anomalies are the periods of maximum warmth for each individual simulation.We define maximum warmth when the temperature is between the upper limit (which is the absolute temperature maximum) and the lower limit (which is the upper limit minus 10 % of the differences between the LIG maximum and minimum values).Note that these calculations are not performed on the 50 yr average-data but on the values from a polynomial which is fitted to the temperature-series in order to get a robust, multi-millennial signal.The individual period of maximum warmth is only shown if the maximum temperature along the polynomial is outside the 1 σ range of the LIG mean.Note that in all simulations shown in this figure a fixed-day calendar was used.

Fig. 3 .
Fig. 3. Simulated 130-115 ka BP January and July MMM surface air temperature anomalies (K; black) and STDEV (1 σ ; grey) for five different latitude bands.The MMMs and STDEVs are calculated using the temperature time series shown in Fig. 2.Note that the number of models used in calculating the MMM and STDEV differs per latitude band and per time interval, because Bern3D has latitudinal grid limits at 76 • N and 76 • S and the MPI-UW and KCM simulations run only from respectively 128-115 ka BP and 126-115 ka BP.Furthermore, the applied orbital acceleration in the CCSM3 and KCM and the applied periodically synchronous integration in MPI-UW impact the size of the STDEV (see also Sect.2.3).The vertical grey lines give the MMM mid-points of the period of maximum warmth.The width of the corresponding grey box gives the spread in the periods of maximum warmth of the individual simulations (1 σ ) while its height gives an indication of the number of simulations used in the calculation (see Fig.2).This model spread is further illustrated by the depicted midpoints of the individual periods of maximum warmth (see Fig.2for the corresponding colour coding).The period of maximum warmth is only shown if for at least four individual simulations we find a maximum temperature along the polynomial which is outside the 1 σ range of the model individual LIG mean.Also shown in this figure are the corresponding 130-115 ka BP insolation anomalies for the mid-latitude of the chosen latitude bands (W m −2 ).The left-hand column gives the December (yellow), January (orange) and February (red) insolation anomalies, while in the right-hand column insolation anomalies for June (yellow), July (orange) and August (red) are given.All values are anomalies relative to pre-industrial values.

Fig. 4 .
Fig. 4. Timing of the simulated January and July LIG temperature maximum (MWT) in the different simulations.The timing (ka BP) is the period for which the highest 50 yr average temperature anomalies are found along a polynomial fitted to the individual temperature time series for every grid cell.Values are only shown if the highest temperatures are outside the 2 σ range of the mean LIG temperature anomaly for the specific grid cell.Note that Bern3D has latitudinal grid limits at 76 • N and 76 • S. Furthermore, the MPI-UW and KCM simulations run from 128-115 ka BP and 126-115 ka BP respectively.Therefore, the colour corresponding to a timing of respectively 128-127 ka BP and 126-125 ka BP should be interpreted as being > 127 ka BP or > 125 ka BP in the case of the MPI-UW and KCM simulations.In all simulations shown in this figure, a fixed-day calendar was used.

Fig. 5 .
Fig. 5. MMM (left) and STDEV (right) of the timing of maximum warmth (MWT) for January and July based on the individual simulations as shown in Fig. 4. The standard deviation gives a measure of the spread in the simulated timing of maximum warmth between the different models (ka; 2 σ ).The MMM and STDEV are only shown if in at least four individual simulations we find a maximum temperature along a fitted polynomial which is outside the 1 σ range of the model individual LIG mean (see Fig. 4).Note that the calculations do not include the Bern3D model north of 76 • N and south of 76 • S in accordance with its latitudinal grid limits.The calculations do include the KCM and MPI-UW simulations.However, since they only run from 126-115 ka BP and 128-115 ka BP respectively, this might slightly offset the results shown in this figure.

Fig. 6 .
Fig. 6.Simulated LIG evolution of the strength of the AMOC in the seven different simulations.The maximum overturning stream function in the North Atlantic (Sv) is used to indicate the strength of the AMOC.All values are averages over 50 astronomical years.The arrows in the top panel indicate the moment when the meltwater flux from the remnant ice sheets into the oceans in the Bern3D simulation ceases (∼ 125 ka BP) and when continental ice sheets on the NH start to expand again (∼ 121 ka BP).