Effects of melting ice sheets and orbital forcing on the early Holocene warming in the extratropical Northern Hemisphere

The early Holocene is marked by the final transition from the last deglaciation to the relatively warm Holocene. Proxy-based temperature reconstructions suggest a Northern Hemisphere warming, but also indicate important regional differences. Model studies have analyzed the influence of diminishing ice sheets and other forcings on the climate system during the Holocene. 15 The climate response to forcings before 9 ky (refers to ky B.P), however, remains not fully comprehended. We therefore, by employing the LOVECLIM climate model, studied how orbital and ice-sheet forcings contributed to climate change and to these regional differences during the earliest part of the Holocene (11.5–7 ky). Our equilibrium experiment for 11.5 ky suggests lower annual mean temperatures at the onset of the 20 Holocene than in the preindustrial era with the exception of Alaska. The magnitude of this cool anomaly varied regionally and these spatial patterns are broadly consistent with proxy-based reconstructions. Temperatures throughout the whole year in northern Canada and northwestern Europe for 11.5 ky were 2–5°C lower than those of the preindustrial era as the climate was strongly influenced by the cooling effect of the ice sheets, which was caused by enhanced surface albedo and 25 ice-sheet orography. In contrast, temperatures in Alaska for all seasons for the same period were 0.5–3°C higher than the control run, which were caused by a combination of orbital forcing and stronger southerly winds that advected warm air from the South in response to prevailing high air pressure over the Laurentide Ice Sheet (LIS). The transient experiments indicate a highly inhomogeneous early Holocene temperature warming 30


Introduction
The early Holocene from 11.5 to 7 kyr BP (hereafter noted as kyr) is paleoclimatologically interesting as it represents the last transition phase from full glacial to interglacial conditions.This period is characterized by a warming trend in the Northern Hemisphere (NH) that has been registered in numerous proxy records and indicated by stacked temperature reconstructions (Shakun et al., 2012).Oxygen isotope measurements from ice cores in Greenland (Dansgaard et al., 1993;Grootes et al., 1993;Rasmussen et al., 2006;Vinther et al., 2006Vinther et al., , 2008) ) and the Canadian High Arctic (Koerner and Fisher, 1990) consistently show an increase in δ 18 O by up to 3-5 ‰, which indicates an approximate warming in the climate system (Vinther et al., 2009).Moreover, this early Holocene warming is also registered in biological proxies.For example, a 4-5 • C warming in western and northern Europe is indicated by chironomid and macrofossil data obtained from lake sediments (Brooks and Birks, 2000;Brooks et al., 2012;Birks, 2015).In addition, this transition is recorded in other high-resolution records from further east in Eurasia, such as in the speleothems from China (Yuan et al., 2004;Wang et al., 2005).Comparable trends have been identified in marine sediment core data, such as sea surface temperature (SST) rise in the North Atlantic reflected by the variation in δ 18 O and in planktonic foraminifera (Bond et al., 1993;Kandiano et al., 2004;Hald et al., 2007).Although these proxy records provide a general view of early Holocene warming, their detailed expression in different regions and the reasons for this spatial variation are poorly known.
The orbitally induced increase in NH June insolation was one of the main external drivers of climate change during the last deglaciation (Berger, 1988;Denton et al., 2010;Abe-Ouchi et al., 2013;Buizert et al., 2014).This increase peaked in the earliest Holocene (Berger, 1978) and resulted in warming over large areas.However, the early Holocene was also characterized by adjustments in components of the climate system that further affected the temperature through various feedback mechanisms.In the cryosphere, the Laurentide Ice Sheet (LIS) and Fennoscandian Ice Sheet (FIS) were melting at a fast rate and eventually disappeared at around 6.8 and 10 kyr, respectively (Dyke et al., 2003;Occhietti et al., 2011), which exerted multiple influences on the climate system (Renssen et al., 2009).First, the surface albedo was much higher over the ice sheets compared to ice-free surfaces, which resulted in relatively low temperatures.Second, the ice-sheet topography could have also influenced the climate through the mechanism of adjustment to the atmospheric circulation (Felzer et al., 1996;Justino and Peltier, 2005;Langen and Vinther, 2009).For instance, a large-scale ice sheet could have generated a glacial anticyclone that locally could have further reduced the temperature (Felzer et al., 1996), but it may also have caused a 2-3 • C warming over the North Atlantic in the Last Glacial Maximum (LGM) (Pausata et al., 2011;Hofer et al., 2012).Third, both modeling and proxy studies have found that the Atlantic Meridional Overturning Circulation (AMOC) was relatively weak during the early Holocene due to the ice-sheet melting, which led to reduced northward heat transport and extended sea-ice cover (Renssen et al., 2010;Roche et al., 2010;Thornalley et al., 2011Thornalley et al., , 2013)).Overall, the net effect of ice sheets on the early Holocene climate can be expected to have tempered the orbitally induced warming at the mid-and high latitudes.Important adjustments in the carbon cycle occurred in the early Holocene, as evidenced by the rise in atmospheric CO 2 levels by 20-30 ppm that contributed to the warming (Schilt et al., 2010).Changes also happened in the biosphere during the early Holocene.Vegetation reconstructions revealed a northward expansion of boreal forest in the circum-Arctic region after the retreat of the ice sheets (MacDonald et al., 2000;Bigelow et al., 2003;CAPE project, 2001;Fang et al., 2013).This expansion of boreal forest into regions that were not previously vegetated or were covered by tundra caused a reduction of the surface albedo and induced a positive feedback to the warming trend (Claussen et al., 2001).
The impact of these forcings on the Holocene climate has been examined in modeling studies.The focus in these studies has been on the influence of the decay of the LIS and Greenland Ice Sheet (GIS) on the climate after 9 kyr relative to other climate forcings (Renssen et al., 2009;Blaschek and Renssen, 2013).Renssen et al. (2009) used transient simulations performed with the ECBilt-CLIO-VECODE model and found that the Holocene climate was sensitive to the ice sheets and that the LIS cooling effects delayed the Holocene Thermal Maximum (HTM) by up to thousands of years.Blaschek and Renssen (2013) applied a more recent version of the same model (renamed to LOVECLIM) and revealed that the GIS melting had an identifiable impact on the climate over the Nordic Seas.However, these Holocene modeling studies only started at 9 kyr.The most important challenges in simulating climate during the initial phase of the early Holocene are the inherent uncertainties in the icesheet forcings in terms of the ice-sheet dynamics and the related meltwater release.Recent deglaciation studies based on cosmogenic exposure dating indicate slightly older ages of deglaciation in some regions than suggested by radiocarbon dating data (Carlson et al., 2014;Clark, unpublished data), primarily because of a large uncertainty in bulk organic sample ages and the possibility of old carbon contamination (Carlson et al., 2014;Stokes et al., 2015).Furthermore, the Younger Dryas stadial ended at 11.7 kyr and may still have influenced the early Holocene climate due to the long response time of the deep ocean (Renssen et al., 2012).Therefore, the climate system's response to forcings before 9 kyr, especially those of the ice sheets, is poorly understood.
We have extended the study of Blaschek and Renssen (2013) back to 11.5 kyr to explore the early Holocene climate response to these key forcings.By employing the same climate model of intermediate complexity LOVECLIM, we first analyzed the impact of forcings on the climate at 11.5 kyr and subsequently investigated the influence of two ice-sheet deglaciation scenarios in transient simulations.The comparison of these different simulations enables us to disentangle how the ice sheets influenced the early Holocene climate.More specifically, we have addressed the following research questions.(1) What were the spatial patterns of simulated temperature at the onset of the Holocene (11.5 kyr)?(2) What were the roles of the forcings, especially ice-sheet decay, in shaping these features?(3) What was the spatiotemporal variability in the simulated early Holocene temperature evolution?

The LOVECLIM model
We conducted our simulations with version 1.2 of the threedimensional Earth system model of intermediate complexity LOVECLIM (Goosse et al., 2010), in which the components of the atmosphere, ocean including sea ice, vegetation, ice sheets and carbon cycle are dynamically included.However, in our version, the components for the ice sheets and the carbon cycle were not activated.Therefore, the ice-sheet evolution and greenhouse gases were prescribed in our present study.The atmospheric component is the quasi-geostrophic model ECBilt, which consists of three vertical layers and has T 21 horizontal resolution (Opsteegh et al., 1998).CLIO is the ocean component, which consists of a free-surface, primitive-equation oceanic general circulation model (GCM) coupled to a three-layer dynamic-thermodynamic sea-ice model (Fichefet and Maqueda, 1997).The ocean model includes 20 vertical levels and a 3 • × 3 • latitude-longitude horizontal resolution (Goosse and Fichefet, 1999).These two core components were further coupled to the biosphere model VECODE, which simulates the dynamics of two main terrestrial plant functional types, trees and grasses, in addition to desert (Brovkin et al., 1997).More details on LOVE-CLIM can be found in Goosse et al. (2010).
The LOVECLIM model is a useful tool to explore the mechanisms behind climate change, and it has made critical contributions to our understanding of the climate history and variabilities observed in proxy records (Renssen et al., 2005(Renssen et al., , 2006(Renssen et al., , 2010)).For example, it has helped with investigations of the potential forcings behind abrupt climate events (Renssen et al., 2002;Wiersma and Renssen, 2006;Renssen et al., 2015) as well as understanding the role of the decaying LIS and GIS in temperature evolution over the last 9 kyr (Renssen et al., 2009;Blaschek and Renssen, 2013).Moreover, the LOVECLIM model simulates a reasonable modern climate (Goosse et al., 2010).It also simulates the meridional overturning stream function reasonably well and reproduces a large-scale structure of atmosphere circulation that agrees with observations and with other models (Goosse et al., 2010).In addition, the model's sensitivity to freshwater perturbation is reasonable compared to that of other models (Roche et al., 2007), and its sensitivity to a doubling of atmospheric CO 2 concentration is 2 K, which is at the lower end of coupled general circulation model (GCM) estimates (Flato et al., 2013).

Prescribed forcings
We included the major climatic forcings in terms of greenhouse gases (GHGs) in the atmosphere, astronomical parameters (orbital forcing, or ORB) and decaying ice sheets.In all simulations, the solar constant, aerosol levels, the continental configuration and bathymetry were kept fixed at preindustrial values.We based the concentrations of CO 2 , CH 4 and N 2 O on ice core measurements for GHG forcing (Loulergue et al., 2008;Schilt et al., 2010).The radiative GHG forcing anomaly (relative to 0 kyr) in W m −2 (Ramaswamy et al., 2001), representing the overall GHG contribution, at first showed a rapid rise with a peak of −0.3 W m −2 at 10 kyr, which was followed by a slight decrease towards a minimum at 7 kyr, and gradually increased towards 0 kyr (Fig. 1).The astronomical parameters (eccentricity, obliquity and longitude of perihelion) determine the incoming solar radiation at the top of atmosphere and were derived from Berger (1978).An example of the resulting change in insolation is shown as the anomaly for June at 65 • N in Fig. 1, which shows the gradual decrease over the course of the Holocene.While the global annual mean insolation stayed at almost the same level (not shown), both changes in obliquity and precession are resulting in insolation variations on the multi-millennial timescale of the Holocene.At the beginning of the Holocene (11.5 kyr), the orbitally induced insolation anomaly in the NH was positive in summer and negative in winter (Fig. S1 in the Supplement).Overall, this setup of GHG and ORB forcing is in line with the PMIP3 protocol ((http://pmip3.lsce.ipsl.fr),except that our simulation excluded the increase in GHG levels during the industrial era (Ruddiman, 2007).Accordingly, the terms preindustrial (era) and 0 kyr are considered equivalent in the present text and indicate modern conditions without anthropogenic impacts.
We took three aspects into account concerning the icesheet forcing, namely their spatial extent, their thickness and their meltwater discharge.The reconstructions of ice-sheet spatial extent were based on the dating of geological features and on the correlation of these geological data sets between different regions (Dyke et al., 2003;Svendsen et al., 2004;Putkinen and Lunkka, 2008).According to these reconstructions, the FIS at 11.5 kyr covered most of Fennoscandia except for southern Scandinavia and eastern Finland (Svendsen et al., 2004;Putkinen and Lunkka, 2008;Clark, unpublished data).The LIS occupied most of the lowland area north of the Great Lakes region and filled the whole Hudson Strait (Licciardi et al., 1999;Dyke et al., 2003;Occhietti et al., 2011).The thickness of LIS was up to 2000 m, and for the FIS this thickness was only about 100 m (Ganopolski et  The early Holocene warming in the extratropical Northern Hemisphere q q q q q q q q q q q q q q q q q q q q 11 10 9 8 7 0 10 20 30 40 50 60 Ice sheet area and max thickness Time (kyr BP) Area of ice sheets (10E+5 km 2 ) q q q q q q q q q q q q q q q q q q q q q q LIS_size FIS_size LIS_tck FIS_tck 2010), which is comparable with the ICE-5G reconstruction (Peltier, 2004).Both the spatial extent of the ice sheets and their thickness were updated every 250 years in our transient experiments, and they decreased rapidly during the earliest Holocene, followed by a more gradual deglaciation rate from 8 kyr onward (Fig. 2a).
We applied the meltwater release for 1200 years in our equilibrium experiments for 11.5 kyr by adding 0.11 Sv (1 sverdrup is 10 6 m 3 s −1 ) of freshwater at the St. Lawrence River, 0.05 Sv at Hudson Strait and Hudson River, 0.055 Sv coming from the FIS, and 0.002 Sv from the GIS (Liccia-rdi et al., 1999;Jennings et al., 2015).The total freshwater volume added to the oceans in our transient experiments was about 1.46 × 10 16 m −3 in the first 4700 years (Fig. 2b), which roughly matches the estimated ice-sheet melting volume during the early Holocene (Dyke et al., 2003;Ganopolski et al., 2010;Clark, unpublished data).The volume of meltwater was slightly lower than the volume of the estimated 60 m sea level rise that took place during the early Holocene (Fig. 2c) (Lambeck et al., 2014), which suggests a coeval Antarctic melting contribution that is not considered here.Given the lack of a direct imprint left by meltwater on terrestrial records and hence the relatively large uncertainty, we used two versions of the freshwater flux (thick dashed lines and solid lines with symbols in Fig. 2b) that represent two possible deglaciation scenarios of the GIS and FIS, named FWF-v1 and FWF-v2.The GIS FWF_v1 scenario is derived from the ICE_5G reconstruction, and FWF_v2 is based on the reconstruction of Vinther et al. (2009), which suggests a faster GIS thinning.The two FIS freshwater flux (FWF) scenarios are based on two estimations of the FIS melting, since the recent cosmogenic dating (FWF_v2) supports a faster melting (Clark, unpublished data) than previously thought (FWF_v1).However, we kept the freshwater discharge from LIS the same as in version 1, since the LIS deglaciation has been relatively well studied and we are more certain about its contribution.

Setup of experiments
We performed two types of experiments: equilibrium and transient simulations.First, the equilibrium experiments of OG11.5 and OGIS11.5 with boundary conditions for 11.5 kyr (Table 1) were designed.The OGIS11.5 experiment included ice-sheet forcing, whereas no ice sheets were included in OG11.5 (Table 2).Each of these experiments was initiated from the model's default modern condition and was run for 1200 years, of which the last 200 years of data was used for the analysis.Renssen et al. (2006) demonstrated that a 1200year spin-up is sufficient to reach a quasi-equilibrium in all components of the model.
The data from the end of the 1200-year equilibrium run were then taken to initialize the transient experiments that covered the last 11.5 thousand years.In the first transient simulation named ORBGHG, both GHG and ORB varied on an annual basis.In the second simulation OGIS_FWF-v1, the ice-sheet topography (Fig. 2a) and FWF_v1 (thick dashed lines in Fig. 2b) were additionally included.A third experiment (named OGIS_FWF-v2) was performed with the freshwater version 2 (solid lines with symbols in Fig. 2b) and with the same ice-sheet topography as in OGIS_FWF-v1 to further investigate the climate response to the relatively uncertain freshwater forcing.Both OGIS_FWF-v1 and OGIS_FWF-v2 were initialized from the OGIS11.5 experiment.A preindustrial simulation (PI) was run for 1200 years from the model's default (representing modern conditions)  with the boundary conditions that are shown in Table 1 and, similarly to the other equilibrium experiments, the results of the last 200 years were used as a reference.These simulations and their forcings are summarized in Table 2.All temperature values in this study are shown as deviations from the PI simulation (indicates the climate at 0 kyr).Temperatures presented here are simulated near-surface temperature values without the environmental lapse rate corrections to the sea level temperature, which imply approximately a 0.5 • C cold bias over ice-sheet covered regions when compared with site-specific proxy records.

Equilibrium experiments at the onset of the Holocene
3.1.1Simulation with only ORB and GHG forcings at 11.5 kyr (OG11.5) In the experiment OG11.5, summer temperatures were 2-4 • C higher over most of the extratropical continents than in the PI simulation, with a maximum deviation of 5 • C in the central parts of the Northern Hemisphere continents (Fig. 3a).The warming over the oceans was about 1.5 • C, and less conspicuous than that over the continents.These warmer conditions were caused by the orbitally induced positive summer insolation anomaly, as all atmospheric greenhouse gas levels were lower at 11.5 kyr than in the preindustrial era (Fig. 1).The most obvious feature of simulated winter temperatures was the marked contrast between high latitudes and areas more to the south (Fig. 3b).For instance, midlatitudes were 1.5-3 • C cooler with the strongest cooling in the central continents, whereas the high-latitude Arctic was clearly warmer with a maximum up to +3 • C than the PI.This latitudinal gradient can be seen in annual mean temperatures as well (Fig. 3c).Annual mean temperatures over the Arctic were about 1-4 • C higher than the PI.The warming was slightly larger in winter than in summer, and this seasonal difference mirrors the Arctic Ocean damping effect on a seasonal signal due to a large heat capacity.Temperatures at lower latitudes were annually unchanged (mostly within ±0.5 • C) with a stronger seasonality (with warmer summers and cooler winters), which is consistent with the insolation change at 11.5 kyr (Fig. S1).
3.1.2Climate response to melting ice sheets at 11.5 kyr (OGIS11.5) Our simulation OGIS11.5 (including the impact of ice sheets) suggests a much cooler climate than that of OG11.5.Most notably, ice sheets induced a strong summer cooling over ice-covered areas, and reduced temperatures up to 5 • C compared to the PI simulation, with the strongest cooling at the center of the LIS (Fig. 3d).Additionally, SSTs were also more than 1.5 • C lower over the North Atlantic Ocean.In contrast, temperatures over the ice-free continents were mostly above the preindustrial level, but still lower than those found in OG11.5, except for the Alaska region.Although the area with colder conditions clearly expanded more in the OGIS11.5 simulation than in OG11.5, the central Arctic was still warmer in OGIS11.5 relative to PI (Fig. 3e).Alaska was the only continental region where winter temperatures exceeded the preindustrial values by up to +3 • C. The strongest cooling effect was present in the regions covered by ice sheets, for instance more than 3 • C cooler over the LIS.Simulated annual mean temperatures in OGIS11.5 clearly showed overall lower values than the PI due to the icesheet impacts (Fig. 3f).The Eurasian continent was mostly 1.5-3 • C cooler, and a maximum temperature reduction of more than 5 • C was found over the LIS.Only two areas were still warmer: Alaska, including the adjacent sector of the Arctic Ocean, and the Nordic Seas.The most distinct feature was thus a thermally contrasting pattern over North America, with simulated temperatures being around 2 • C higher than those the PI for Alaska, whereas over most of Canada temperatures were more than 3 • C lower.

Transient simulation for the Holocene
It is clear in our analysis of Sect.3.1 that the climate showed different responses in the following areas: the Arctic, northwestern Europe, northern Canada, Alaska and Siberia (marked in Fig. S2).Therefore, these areas were selected for special examination and the temperature evolutions of these regions will be shown.Our major focus was on millennialscale temperature trends; therefore, we applied a 500-year running mean to our simulated time series that effectively filtered out high-frequency variability.

Temperature evolution in the Arctic
Arctic summer temperatures in ORBGHG continuously decreased, which resulted in a total cooling of 2 • C during the Holocene (Fig. 4).Winter temperatures showed an even stronger overall cooling of −3 • C.
The simulation OGIS_FWF-v1 with full forcings reveals a more complicated Arctic climate evolution compared to that of ORBGHG.The effect of ice sheets at the onset of the Holocene caused temperatures in both summer and winter to be more than 2 • C lower than those indicated by ORBGHG.The final deglaciation of the FIS happened at by a more gradual warming toward the maximum anomaly of 1 • C warmer than PI at about 7.5 kyr (Fig. 4a).Simulated winter temperatures stayed at a level of 2 • C lower than that of ORBGHG before 7 kyr, which was followed by a rapid increase of about 1.5 • C within a 500-year period, and then reached a temperature peak of about 1.5 • C warmer than the PI (Fig. 4b).Simulated annual mean temperatures showed a relatively stable rise until 6.5 kyr, which reached a maximum of about 1.5 • C warmer than that of PI (Fig. 4c).The simulation OGIS_FWF-v2 gives similar results for the Arctic but had an even cooler climate before 9 kyr than in OGIS_FWF-v1, with the maximum cooling of up to 0.3 • C for all seasons at 10.5 kyr.

Temperature evolution in northwestern Europe
The ORBGHG simulation indicates smaller climate variability in northwestern Europe than in the Arctic.Temperatures declined by around 1.5 • C through the entire period in sum- mer and less than 0.5 • C for annual mean, and rose by 0.5 • C in winter (Fig. 5), which implies a decreasing seasonality toward the preindustrial era.This contrasts markedly with the clear cooling of climate in each season in the Arctic.
The OGIS_FWF-v1 simulation shows an overall cooler climate in northwestern Europe at the onset of Holocene, with temperature anomalies of −1.5 • C in summer, −3 • C in winter and −2.8 • C in annual mean compared to the PI simulation.Temperatures increased from this point (11.5 kyr) toward 6 kyr at an overall rate of 0.28, 0.48, and 0.54 • C kyr −1 for summer, winter and annual mean, respectively.The most important feature in summer was a sharp temperature rise from a negative anomaly (−1.5 • C) to a positive one (+1 • C) by 10 kyr, when the first peak was reached.Subsequently, a slight cooling was noted before 8 kyr, followed by another temperature increase, which led to a second warming peak at 7.4 kyr.The climate in winter showed a relatively stable warming by 6.5 kyr with no identifiable warm peak.Annual temperatures reflected the same phases of warming as in summer, one before 10 kyr and another before 7.5 kyr, but without a clear early temperature peak.Temperatures in all seasons from around 7 kyr followed the ORBGHG simulation.It is worth noting that the OGIS_FWF-v2 simulation indicates a further cooling in summer between 11.5 and 9 kyr compared to OGIS_FWF-v1, which was also reflected in annual mean temperatures.As a result, there was only one  clear thermal maximum in summer for northwestern Europe, which peaked at around 7.4 kyr.

Temperature evolution in northern Canada
Simulated temperatures in ORBGHG for northern Canada decreased by 2.2 • C in summer, by 0.6 • C in winter and by 1.1 • C for annual average during the Holocene (Fig. 6).The stronger cooling in summer than in winter reflected a strong early Holocene seasonality, which decreased over the whole period.
The OGIS_FWF-v1 simulation describes a much cooler climate in northern Canada during the early Holocene than that indicated by ORBGHG.This cooling was up to 5 • C for all seasons at the onset of Holocene.The climate dramatically warmed with an overall high rate of more than 1 • C kyr −1 in both winter and summer during the early Holocene, which was due to the impact of the decaying LIS.The early Holocene warming was, however, not linear because an initial phase with more rapid warming was followed by a more gradual temperature increase.In summer, this warming resulted in a thermal peak at around 7.4 kyr, which was about 1.5 • C warmer than the PI.From 7.4 kyr onwards, the climate experienced a gradual cooling that was very similar to that of ORBGHG.Simulated temperatures in winter and annual mean did not show such a clear warm peak in comparison to summer.The results of OGIS_FWF-v2 only indicated marginal differences relative to OGIS_FWF-v1 for all seasons.Overall, the most significant feature of simulated temperatures in northern Canada was the strong warming that took place in the early Holocene.

Temperature evolution in Alaska
The ORBGHG simulation shows an overall cooling in Alaska for all seasons.Simulated summer and annual mean temperatures experienced a decrease of more than 2 • C throughout the whole period.Winter temperatures had slightly increased by 10 kyr, and then stayed about 2 • C higher for a period of 800 years, which was followed by a constant decrease toward the preindustrial value.
In contrast to other areas, both summer and winter temperatures in OGIS_FWF-v1 showed an overall cooling trend in Alaska during the entire Holocene (Fig. 7), which was slightly higher than in our ORBGHG simulation.The OGIS_FWF-v1 simulation indicates a 2 • C decline in summer temperature over the whole period, with a slightly faster rate between 7 and 6.5 kyr.Simulated winter temperatures decreased by 3.5 • C during the early Holocene, with two small declines at 9.5 and 6.7 kyr.Annual temperatures in the OGIS_FWF-v1 simulation reflected a 2.3 • C cooling during the Holocene.The OGIS_FWF-v2 simulation represents an Alaskan temperature trend that is rather similar to that of OGIS_FWF-v1.

Temperature evolution in Siberia
The ORBGHG simulation describes an almost 2 • C decline of summer temperatures over Siberia during the last 11.5 thousand years (Fig. 8).Simulated winter temperatures showed a smaller variation, as it decreased by less than 1 • C, and annual mean temperatures decreased by around 1 • C over the course of the Holocene.The evolution of simulated temperatures in ORBGHG over Siberia was on a similar scale to that of northwestern Europe.
The difference of simulated Siberian temperatures between ORBGHG and OGIS_FWF-v1 varied in summer and winter.On the one hand, simulated summer temperatures in OGIS_FWF-v1 were generally similar to that in OR-BGHG with the exception of a small warming of 0.7 • C before 10 kyr.On the other hand, winter temperatures in the OGIS_FWF-v1 simulation were around 2 • C lower than in ORBGHG before 7 kyr, followed by a rapid increase over the next 500 years, after which it followed the ORBGHG simulation.Consequently, simulated early Holocene warming lasted much longer in winter than in summer.Simulated Siberian temperature evolution in OGIS_FWF-v2 generally followed that of OGIS_FWF-v1.

Discussion
We will evaluate our results by briefly comparing the simulations with proxy-based reconstructions, which will be followed by an analysis of the mechanism behind the simulated temperature patterns.The impact of freshwater forcing will also be discussed based on the two FWF scenarios.

Comparison of simulations with proxy records
At the onset of the Holocene, the overall cool climate indicated by the reconstructions generally matches that of our OGIS11.5simulation, which shows lower annual temperatures at 11.5 kyr than the PI.Climate reconstructions based on proxy data generally show a cooler early Holocene over northern Europe than at 0 kyr both in the summer and winter (Heiri et al., 2104;Mauri et al., 2015).Terrestrial and ocean sediment data also suggest a cooler early Holocene climate over eastern Siberia (Klemm et al., 2013;Tarasov et al., 2013) and slightly lower SSTs over the North Atlantic Ocean (Came et al., 2007;Berner et al., 2008).Cooler conditions over the Barents Sea and Greenland are also indicated by multiple proxies (Peros et al., 2010;de Vernal et al., 2013;Vinther et al., 2008).Therefore, these proxy data agree with simulated lower temperatures over these areas.
However, there is less agreement with proxies in places where the reconstructions are sparse.The only available pollen-based reconstruction from the western side of the Ural Mountains suggests similar early Holocene summer temperatures (within 1 • C anomaly) compared to the preindustrial era (Salonen et al., 2011), whereas OGIS11.5 indicates that summer temperatures were slightly higher at 11.5 kyr over most areas.At high latitudes, the sea-ice cover reconstructions serve as an indirect paleotemperature proxy due to the scarcity of temperature records, and reveal an inconclusive temperature signal over the Canadian Arctic (de Vernal et al., 2013), whereas our simulation reflects an overall warmer climate in the west and cooler conditions in the east.
Proxies indicate significantly different climate patterns over the east and the west of northern America.The later initiation and termination of HTM over northern Canada imply lower temperatures during the early Holocene in the east (Kaufman et al., 2004).However, the higher-thanpresent early Holocene temperatures over central Beringia and Alaska are reflected by peat accumulation and by northward expansion of animal species (Kaufman et al., 2004;Jones and Yu, 2010).This thermal contrast agrees with those simulated patterns in the OGIS11.5 simulation, which indicates warmer temperatures for Alaska and a much cooler climate over Canada.However, this interpretation of high temperature was recently challenged by Kaufman et al. (2016), who argued that the highest summer temperature in Alaska occurred as late as 8-6 kyr.Hence our simulation agrees better with the interpretation of Kaufman et al. (2004).In general, our simulation with full forcings was able to capture main temperature features indicated in proxy-based reconstructions.Shakun et al. (2012) and Marcott et al. (2013) stacked multiple proxies to construct a record of temperatures since the LGM.Both above stacked reconstructions and our simulation OGIS_FWF-v2 show that the Holocene was generally characterized by an initial warming and subsequent Holocene warm period over the NH extratropics, which indicates the broad consistency between simulation and proxy data.However, there are some disagreements related to seasonality (Fig. 9).Marcott et al. (2013) interpreted the stacked temperature reconstruction as representative of the annual mean climate, whereas it shows a better agreement with our simulated summer temperature than with annual mean value (Fig. 9).One potential explanation for this seasonal mismatch is that some proxy records have seasonal bias toward summer conditions, as has been suggested recently for many marine-based SST reconstructions from high latitudes (Lohmann et al., 2013).Further region-by-region comparisons of these warming rates with proxy records are beyond the focus of this work and will be dealt with in a future publication.

Mechanism of climate response to forcings
It is clear from our data that the spatial patterns of climate response at the onset of the Holocene can be attributed to the variation in the dominant forcings prevailing in the different areas.Orbital-scale insolation variations are important driving factors for the early Holocene climate.For instance, higher temperatures in Alaska could be attributed to the orbitally induced positive insolation anomaly in combination with an anomalous atmospheric circulation caused by the remnant LIS.The air descended over the cold LIS surface, which created a high surface pressure anomaly that produced a clockwise flow anomaly at the surface, as indicated by the 800 hPa geopotential height (Fig. 10).This induces stronger southerly winds over Alaska, which advected relatively warm air from the south.A potentially different early Holocene atmospheric circulation near the North Atlantic was also found in a proxy record of Steffensen et al. (2008), who reported an abrupt transition of deuterium excess that indicates a temperature change of precipitation moisture sources and is thus indirectly connected to atmospheric circulation changes.
The strong influence of the ice sheets on early Holocene temperatures has been found in previous studies (Renssen et al., 2009(Renssen et al., , 2012)).Simulated lower summer temperatures over northern Canada and northwestern Europe in our OGIS11.5simulation were the result of such ice-sheet-induced cooling, which would have fully overwhelmed the warming effect of the positive summer insolation anomaly.The ice-sheet cooling effect could partly be explained on a local scale by the enhanced albedo over the ice sheets and by the climate's high sensitivity to albedo change (Romanova et al., 2006).Indeed, the summer surface albedo over the ice sheets was much higher (up to 0.8) than over ice-free surfaces, where the values varied from only 0.1 to 0.5, depending on the vegetation type and the fractional snow cover (Fig. 11).Temperatures could be further reduced by the ice-sheet orography impact.The elevation of ice sheets introduced descending air over the ice-sheet surface, which caused locally cooler conditions.There was also an approximate 0.5 • C cold bias induced by the lapse rate effect when compared with the sitebased records.
Changes in vegetation and land cover during the early Holocene contributed to climate change as well, especially over ecotonal regions.Modeling studies suggest that deforestation in boreal regions could decrease regional temperatures by up to 1 • C due to an increase in surface albedo and related positive feedbacks (Levis et al., 1999;Claussen et al., 2001;Liu et al., 2006).Taking Siberia as an example, the insolation-induced warming was partially offset by the overall higher summer albedo (Fig. 11) induced by the southward expansion of the tundra and/or bare ground and related feedbacks at 11.5 kyr, resulting in a minor warming in summer.The albedo-related feedbacks and the smaller annual insolation anomalies jointly result in a 0.5-2 • C cooler in annual climate at 11.5 kyr.We are aware of the potential role of permafrost at high latitudes; however, the discussion of the impact of permafrost thaw is hindered by the fact that our model version did not include a dynamic permafrost module.A version of LOVECLIM that is coupled to a permafrost module (VAMPERS) is currently in development (Kitover et al., 2015), and should enable us to quantify the role of permafrost in a future study.
Meltwater release and sea-ice-related changes also had a footprint in the early Holocene climate.The OGIS11.5 simulation produces a sluggish AMOC in the North Atlantic  with the largest decrease being more than 3 Sv.It was also reflected in a shallower overturning circulation at 11.5 kyr compared to the PI simulation as a response to meltwater release (Fig. 12).This slowdown also coincides with the foraminifera data from the Arctic Ocean and the Fram Strait that suggest a reduced northward oceanic heat transport (Thornalley et al., 2009).The slowdown and reduced heat transport led to slightly lower temperatures at high latitudes (western Arctic Ocean) at 11.5 kyr than that at 0 kyr.Likewise, after the meltwater fluxes of the LIS diminished around 7 kyr, strong intensification of the AMOC followed.This sudden intensification of AMOC would explain the rapid Arctic temperature increase that occurred at this time (Fig. 4).However, it is important to note that the temperature decrease was not simply inversely linear with the amount of northward transport of heat since the sea-ice feedbacks further reinforce this change (Roche et al., 2007).In fact, sea-ice coverage in the OGIS11.5 simulation was much more extensive over the Davis Strait (northern Labrador Sea) than the corresponding value in OG11.5 (Fig. 13).This extended seaice cover in this region was stronger than the direct cooling effect of the reduced oceanic heat transport.Such an anomaly might be explained by positive feedbacks involving sea ice being active (Renssen et al., 2005).The Greenland Sea warming could be attributed to enhanced convective activity that releases more oceanic heat into the atmosphere.This enhanced convective activity was caused by the shift of deep water formation from the eastern Greenland Sea to the west, which was initially induced by the freshwater discharge from ice-sheet melting.The net response of the climate reflects the impact of a combination of forcings and feedbacks, which showed a high temporal-spatial variability.

Early Holocene warming and climate-ocean system response to freshwater
The simulation of OGIS_FWF-v2 indicates a stronger cooling (before 9 kyr) in the Arctic and northwestern Europe than those found in the OGIS_FWF-v1 with the strongest temperatures reduction at around 10 kyr.The enhanced freshwater influx from the GIS and the redistributed meltwater from the FIS caused an alteration in the surface ocean freshening in the Nordic Seas, which reduced convective activity (Renssen et al., 2010;Blaschek and Renssen, 2013).Indeed, this reduction led to a further slight reduction of the northward heat transport by the Atlantic Ocean, which was associated with a further AMOC weakening (by 1-2 Sv than in FWF_v1): this in turn produced a slightly stronger cooling  at 10 kyr (Fig. S3) and a sea-ice expansion over the Denmark Strait (Figs. 14 and S4).However, the efficiency of the above meltwater flux freshening effect is determined by multiple aspects.The most important factor is the maximum flux of meltwater that was added to the ocean, while the total freshwater amount had only a second-order effect (Roche et al., 2007).Numerous investigations on the behavior of the coupled atmosphere-ocean system suggest that the application of freshwater will not lead to a disruption of the North Atlantic Deep Water production (NADW) as long as a certain threshold is not crossed (Ganopolski et al., 1998;Rahmstorf et al., 2005).Apart from the intensity and duration, the ocean circulation response to freshwater also depends on the location where this freshwater is released.For instance, it is more sensitive to the release of freshwater in the eastern Norwegian Sea than at the St. Lawrence River outlet since the former is closer to the main site with NADW formation (Roche et al., 2010).This is consistent with a previous study by Blaschek and Renssen (2103), who found that freshwater from the GIS did have a tangible impact on the Nordic Seas, even though the total amount was minor.
Since the second freshwater scenario (OGIS_FWF-v2) includes a slightly larger FWF from the GIS (compared to that in OGIS_FWF-v1) and the FWF was released in a sensitive area, the location-dependent sensitivity could also partially explain further AMOC weakening in the OGIS_FWF_v2 simulation compared to OGIS_FWF-v1.
The OGIS_FWF-v1 simulation indicates two peaks in the temperature evolution over northwestern Europe, at around 10 and 7 kyr.High temperatures at 7 kyr are recorded in proxy-based reconstructions as well.However, no warm peak at 10 kyr was observed in pollen-based reconstructions, which actually suggests a cooler climate prevailed at 10 kyr than in the preindustrial Europe (Mauri et al., 2015).In contrast to the climate simulated in OGIS_FWF-v1, the simulation with updated freshwater (OGIS_FWF-v2) produced a warming trend that is consistent with a highest temperature around 7 kyr.Moreover, the OGIS_FWF-v1 produced a temperature decrease between two peaks, whereas the proxies indicated a rapid temperature increase at the beginning followed by a more gradual warming (Brooks et al., 2012  northwestern Europe, the OGIS_FWF-v2 represented a more realistic climate than OGIS_FWF-v1 did, which implies that the existing uncertainties in the reconstructions of ice-sheet dynamics can be evaluated by applying different freshwater scenarios.Further comparison with proxy data and with other model transient simulations will be conducted in a future paper.

Conclusions
We performed both equilibrium and transient simulations by employing the LOVECLIM climate model to explore the spatial patterns of the climate response to forcings at the onset of the Holocene and temperature evolution over the last 11.5 thousand years.We focused on three research questions in our analysis, which are outlined below with the main finding: 1. What were the spatial patterns of simulated temperature at the onset of the Holocene?
The temperature anomalies relative to PI at 11.5 kyr were regionally heterogeneous, which are shown as a range of annually negative anomalies over many areas but which were positive in Alaska.The climate in eastern northern Canada and northwestern Europe was much cooler than in other regions, with temperature anomalies of −2 to −5 • C relative to 0 kyr throughout the year.The climate over the northern Labrador Sea and the North Atlantic was also 0.5-3 • C cooler.Temperatures in Siberia were 0.5-3 • C and 1.5-3 • C lower in winter and annually, respectively, and summer temperatures showed only a small deviation (between +0.5 and −1.5 • C) compared to 0 kyr.Simulated summer temperature anomaly in the eastern Arctic Ocean was also small (between +0.5 and −0.5 • C), and annual temperatures were 0.5-2 • C lower.In contrast to cooler conditions in other areas, temperatures in Alaska were 1.5-3 • C higher than the preindustrial period for all seasons.
2. What were the roles of forcings, especially ice-sheet decay, in shaping these features?
The ice-sheet cooling effect in northern Canada and northwestern Europe overwhelmed the warming impact of the positive insolation anomaly, which caused the relatively cold climate at 11.5 kyr.In particular, the enhanced surface albedo over the ice sheets and the orographic effect were important in promoting these cold conditions.The cooler climate over the northern Labrador Sea and the North Atlantic was related to both reduced northward heat transport and enhanced sea-ice feedbacks.A small summer temperature anomaly was found in Siberia, where the positive insolation anomaly was partially offset by the cooling effect of the higher albedo associated with the relatively extensive tundra cover in the early Holocene.Overall, lower winter and annual temperatures at 11.5 kyr over central Siberia can be attributed to both vegetation-related albedo feedbacks and to the relatively small negative insolation deviation compared to the preindustrial level.
The dominant factors driving the climate in eastern Arctic Ocean climate were the amount of northward heat transport associated with the strength of ocean circulation and the orbitally forced insolation variation.Annual mean temperatures at 11.5 kyr were lower than at 0 kyr because the cooling effect of a reduced northward oceanic heat transport (induced by weakened ocean circulation) was larger than the insolation-induced warming.During summer, these two factors were of similar magnitude and temperatures were similar to those of the preindustrial era.Temperatures in Alaska were higher for all seasons in response to the dominant positive insolation anomaly and the enhanced southerly winds induced by the LIS, which advected relatively warm air from the south.Therefore, this regional heterogeneity is the result of the climate response to a range of dominant forcings and feedbacks.
3. What was the spatiotemporal variability in the simulated early Holocene evolution?
Above, geographical variability is also reflected in the Holocene temperature evolution, especially in the early Holocene warming.In Alaska, the climate was constantly cooling throughout the Holocene due to the decreasing insolation and atmospheric circulation variability.In contrast, northern Canada experienced a strong warming with an overall warming rate over 1 • C kyr −1 , and this warming lasted until 7 kyr.Although different forcings and mechanisms played different roles in northwestern Europe, the Arctic and Siberia, the overall warming effect was similar for these regions, with a rate of around 0.5 • C kyr −1 .In addition, the comparison of early Holocene temperatures over northwestern Europe with proxy records suggests that the OGIS_FWF-v2 represented a more realistic climate condition than the OGIS_FWF-v1 does, and it implies that the uncertainties with regard to the ice-sheet decay can potentially be constrained by applying different deglaciation scenarios and comparing then with networks of proxy records.Overall, our results demonstrated a large spatial variability in the climate response to diverse forcings and feedbacks, both for the early Holocene temperature distribution and for the early Holocene warming, and this data-model comparison also helps in understanding the difference between proxy records.
The Supplement related to this article is available online at doi:10.5194/cp-12-1119-2016-supplement.

Figure 1 .
Figure1.Evolution of greenhouse gas concentrations (GHGs) shown as the radiative forcing's deviation from the preindustrial level (with solid lines corresponding to the left axis), and June insolation at 65 • N derived from orbital configuration (with red line and the axis on the right).

Figure 2 .
Figure 2. The prescribed ice-sheet forcing during the early Holocene.(a) Variation in ice-sheet extent (km 2 ) displayed as the black lines with the axis on the left and their maximum thickness (m) indicated by the green lines with the axis on the right.A relatively minor change in GIS is not shown due to its small scale.(b) Two freshwater flux scenarios (in mSv), FWF-v1 (thick dashed lines) and FWF-v2 (solid lines).(c) Total meltwater discharge in equivalent sea level (m).

Figure 3 .
Figure 3. Simulated temperatures for 11.5 kyr, shown as deviation from the PI.Left column shows the simulation with only GHG and ORB forcings (OG11.5).For the right column, the ice-sheet forcing is included (OGIS11.5).Upper, middle and lower panels present summer (JJA), winter (DJF) and annual mean temperatures, respectively.

Figure 4 .
Figure 4. Simulated temperature evolution, shown as the anomalies compared to PI, since the early Holocene at high latitudes (north of 70 • ).Panels (a), (b) and (c) represent the summer, winter and annual values, respectively.The slope indicates the overall warming rate and is based on the least-squares regression over the period from the 11.5 to 6 kyr, as from 6 kyr the temperature starts to decrease.It is only a general estimation, and thus uncertainty ranges are not provided.The warmest peak is marked by a shaded bar and represents the simulated peak during which the temperature was over 1 • C higher than the PI.Both slope calculation and warm peak are based on the OGIS_FWF-2 simulation.

Figure 9 .
Figure 9. Model-data comparison over the latitudinal band of 30-90 • N, shown as a deviation from the PI.The stacked temperature reconstruction with 1δ uncertainty (grey band) is based on Marcott et al. (2013).

Figure 10 .
Figure 10.Geopotential height anomalies from global mean (in m 2 s −2 ) at 800 hPa in the extratropical Northern Hemisphere.Panels (a), (b) and (c) show the control condition PI and the simulations OG11.5 and OGIS11.5, respectively.

Figure 12 .
Figure 12.Meridional overturning stream function (Sv) in the Atlantic Ocean basin.Panels (a), (b) and (c) indicate the control run (PI) and the simulations OG11.5 and OGIS11.5, respectively.On the left-hand side, depth is indicated in kilometers.Positive values indicate a clockwise circulation.Maximum AMOC strength value was 22 Sv (reached at about 1200 m depth) in the PI and OGIS simulation, while it was only about 14 Sv (reached at 600-700 m depth) in OGIS11.5.

Figure 14 .
Figure 14.Response of the ocean variables (shown as a 100-year average) to forcings during the Holocene.(a) Maximum meridional overturning stream function (Sv) in the North Atlantic.(b) Sea-ice area (10 12 m 2 ) in the Northern Hemisphere.

Table 1 .
Boundary conditions for 11.5 kyr and the preindustrial (PI) era